11institutetext:  Laboratoire d’astrophysique de Bordeaux, Univ. Bordeaux,
CNRS, B18N, allée Geoffroy Saint-Hilaire, 33615 Pessac, France
22institutetext: Laboratoire de Météorologie Dynamique (IPSL), Sorbonne
Université, Centre National de la Recherche Scientifique, École Polytechnique,
École Normale Supérieure, Paris, France 33institutetext: LESIA, Observatoire
de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, 5
place Jules Janssen, 92195 Meudon, France 44institutetext: LATMOS/IPSL,
Sorbonne Université, UVSQ Université Paris-Saclay, CNRS, Paris, France

# A 3D picture of moist-convection inhibition
in hydrogen-rich atmospheres: Implications for K2-18 b

Jérémy Leconte 11 Aymeric Spiga 22 Noé Clément 11 Sandrine Guerlet 2233 Franck
Selsis 11 Gwenaël Milcareck 22 Thibault Cavalié 1133 Raphaël Moreno 33
Emmanuel Lellouch 33 Óscar Carrión-González 33 Benjamin Charnay 33 Maxence
Lefèvre 44

(Submitted 14 December 2023 / Accepted 9 January 2024)

While small, Neptune-like planets are among the most abundant exoplanets, our
understanding of their atmospheric structure and dynamics remains sparse. In
particular, many unknowns remain on the way moist convection works in these
atmospheres where condensable species are heavier than the non-condensable
background gas. While it has been predicted that moist convection could shut-
down above some threshold abundance of these condensable species, this
prediction is based on simple linear analysis and relies on some strong
assumptions on the saturation of the atmosphere. To investigate this issue, we
develop a 3D cloud resolving model for hydrogen-dominated atmospheres with
large amounts of condensable species and apply this model to a prototypical
temperate Neptune-like planet — K2-18 b. Our model confirms the shut-down of
moist convection above a critical abundance of condensable vapor and the onset
of a stably stratified layer in the atmosphere of such planets, leading to
much hotter deep atmospheres and interiors. Our 3D simulations further provide
quantitative estimates of the turbulent mixing in this stable layer, which is
a key driver of the cycling of condensables in the atmosphere. This allows us
to build a very simple, yet realistic 1D model that captures the most salient
features of the structure of Neptune-like atmospheres. Our qualitative
findings on the behavior of moist convection in hydrogen atmospheres go beyond
temperate planets and should also apply to the regions where iron and
silicates condense in the deep interior of hydrogen-dominated planets.
Finally, we use our model to investigate the likelihood of a liquid ocean
beneath a H2 dominated atmosphere on K2-18 b. We find that the planet would
need to have a very high albedo ($A>0.5-0.6$) to sustain a liquid ocean.
However, due to the spectral type of the star, the amount of aerosol
scattering that would be needed to provide such a high albedo is inconsistent
with the latest observational data.

††offprints<EMAIL_ADDRESS>
##  1 Introduction

Convection, the process by which an unstably stratified fluid transports
energy upward to restore neutrality, shapes the thermal structure of all the
deep planetary atmospheres in our Solar System. Inherently caused by the
difference in the radiative opacity of the atmosphere between the wavelengths
at which it is heated by the star and those at which it can cool by thermal
radiation, convection usually happens wherever radiative processes are
insufficient to carry energy out. This leads to atmospheres consisting of two
main layers: a deep troposphere overlain by a stratosphere. This structure
stems from such basic principles that it is envisioned to hold true on most
exoplanets with a substantial atmosphere (Robinson & Catling, 2014), although
the intense irradiation they receive can mean that the stable radiative zone
can extend very deep (Guillot & Showman, 2002).

But for this balance to hold, convection itself needs to be able to develop,
for it can be hindered by other dynamical processes. For example, when there
is a compositional gradient in the atmosphere - whether it is caused by
condensation or chemical reactions - the resulting mean molecular weight
gradient in the gas can affect the thermal gradient needed to initiate
convection (Nakajima et al., 2000; Garaud, 2018; Daley-Yates et al., 2021;
Habib & Pierrehumbert, 2023). In a more extreme fashion, when the atmosphere
contains enough of a condensable species (hereafter referred to as vapor) that
is heavier than the non-condensable background gas (hereafter referred to as
air), condensation can completely suppress convection, whatever the thermal
gradient (Guillot, 1995; Leconte et al., 2017; Markham et al., 2022).
Intuitively, this is due to the fact that when we consider saturated moist air
following the Clausius-Clapeyron law, there is a threshold specific
concentration of condensable vapor,

$\displaystyle
q_{\mathrm{cri}}\equiv\frac{1}{M_{v}-M_{d}}\frac{R^{\star}T}{L},$ (1)

above which a change in vapor abundance due to condensation affects buoyancy
more (and in the opposite direction) than the temperature change that led to
it. In this equation, $M_{v}$ and $M_{d}$ are the molar masses for the vapor
and dry air respectively, $L$ is the specific latent heat, $T$ the
temperature, and $R^{\star}$ the universal molar gas constant. Above that
threshold, the effective thermal expansion coefficient of the fluid becomes
negative (Markham et al., 2022). A super-adiabatic region thus becomes
perfectly stable to convection because parcels of fluid become denser than
their environment as they rise.

Leconte et al. (2017) hypothesized that in hydrogen-dominated atmospheres, a
stable layer would form near the cloud deck of any species for which the deep
abundance would be higher than the critical concentration $q_{\mathrm{cri}}$.
This critical concentration is probably reached for water in Saturn (Li &
Ingersoll, 2015), and for water and methane in Uranus, Neptune, and most exo-
Neptunes. Markham et al. (2022) and Misener & Schlichting (2022) independently
extended this argument to iron at the core-envelope boundary in small neptune-
like exoplanets. The main consequence of this stable, super-adiabatic layer is
that the deep atmosphere and interior can be much hotter than predicted by
standard models for the same effective temperature (hence the same temperature
of the upper troposphere and stratosphere), with implications for the past
evolution and chemistry of our ice giants (Cavalié et al., 2017, 2020; Markham
& Stevenson, 2021).

However, the process of convection inhibition has been mostly studied
analytically with unidimensional arguments, and many caveats remain:

  * $\bullet$

The analytical theory of Leconte et al. (2017) assumed a fully-saturated
medium, whereas we know that on Earth, moist convective regions are all but
fully-saturated. Would a stable layer really form in a realistic, subsaturated
atmosphere?

  * $\bullet$

If the stable layer is indeed devoid of any turbulence, as has been predicted,
it should prevent any upward transport of vapor, entailing a dry upper
troposphere without moist convection. However, this contrasts with the
observation of moist convective systems in the atmosphere of ice giants. Why
is that?

  * $\bullet$

How does convection inhibition affect the atmosphere of neptune-like
exoplanets and is this observable with our current methods and observatories?

To answer these questions, we develop a 3D cloud-resolving model that is able
to simulate hydrogen-dominated atmospheres with large amounts of any
condensable volatile (see Sect. 2). Because the atmosphere of Solar System ice
giants have radiative timescales on the order of decades to centuries, we
focus on simulating the atmosphere of a prototypical temperate sub-neptune:
K2-18b (Cloutier et al., 2017). While enabling the equilibration of the
simulation with reasonable computing resources, this also helps us shed light
on the atmospheric structure of such temperate, hydrogen-dominated planets
(Sect. 3). We argue that many features exhibited by convection in this regime
should be fairly general and apply to other systems with moist-convection
inhibition. In particular, we will show that, while being more realistic, the
3D simulations largely support the thermal structure envisioned in the earlier
study of Leconte et al. (2017). In addition, a detailed study of the rich
energy and volatile cycles in the 3D simulations allows us to build a much
more consistent picture of these atmospheres and to develop a brand new 1D
model that can be used to explore their diversity (Sect. 4). Finally, in Sect.
5, we discuss how the chemical composition of temperate Neptune-like planets
is affected by convection inhibition, potentially acting as a tracer of this
process.

Interestingly, during the writing of this manuscript, new near- and mid-
infrared JWST transit observations of K2-18 b have been released (Madhusudhan
et al., 2023). One of the most intriguing feature of these observations is the
potential non-detection of NH3 that has been interpreted as the sign of a
shallow atmosphere overlying a liquid ocean. Indeed, if a deep atmosphere were
present, the thermochemical conditions at the bottom of such an atmosphere
should replenish the upper atmosphere in ammonia (Hu et al., 2021). Although
this non-detection will probably need to be confirmed when our knowledge of
instrumental systematics improve, we have tried to quantify in more details
whether or not the presence of a liquid ocean can be consistent with both
these new observations and our improved knowledge of the thermal structure of
these temperate sub-neptunes. For that purpose, in Sect. 6, we use our updated
1D model to put limits on the conditions necessary to sustain a liquid ocean
on K2-18 b and show that this requires planetary Bond albedos in excess of
0.5-0.6, that are relatively hard to achieve around a late-type star. We
further demonstrate that the kind of aerosol properties that we need to create
such a high albedo are incompatible with the current transit observations that
exhibit relatively deep near-infrared methane bands.

##  2 3D hydrodynamical cloud-resolving model

The mechanism for moist-convection inhibition rests on a fairly basic
principle: when the mean molecular weight of the vapor is larger than the one
of the air, a parcel of fluid becomes denser as its vapor concentration
increases. The ingredients for convection inhibition are therefore inherently
present in the basic equations of the hydrodynamics of a moist atmosphere.
However, in many terrestrial atmospheres, the condensable species are found in
trace amounts and neglecting their contribution to the density of the air and
to the global mass of the atmosphere is a rather valid assumption, used in
many numerical models. Such models thus cannot exhibit convection inhibition.

Fortunately, the ”trace gas” assumption is challenged when modeling strong
convective events in hot and moist regions of the Earth—even if on Earth, the
presence of water vapor facilitates convection as it is lighter than molecular
nitrogen. For this reason, many small-scale, non-hydrostatic models of the
atmosphere – the so-called cloud-resolving models – incorporate the mass-
loading effect of water vapor (Bryan & Fritsch, 2002; Skamarock et al., 2019).

To model the convection in hydrogen-rich atmospheres, we thus used the
dynamical core of the 4th version of the Weather Research and Forecast model
(hereafter WRF V4) described in Skamarock et al. (2019) that we coupled with
the physical parametrizations of the LMD Generic Planetary Climate Model
(Generic PCM, the former Generic LMD GCM; Wordsworth et al. 2011; Leconte et
al. 2013). The following sections describe the details of each code and their
coupling as well as the specific developments and the numerical setup used.

###  2.1 Equations

We simulate the compressible, non-hydrostatic Euler equations using the
philosophy described in Laprise (1992) that uses a terrain-following
hydrostatic-pressure vertical coordinate defined as:

$\displaystyle\eta_{d}\equiv\frac{\pi_{d}-\pi^{t}}{\pi_{d}^{s}-\pi^{t}}\equiv\frac{\pi_{d}-\pi^{t}}{\mu_{d}},$
(2)

where $\pi^{t}$ and $\pi_{d}^{s}$ are the dry hydrostatic pressure at the top
of the model and at the surface (it is assumed that there is negligible
moisture above the model top), and dry mass in the model column is
$\mu_{d}\equiv\pi_{d}^{s}-\pi^{t}$.

The hydrostatic pressure, $\pi$, is the local pressure the fluid would have if
in hydrostatic equilibrium with the same mass above the point considered. For
a constant gravity, it is equal to the mass of fluid above a given point
divided by gravity. Defining the geopotential, $\phi\equiv gz$, and the
specific volume of the fluid, $\alpha$, one can write the equivalent of the
hydrostatic equation

$\displaystyle\frac{\partial\phi}{\partial\pi}=-\alpha,$ (3)

except this equation is valid even when there is no equilibrium. When there is
moisture in the atmosphere, we define $\pi_{d}$ as the mass of dry air above a
given point divided by $g$. It can be easily shown that

$\displaystyle\frac{\partial\pi}{\partial\pi_{d}}=\frac{\alpha_{d}}{\alpha},$
(4)

so that

$\displaystyle\frac{\partial\phi}{\partial\pi_{d}}=-\alpha_{d}.$ (5)

The mass-weighted, prognostic atmospheric variables in the flux-form dynamical
equations are

$\displaystyle\mathbf{V}\equiv(U,V,W)\equiv\mu_{d}\mathbf{v},\
\Omega\equiv\mu_{d}\dot{\eta}_{d},\ \Theta\equiv\mu_{d}\theta,\
\chi_{\mathrm{i}}\equiv\mu_{d}r_{\mathrm{i}},$ (6)

and the geopotential $\phi$ that is not written in flux form as it is not a
conserved quantity. The usual velocities are denoted $\mathbf{v}=(u,v,w)$ and
we define the horizontal velocity for sake of compacity
$\mathbf{\bar{V}}\equiv\mu_{d}\bar{\mathbf{v}}\equiv\mu_{d}(u,v)$.
$\theta\equiv T(p_{0}/p)^{R/c_{p}}$ is the potential temperature, $p_{0}$ is
an arbitrary reference pressure, $R$ the specific constant of the gas, and
$r_{\mathrm{i}}$ is the mass mixing ratio of the various tracer species with
respect to dry air (in particular water vapor with index $v$ and condensed
water with index $c$). For the sake of compactness of notation, the specific
concentration of the tracers, $q_{\mathrm{i}}=r_{\mathrm{i}}/(1+r_{v})$, will
sometimes be used when appropriate.

Adapting the equations of Kasahara (1974) and Laprise (1992) to a moist
atmosphere yields the following prognostic equations in flux form:

$\displaystyle\partial_{t}U+(\nabla\cdot\mathbf{V}u)-\mu_{d}\alpha\
\partial_{x}p+(\alpha/\alpha_{d})\ \partial_{\eta_{d}}p\ \partial_{x}\phi$
$\displaystyle=\mu_{d}F_{x},$ (7)
$\displaystyle\partial_{t}V+(\nabla\cdot\mathbf{V}v)-\mu_{d}\alpha\
\partial_{y}p+(\alpha/\alpha_{d})\ \partial_{\eta_{d}}p\ \partial_{x}\phi$
$\displaystyle=\mu_{d}F_{y},$ (8)
$\displaystyle\partial_{t}W+(\nabla\cdot\mathbf{V}w)-g\left[(\alpha/\alpha_{d})\
\partial_{\eta_{d}}p-\mu_{d}\right]$ $\displaystyle=\mu_{d}F_{z},$ (9)
$\displaystyle\partial_{t}\mu_{d}+\left(\nabla\cdot\mathbf{V}\right)$
$\displaystyle=0,$ (10)
$\displaystyle\partial_{t}\phi+\frac{1}{\mu_{d}}\left[\left(\mathbf{V}\cdot\nabla\phi\right)-gW\right]$
$\displaystyle=0,$ (11)
$\displaystyle\partial_{t}\Theta+\left(\mathbf{V}\cdot\nabla\theta\right)$
$\displaystyle=\mu_{d}\frac{\theta}{T}\frac{Q}{c_{p}},$ (12)
$\displaystyle\partial_{t}\chi_{\mathrm{i}}+\left(\mathbf{V}\cdot\nabla
r_{\mathrm{i}}\right)$ $\displaystyle=\mu_{d}S_{\mathrm{i}},$ (13)

where the various differential operators are defined as follows

$\displaystyle\left(\nabla\cdot\mathbf{V}\right)\equiv\partial_{x}U+\partial_{y}V+\partial_{\eta_{d}}\Omega,$
(14)
$\displaystyle\left(\nabla\cdot\mathbf{V}a\right)\equiv\partial_{x}Ua+\partial_{y}Va+\partial_{\eta_{d}}\Omega
a,$ (15) $\displaystyle\left(\mathbf{V}\cdot\nabla a\right)\equiv
U\partial_{x}a+V\partial_{y}a+\Omega\partial_{\eta_{d}}a.$ (16)

The remaining quantities are $(F_{x},F_{y},F_{z})$ the external forces per
unit mass exerted on the fluid, $Q$ the specific, diabatic heating rate, and
$S_{\mathrm{i}}$ the specific source/sink rates for tracer $\mathrm{i}$ (in
our case, the vapor and condensed phases).

Our set of equations differ from either Daley-Yates et al. (2021) or Habib &
Pierrehumbert (2023) who chose to solve an equation for the conservation of
energy directly, whereas we chose to formulate the problem with a potential
temperature. The effect of this numerical choice should probably be
investigated in more depth by simulating the same setup with these different
codes.

To close the system, we provide an equation of state for the ideal mixture of
vapor and dry gas that writes

$\displaystyle
p=p_{0}\left(R_{d}\theta\left[1+(R_{v}/R_{d})r_{v}\right]/p_{0}\alpha_{d}\right)^{\gamma},$
(17)

where $R_{v}$ and $R_{d}$ are the specific gas constant of the vapor and dry
air, and $\gamma=1/(1-R/c_{p})$ is kept constant. Notice that the ratio
$R/c_{p}$ for the whole gas needs to be constant to express the entropy
equation in the form of an equation on a potential temperature (Eq. (12)).
However, because $R=q_{d}R_{d}+q_{v}R_{v}$ will change with the vapor
concentration, the $c_{p}$ used in the right-hand side of Eq. (12) must be
varied accordingly, i.e. following $c_{p}=q_{d}c_{p,d}+q_{v}c_{p,v}$. This
corresponds to a situation where the vapor and the dry gas would have the same
molar heat capacity.

Apart from its effect on the equation of state, the vapor and condensates
affect the dynamics through the $\alpha_{d}/\alpha\equiv
1+\sum_{\mathrm{i}}r_{\mathrm{i}}$ terms that account for the mass loading
effect of both phases.

###  2.2 Dynamical core and numerical implementation

To solve these equations, we use the 4th version of the Weather Research and
Forecast model (hereafter WRF V4) described in Skamarock et al. (2019). We
refer the reader to this technical note for all the details of the general
numerical implementation and of the various schemes and options mentioned
hereafter.

Our baseline simulations are performed on a 64x64x256 grid with a 2 km
horizontal resolution. The physical vertical resolution is variable as the
grid is based on fixed $\eta_{d}$ levels but varies around 400 m in most of
the domain with thinner layers near the surface. Despite the very large scale
height of the modeled atmospheres compared to Earth, test simulations with
coarser vertical grids seem to indicate that such thin layers are necessary to
satisfactorily resolve the convection. The horizontal boundaries are periodic,
the surface is a rigid lid and the top of the model is a fixed pressure
boundary around 3000 Pa. A damping layer extends over the top 20 km to damp
upward propagating gravity waves following the vertical-velocity implicit
Rayleigh damping scheme.

The dynamical timestep was fixed for each simulation to fulfill the CFL
stability condition and was on the order of 5 s. Physical parametrization
described in the next section are called every 75 s and radiative transfer
calculations are performed once every ten minutes.

As one of the goals of the study is to identify the dynamical processes that
can transport energy and tracers in the modeled atmospheres, we turn off all
the parametrizations of subgrid-scale turbulent transport usually included in
cloud-resolving models. This enables us to know that any transport observed in
the simulations is the result of resolved dynamical motions and not due to ad-
hoc parametrizations (Parmentier et al., 2013).

The various tracers are transported by the dynamics with the 5th order
monotonic advection scheme recently implemented in WRF V4. This proved crucial
in providing physical results in the presence of a sharp vapor gradient
arising in the stable layer above convective motions. In particular, using
only the positive definite scheme (as implemented in WRF V3) creates a
spurious local minimum of vapor concentration above the stable layer that
dried the troposphere above, suppressing moist convection there. This is in
accordance with the known tendency of non-monotonic advection schemes to
create spurious minima around sharp edges (Skamarock, 2006), but here the
manifestation goes well beyond a small inaccuracy, warranting the use of a
monotonic scheme.

###  2.3 Physical parameterisations

The various physical source terms in the equations of motion, in particular
the diabatic heating and the tracer source terms, $Q$ and $S_{\mathrm{i}}$
respectively, are computed in each column of the model using the physical
parametrizations from the Generic PCM. This strategy of coupling the WRF
dynamical core with physical parameterizations for planetary atmospheres has
been developed for martian applications in Spiga & Forget (2009) and was later
extended to Venus (Lefèvre et al., 2018) and temperate exoplanets (Lefèvre et
al., 2021).

In the present study, the coupling has been improved in several ways:

  * •

We transitioned toward the 4th version of WRF to take advantage of the new
features of the dynamical core. The inclusion of a monotonic advection scheme,
in particular, was instrumental in properly modeling the steep gradients in
the stable layer.

  * •

As we model water vapor rich atmospheres, a particular attention has been
devoted to improving the water conservation in the interface between the
dynamical core that uses tracer mixing ratios ($r_{\mathrm{i}}$) and the
physical parametrizations in PCM that use tracer specific concentration with
respect to the whole gas ($q_{\mathrm{i}}=r_{\mathrm{i}}/(1+r_{v})$) – the
difference between those two variables, which was neglected in previous
versions, is no longer negligible in the case studied here.

  * •

As the usual PCM global climate model assumes hydrostatic equilibrium, the
total thermodynamic pressure at level interfaces was used to calculate the
mass of gas in model layers. This led to a non-conservation of the total mass
of the atmosphere. The hydrostatic pressure ($\pi$) provided by the WRF
dynamical core is now used instead.

In terms of physical processes, many parametrizations in the physical part of
the Generic PCM have been developed to be used with a hydrostatic, global
dynamical core: dry and moist convective adjustment, subgrid-scale humidity
distribution, etc. These are turned-off when running the Generic PCM in a
cloud-resolving mode because these processes are resolved in the simulations.
The main remaining parametrizations are the radiative transfer and the
formation and evaporation of precipitations.

For the radiative transfer, use the generic two-stream module of the PCM with
correlated-$k$ coefficient tables. For these specific simulations, we used
Exo_k (Leconte, 2021) to bin-down and combine the correlated-$k$ table
opacities from Blain et al. (2021) into tables with 21 channels for the
thermal emission and 9 for the stellar radiation with 25 pressure points
uniformly distributed in log $P$ between 0.1 and 107 Pa, 15 temperature
points, and 9 values of the water vapor mixing ratio that is variable in the
simulations. The abundances of the other species depend on the chosen
metallicity and will be detailed hereafter. The CIA for H2-H2, H2-He are taken
from the HITRAN database (Richard et al., 2011). The H2O-H2O and H2O-air
continua are based on MT-CKD 3.3 (Mlawer et al., 2012). The stellar spectrum
is approximated by a blackbody at the effective temperature of the star and
the domain receives a constant insolation with an effective solar zenith angle
of 60∘.

For the precipitations, we have opted for a very simple scheme to isolate the
contribution of the dynamics and be able to later identify how more realistic
microphysics impact the atmospheric structure. Parcels of air that are
supersaturated are instantaneously brought back to liquid-vapor equilibrium
iteratively to account for the latent heat effect. Condensates are assumed to
precipitate whenever their mixing ratio exceeds a threshold that we keep
arbitrarily small in this first study to avoid the radiative feedback of
clouds. Precipitations fall instantaneously. Because the planets we model are
envisioned to have a thick atmosphere that reaches deeper that what we can
model, depositing the precipitations at the model surface would be rather
unphysical. Consistently with our choice of a simple microphysics, we want to
keep a physically motivated model without any free parameter. We thus assume
that evaporation is inefficient until droplets reach the boiling level where
bubbles of vapor would form inside the droplets. It is indeed improbable that
rain drops would fall much deeper than that, so that this is clearly an upper
limit on the pressure of the reevaporation level. An other advantage of this
scheme is that it should be conservative in the sense that it will favor the
formation of unsaturated layers. This is very appropriate considering that we
want to test whether saturation in the stable layer is required to suppress
convection.

Because the WRF dynamical core takes into account the weight of the vapor,
part of the thermal energy absorbed by the atmosphere is converted to
potential energy when this vapor is transported aloft. This is reconverted
back to heat when precipitations fall and dissipate their potential energy
through friction (Forget et al., 2006; Ding & Pierrehumbert, 2016). Therefore,
to close the energy budget, we assume that falling precipitations reach their
equilibrium velocity instantaneously so that the potential energy liberated by
condensates crossing a layer is deposited directly in that layer (See Forget
et al. (2006) for details).

To simulate the fact that the deep atmosphere can act as an infinite reservoir
of vapor, the vapor concentration in the first simulation level is always
restored to the imposed internal vapor concentration, $q_{\mathrm{int}}$, at
each physical timestep. We keep the bottom of the model below the
reevaporation level and verify that this pseudo exchange of vapor with the
interior is negligible when the simulation is equilibrated.

###  2.4 Energy conservation

We observe that there is a net deficit of thermal emission of the atmosphere
compared to incoming radiation, even when the model is equilibrated. In our
baseline simulation, this deficit is about 1-2% of the incoming flux, which
seems reasonable considering the level of accuracy sought. Yet, we here try to
identify what are the possible sources of such energy losses and to identify
possible areas of improvement for the model in the future.

Let us first remind the reader that, although our equations do conserve energy
in the dry gas regime (Laprise, 1992), our dynamical core is not formulated
specifically in an energy-conserving way. This can be seen in Eq. (12), which
uses potential temperature (entropy) as its conserved variable. So some
numerical losses of energy are inevitable. Apart from that, we think the
various losses come from the following sources

  * $\bullet$

Gravity waves, which are launched by updrafts in the dry and moist convective
zones carry away mechanical energy upward. When they reach the top sponge
layer, this energy is dissipated without being reconverted to heat.

  * $\bullet$

To remove the well-known problem of the build-up of energy at the grid scale
caused by the turbulent cascade, the dynamical core implements various filters
that are supposed to dissipate this energy. Here again, the dynamical core has
not been designed to convert this dissipation back to heat, creating an energy
sink.

  * $\bullet$

Although using a potential temperature requires only the ratio $R/c_{p}$ to be
constant, the demonstration of energy conservation in Laprise (1992) relies on
the specific heat capacity of the whole gas also being constant. This is
impossible when the molar mass (hence the specific gas constant) of the dry
gas and vapor are different and that vapor concentration is allowed to change.
Even if we have taken into account the variation of $c_{p}$ in Eq. (12), this
variation should also entail a slight non-conservation, especially when the
mixing ratio of vapor is not negligible.

Table 1: Parameters used in the baseline 3D simulation. Surface gravity | $g$ | [m s-2] | 12.41
---|---|---|---
Average insolation | $F_{\star}$ | [W m-2] | 175.
Deep vapor concentration | $q_{\mathrm{int}}$ | [kg.kg-1] | 0.45
Specific heat capacity | $c_{p}$ | [J K-1kg-1] | 5470.
Molar mass of air | $M_{d}$ | [kg mol-1] | 5.42$\times 10^{-3}$
Surface albedo | $A_{\mathrm{s}}$ |  | 0.

##  3 Atmospheric dynamics and thermal structure of temperate Neptunes

Figure 1: Atmospheric structure of the baseline simulation. From left to
right: Temperature, Water vapor specific concentration, vertical velocity (in
m/s), and relative humidity. The two first panels show horizontal and temporal
averages. The black, dash-dotted line in the first panel show the standard
moist-adiabat profile for the same conditions. In the second panel, the black
dotted and dashed lines show the value of the critical inhibition vapor
concentration (Eq. (1)) and the saturation concentration respectively. The two
last panels show snapshots along vertical slices that go through a moist
convective plume. From bottom to top, the atmosphere exhibits a dry
troposphere, a stable layer where vertical motions are strongly suppressed, a
moist troposphere, and a stratosphere. Horizontal dashed lines are plotted at
the boundaries between these zones to allow an easier comparison of the
altitudes of the various features. The rising moist plume (with maximal
velocities of around 8 m/s) is mirrored by a descending cold plume in the dry
region (-15 m/s) caused by the reevaporation of rains at the bottom of the
stable layer.

We now present our simulations of a prototypical temperate Neptune-like
exoplanet: K2-18b (Cloutier et al., 2017). Indeed, compared to our Solar
System ice giants, the radiative timescale of the atmosphere is on the order
of days to months, instead of decades to centuries. Running a full cloud-
resolving model to thermal equilibrium thus becomes feasible. Despite this
choice, the very fine vertical discretization and the small dynamical step
that it entails make these simulations relatively long and expensive to run:
several months on 20 CPUs for several years of simulated time.

To limit the carbon budget of this study, we thus decided to severely limit
the number of configurations that we simulated in 3D and use these to develop
a 1D model to explore the parameter space as will be presented in Sect. 4. In
this section, we thus delineate the salient features of the thermal,
compositional, and dynamical structure of these atmospheres that need to be
incorporated in the simpler 1D model for it to be realistic.

###  3.1 Simulation setup

The parameters used for our baseline simulation are summarized in Table 1. The
atmospheric composition has been computed for a metallicity of 300$\times$
solar, which is consistent with the observations (Blain et al., 2021). In
principle, the precise atmospheric composition should depend on the
temperature at depth, which itself depends on the modeled convection, hence on
the composition through the mean molecular weight. We thus decided to use a
simple approach to convert metallicity into molecular abundances: we assume
that the quenching in our temperate planet will occur near 10 bar and 1000 K
and use the chemical code of Venot et al. (2020) to compute the abundances.
With this approach, the volumic concentration of the main absorbers are
$x_{H_{2}O}$=2.$\times 10^{-1}$, $x_{CO_{2}}$=1.5$\times 10^{-2}$,
$x_{CH_{4}}$=7.4 $\times 10^{-2}$, $x_{CO}$=2.7$\times 10^{-2}$,
$x_{NH_{3}}$=2.9$\times 10^{-4}$, $x_{He}$=1.5$\times 10^{-1}$, and
$x_{H_{2}}$=5.2$\times 10^{-1}$. For water, this corresponds to a deep
specific concentration $q_{\mathrm{int}}=$0.45 kg/kg, but the concentration in
each cell is traced and computed by the model. The other components are
assumed to remain in same proportions everywhere and to form the so-called dry
air, with a molar mass of 5.42$\times 10^{-3}$ kg mol-1.

Because K2-18b receives an insolation that is very close to the runaway
greenhouse threshold (see discussion in Sect. 6), the extreme climate
sensitivity around this transition makes equilibration of the model very long
for the observed insolation. For this reason, and since the goal of our study
is to understand the general behavior of such atmospheres, we use an average
insolation of 175 W/m2 to run our 3D simulations. This would correspond to an
effective bond albedo of $\approx 0.5$, which allows us to be in the right
regime to test the conclusions of Madhusudhan et al. (2023) in Sect. 6.

###  3.2 Thermal structure

The equilibrated atmospheric structure is depicted in Fig. 1. One of our main
findings is that, despite a more complex humidity distribution in our 3D
simulations, we confirm the structure predicted by Leconte et al. (2017) using
1D idealized simulations: a stable layer forms between the level where
$q=q_{\mathrm{cri}}$ and the dry troposphere below. This stable layer is
clearly noticeable due to i) its superadiabatic thermal gradient and ii) its
very low convective velocities. Indeed, no convective plume, either dry or
moist, does penetrate this layer, so that there is no local source of gravity
waves or turbulence. It can also be noted that the static stability in this
layer is very high so that gravity waves have a much smaller vertical
wavelength compared to the stratosphere.

A new feature of the 3D simulation is that it predicts the existence of a very
thin, dry boundary layer between the stable layer below and the moist
troposphere above. This is reminiscent of the surface boundary layer on Earth
where small-scale dry plumes carry humidity from the surface to the
condensation level where moist convection can occur.

As a comparison, we also show in Fig. 1 the thermal structure for a standard
moist-adiabatic atmosphere with the same insolation and parameters (black
dash-dotted curve). The tropopause exhibits similar pressure levels and
temperatures in both cases, which is to be expected as they have the same flux
to output and similar compositions. However, the temperature at depth varies
dramatically—the moist adiabat being more than 150 K colder at the 1 bar
level. The reason for this is twofold: first, the stable layer created a huge
temperature jump in a narrow vertical region and, second, because the
atmosphere below the stable region is dry, the thermal gradient follows the
dry adiabat, which increases the temperature faster with depth than the moist
adiabat. Because the moist region extends much deeper in the standard model,
the temperature difference at depth will be even larger, as will be shown in
later sections.

###  3.3 Vapor cycle and energy budget

We find that the thermal structure is strongly linked to the vapor cycle of
the condensable species. The most apparent manifestations of this cycle in the
simulations are the moist convective plumes (See Fig. 1) that transport vapor
upward in the moist troposphere. These plumes transport sensible and latent
heat, as can be seen in Fig. 2. This upward motion is well-known for creating
large amounts of condensates and precipitations. However, note that a sizable
fraction of the vapor directly condenses at the top of the stable layer,
forming a thin, horizontal cloud deck. This is evidenced by the sudden drop in
the latent heat flux just at the top of the stable layer in Fig. 2.

Rainfalls then transport the condensable species back down in its condensed
phase. The average rainfall rate can be estimated from the latent energy flux
through the stable layer divided by the latent heat of vaporization and is
about 10-5kg/s. As could be expected, the reevaporation of these rainfalls
occurs below the base of the convective plumes. Indeed, a first requirement
for reevaporation is to reach an unsaturated region below the updraft.

Figure 2: Net upward energy fluxes due to sensible heat dynamical transport
(solid) and latent heat (dashed). In the stable region, latent heat
transported by vapor dominates the direct sensible heat flux due to
turbulence, although the latter is not nil. See Appendix A for details of the
calculation method.

An interesting finding is that the level of reevaporation seems to set the
bottom of the stable layer. Indeed, below this level, condensation is
impossible, so there are no vapor sources/sinks. As a result, any vapor
gradient in this lower region – called dry troposphere in Fig. 1 – would be
short-lived and stellar energy deposited at depth is sufficient to trigger
standard, dry convection and mix efficiently the atmosphere. But why cannot
the dry troposphere extend above the reevaporation level? Let us remember that
in a standard atmosphere, what drives convection in the first place is the
fact that thermal radiation is not sufficient to carry the absorbed stellar
radiation away, and convection stops at levels where this radiative cooling
becomes efficient enough. Here, it takes a lot of energy to evaporate all the
falling precipitations. This efficiently shuts down convection at the
reevaporation level and energy is mostly carried upward by vapor in the form
of latent heat above it (See Fig. 2).

To close the cycle, vapor needs to find its way back toward the moist
troposphere through the stable layer. This was a shortcoming of the model of
Leconte et al. (2017): they did not consider that the precipitations formed in
the moist convective layer would necessarily reevaporate below the stable
layer, so that the vapor would need to diffuse upward to maintain moist
convection. In particular, this would have been in contradiction with their
findings that no double-diffusive instability would develop in the stable
layer, so that energy transport would be purely radiative.

Our simulations shed a completely different light on this issue. Although we
cannot test directly the ”no double-diffusive instability” hypothesis because
it would require a much finer resolution (Rosenblum et al., 2011), this issue
is rendered rather moot by the other sources of turbulence in the system.

###  3.4 Turbulent mixing

Even at the $\mathcal{O}(500\,m)$ scale resolved in our simulations,
turbulence spontaneously appears111Turbulence here specifically refers to
small-scale, but resolved motion that appears in stratified, stable regions..
Even though the velocities involved are rather small (0.1 m/s; see Fig. 3),
transport is still significant thanks to the steepness of the vapor and
potential temperature gradient and the rather small vertical extent of the
stable zone.

To quantify this, we estimate the equivalent mixing coefficient (the so-called
eddy diffusivity or $K_{zz}$) in the simulations using a passive tracer whose
concentration ($r_{\mathrm{tra}}$) is fixed at the surface and that undergoes
advection by the flow and local exponential decay with a timescale
$\tau_{\mathrm{tra}}$. Then we use two different methods. In the first one, or
eddy flux approach, we just compute the average turbulent flux of tracer in
the simulations and assume that:

$\displaystyle K_{zz}^{\mathrm{Eddy\,Flux}}\equiv\frac{\langle\rho
r_{\mathrm{tra}}w\rangle}{\rho\langle\partial_{z}r_{\mathrm{tra}}\rangle},$
(18)

where $\langle\rangle$ denotes temporal and horizontal averaging. In the
second method, or integral approach, we use the fact that if the mixing
exhibits a diffusive-like behavior, the steady-state tracer profile should
obey the following law

$\displaystyle\frac{1}{\rho}\partial_{z}\left(\rho
K_{zz}\partial_{z}r_{\mathrm{tra}}\right)=\frac{r_{\mathrm{tra}}}{\tau_{\mathrm{tra}}}.$
(19)

Integrating from the top of the atmosphere where we know that there is no
vertical flux, one can get a second estimate

$\displaystyle
K_{zz}^{\mathrm{Integral}}\equiv\frac{1}{\tau_{\mathrm{tra}}}\frac{\int\rho
r_{\mathrm{tra}}\mathrm{d}z}{\rho\langle\partial_{z}r_{\mathrm{tra}}\rangle}.$
(20)

Note in passing that because the transport is not perfectly diffusive
(especially in convective regions), this approach can yield negative $K_{zz}$
so that we show the absolute value. With this important caveat in mind, both
estimates (shown in Fig. 3) exhibit a relatively good agreement and the
differences inform us on the uncertainty that can be attributed to this
parameter.

Figure 3: Left: Vertical profile of the equivalent vertical mixing coefficient
derived from the simulation using both the Eddy flux (Eq. (18)) and Integral
(Eq. (20)) methods . The dotted line shows the $1/\sqrt{\rho}$ trend for
comparison. The horizontal lines depict the same levels as in Fig. 1. Right:
Profile of the vertical velocity showing the root-mean-square (solid), maximum
upward (dotted) and downward velocity (dashed). Averages and maximum values
are computed over temporal and horizontal dimensions.

The profile of eddy diffusivities is fully consistent with the thermal profile
found: stable layers exhibit low diffusivities whereas convective ones are
more strongly mixed. The values of the Eddy diffusivity in the stable regions
is on the order of 0.1-10 m.s-2, which is comparable to values inferred for
the lower stratosphere on Earth (Allen et al., 1981). As we will discuss in
more details in Sect. 4, the vapor profile shown in Fig. 1 seems to favor eddy
diffusivities that are at the lower end of this range and decreasing slightly
with altitude in the stable region. This is evidenced by the vapor gradient
steepening with altitude as we go from the bottom to the top of the stable
layer. Indeed, there is no condensation or reevaporation in that region so
that, in steady state, the upward vapor flux ($\propto
K_{zz}(\mathrm{d}q/\mathrm{d}z)$) needs to be constant. So a steepening
gradient means a decreasing mixing.

This turbulence seems to be produced by the dry updrafts hitting the bottom of
the stable layer (Lane et al., 2003). This both creates upward propagating
gravity waves and directly stirs the medium. We believe that gravity waves
themselves do not participate much in the mixing in the stable region because
i) they do not transport matter in the linear regime and ii) their amplitude
is too low at this level to break and induce subsequent mixing.

###  3.5 Dynamics

The velocity distributions shown in the right panel of Fig. 3 are very
different between the two convective zones. In the bottom, dry convective
region, we observe a permanent overturning circulation with velocities on the
order of 1 m/s. At this depth, this is largely driven by stellar radiation
that still penetrates efficiently. On the contrary, in the moist convective
region, there is very little background motion apart from upward propagating
gravity waves that still have a low amplitude as they are not far from their
launching region (a few 10-1 m/s). The only exception occurs in the thin
boundary layer just above the stable layer below, around $2\times 10^{4}$ Pa.
We can see a local maximum of both the RMS velocity and the eddy diffusivity
in Fig. 3. In this region, the atmosphere is locally unstably stratified,
driving small-scale dry convection which transports humidity upward until
saturation induces moist convection.

Then, as can be seen in Fig. 1, moist convective plumes form episodically to
release the energy and vapor that has built-up at the top of the boundary
layer. As is common on Earth, these moist convective plumes are much narrower
and reach faster speeds than their dry counterparts. As also seen on Earth,
the rising plumes are mirrored by cold downdraft created by reevaporation
below the convective cloud: the cold pools. One big difference is that on
Earth, the moist updrafts reach higher velocities than the cool downdrafts.
This is because in our atmosphere, updrafts are powered by the combined power
of latent heat and compositional buoyancy (remember that water vapor is
lighter than molecular nitrogen) whereas these two effects compete in
downdrafts, with the thermal effect of latent heat dominating.

In a hydrogen-dominated atmosphere like K2-18 b’s, vapor is usually heavier
than dry air. So latent heat release by condensation needs to fight against
the stabilizing effect of the mean molecular weight gradient to power the
rising, moist plumes. In fact, it can win only when the vapor concentration is
below the critical ratio, which is the very reason why there is a stable
region in the first place. But even in the moist region, this competition
leads to relatively sluggish upward motion. In comparison, the downdrafts that
form below the rising plumes are much more vigorous, with a factor
$\times\,$2-3 in velocity, because the vapor loading helps the cooling due to
reevaporation (that can reach several K) in accelerating the plume downward.

Finally, one can see chevron-like structures around the rising plume in the
third panel of Fig. 1. These are typical of convection induced gravity waves
that can propagate as the moist region is stably stratified (in the sense that
it has a positive Brunt-Väisälä frequency).

###  3.6 Stabilisation and subsaturation

An intriguing feature of our simulations is that the atmosphere is
subsaturated almost everywhere but for a thin region near the top of the
stable layer and in the core of ascending moist plumes. Yet, the atmosphere
exhibits a stable layer where vapor concentration exceeds the critical
inhibition fraction as predicted by Leconte et al. (2017). These two
statements seem in contradiction as the analytical theory invokes saturation
to suppress convection. However, as discussed in Sect 3.4 of Leconte et al.
(2017), saturation is invoked only to suppress moist convection, and
saturation indeed needs to occur to form moist convective updrafts as visible
in the third panel of Fig. 1.

However, the medium does not need to be saturated for the vapor gradient to
have a stabilizing effect—thermohaline convection being a perfect example
(Ledoux, 1947; Stern, 1960). So saturation needs to happen somewhere in the
convective parts of the atmosphere to drive the vertical vapor concentration
gradient between a low value above the cloud deck and a high internal value.
But our simulations confirm that the stable layer can be largely subsaturated
and remain stable.

Interestingly, we find that the saturated layer that coincides with the top of
the stable layer happens exactly at the critical concentration value. Again,
this makes sense because above that level, moist convection can occur. Moist
convective regions must be subsaturated on average (see Fig. 1) because of the
dry subsidence regions that appear to compensate for the upward mass flux in
the convective plumes. In the stable region, where there is very little
motion, the vapor concentration is much more horizontally and temporally
uniform. Hence it is always close to saturation at the layer top because this
is where moist plumes originate.

##  4 A simplified 1D framework for fast modeling

Our 3D simulations are too expensive to be carried over a large diversity of
conditions. Now that we have outlined the most important features of the
structure of Neptune-like atmospheres, we present a 1D model that is able to
reproduce these features for only a tiny fraction of the computational burden.
This model is based on the Exo_k library (Leconte,
2021)222https://perso.astrophy.u-bordeaux.fr/~jleconte/exo_k-doc/index.html
that has been recently updated with a full-fledged time-stepping atmospheric
evolution package described in Selsis et al. (2023).

This atmospheric evolution package has the advantage of being extremely
flexible while using some computational tricks to remain very fast. Hereafter,
we will focus mainly on the new features of the model that pertain to the
inhibition of convection.

Figure 4: Atmospheric structure of the 1D model for the baseline case. From
left to right: Temperature, Water vapor specific concentration, dynamical
fluxes, and radiative fluxes vertical profiles. The dashed curve in the
temperature panel shows the average 3D thermal profile for comparison. In the
second panel, the dashed line shows the value of the saturation concentration
ratio and the red curve shows the profile of eddy diffusivity. Dots show
convective zone where mixing is not diffusive so that the value of the
diffusivity is irrelevant there. The third panel can be compared to Fig. 2.
The last panel shows (minus) the net upward incoming stellar flux (solid) and
the net outgoing thermal flux emitted by the atmosphere (dashed). This latter
panel shows that most of the radiative cooling to space occurs in the stable
layer and above.

###  4.1 Criteria for convection inhibition

To allow for relatively large timesteps, Exo_k treats both moist and dry
convection using standard adjustment schemes that identify sets of adjacent
layers that are unstable and compute the energy fluxes necessary to bring back
these layers to the relevant adiabat.

For dry convection, unstable layers were identified by regions of decreasing
potential temperature ($\theta$). This is equivalent to the Schwarzschild &
Härm (1958) criterion. When one wants to account for the mean molecular weight
effect of the vapor, one needs to use the Ledoux (1947) criterion instead. But
as is well known in terrestrial meteorology (see also Leconte et al. 2017),
this can simply be recast as identifying regions of decreasing _virtual_
potential temperature,

$\displaystyle\theta_{\mathrm{v}}\equiv\left(1-\frac{M_{v}-M_{d}}{M_{v}}q_{v}\right)\
\theta,$ (21)

inducing minimal changes to the code. Any set of unstable layers is brought
back to a neutral state with a uniform composition and potential temperature
(hence uniform virtual potential temperature) in a single timestep and
conserving the total enthalpy of the layer. Layers above or below the newly
adjusted layer that have been destabilized by the adjustment are themselves
adjusted iteratively. This treatment is a bit simpler than the one proposed in
Habib & Pierrehumbert (2023) for dry compositional convection. Our scheme
nonetheless performs very well in the bottom dry-convective region, as shown
hereafter. This is probably due to the fact that mixing in this layer is
continuously driven by thermal and compositional effects and is thus
relatively efficient.

In the previous version of Exo_k (Selsis et al., 2023), moist convective
adjustment was triggered when the temperature gradient was larger than the
moist adiabat and when the medium was saturated (Manabe & Wetherald, 1967).
Following the analysis in Leconte et al. (2017), moist-convection inhibition
has been accounted for by simply suppressing the adjustment in any layer where
$q_{v}>q_{\mathrm{cri}}$. In the convective regions, the excess vapor
condensed during the adjustment process is assumed to instantaneously rain
down to the reevaporation level.

In essence, these changes are sufficient to naturally force a stable layer in
our unidimensional model. This however requires the ability to trace the cycle
of both the vapor and the condensates. For this we use two tracers that are
mixed by convection as described above. For the thermodynamics and the
microphysics of clouds and precipitations, we use exactly the same
parametrization as in our 3D model, which is described in Sect. 2.3, ensuring
that any difference will be due to the dynamics.

###  4.2 Turbulent mixing in stable layers

While it is not needed to create a stable layer, we have seen in Sect. 2 that
turbulent mixing is an important ingredient in determining the strength of the
vapor cycle and, to a lesser extent, in transporting sensible heat through the
stable layer. To incorporate that, we added a diffusive flux of tracers and
entropy of the form $K_{zz}(\mathrm{d}/\mathrm{d}z)$, where the eddy
diffusivity, $K_{zz}$, is a free parameter. In our simplest model, hereafter
called the constant diffusivity case, this parameter is constant throughout
the atmosphere and calibrated to yield the proper flux of vapor through the
stable layer. This yields a value of about $0.3$ m2/s, which is representative
of the values found in the stable layer in Fig. 3. Remember that tracers are
already fully mixed in the convective regions, which emulates an infinite eddy
diffusivity. So there is no need to increase our $K_{zz}$ there.

Yet, we find that if one wants to model very accurately the shape of the
thermal profile in the stable zone, a constant $K_{zz}$ provokes too abrupt a
transition at the top of the dry convective region, where convective plumes
can overshoot. To model this, we implement an alternative formulation,
referred to as the baseline scenario, where the eddy diffusivity assumes a
larger value ($K_{zz}^{\mathrm{max}}$) just at the top boundary of the dry
convective region (also called the radiative-convective boundary, whose
pressure is $p_{\mathrm{RCB}}$) and drops off very rapidly as some high power
$\alpha$ of the pressure before settling to a lower constant value
($K_{zz}^{\mathrm{min}}$) higher up:

$\displaystyle
K_{zz}=\left\\{\begin{array}[]{ll}K_{zz}^{\mathrm{max}}&p>p_{\mathrm{RCB}}\\\
\mathrm{Max}(K_{zz}^{\mathrm{max}}(p/p_{\mathrm{RCB}})^{\alpha},K_{zz}^{\mathrm{min}})&p<p_{\mathrm{RCB}}.\end{array}\right.$
(24)

The value of $\alpha$ is set to reproduce the sharp decrease of the turbulence
above the convective region seen in Fig. 3. This yields $\alpha\approx 13$.
Then, $K_{zz}^{\mathrm{min}}$ and $K_{zz}^{\mathrm{max}}$ are tuned to
reproduce the vapor flux and the thermal gradient just above the radiative-
convective boundary, respectively. This yields $K_{zz}^{\mathrm{min}}\approx
0.08$ m2/s and $K_{zz}^{\mathrm{max}}\approx 3$ m2/s, which is also consistent
with the numerical values from the 3D simulation.

###  4.3 1D/3D comparison

In Fig. 4, we present the atmospheric structure for the baseline scenario.
Notice that because there are some energy losses in the dynamics of the 3D
model (see Sect. 2.4), the baseline case is ran with a decrease of the
incoming stellar flux of the corresponding amount to allow for a proper
comparison.

The agreement between the 1D and 3D models is rather striking. With very few
free parameters, not only both thermal structures are very close, but the
shape and magnitude of the various energy fluxes are reproduced as well. We
even recover the thin, dry boundary layer between the top of the stable layer
and the moist convective layer, which is evidenced by the small layer with
constant vapor concentration near $2\times 10^{4}\,$Pa in the second panel of
Fig. 4. This shows that this boundary layer is not created by dynamical
requirements alone. The most notable discrepancy is the fact that the moist
convective layer is fully saturated in our 1D model, which is to be expected
because saturation is a prerequisite to the onset of convection in our moist
adjustment scheme. Our 3D simulations are closer to what happens in reality
where saturation is only required in the rising plumes and dry subsident
regions force the convective region to be subsaturated on average.

Fig. 5 shows how the thermal profile is modified when changing the profile of
the turbulent eddy diffusivity. As expected, the change in the slope of the
thermal gradient is more abrupt at the top of the dry convective zone in the
case of a constant eddy diffusivity. This changes the 1 bar temperature by
about 20 K. We also ran the baseline model but with the same actual input flux
as the 3D simulations to quantify the potential effect of the dynamical losses
of the latter. This also causes an increase of the deep adiabat of about 20 K.
Although not negligible, we have to bear in mind that these differences are
much smaller than the effect of the inhibition itself, which raises the
temperature at the 1 bar level by $\sim$200 K.

Figure 5: Sensitivity analysis showing the impact of various assumptions on
the equilibrated temperature profile. The solid curve is the baseline
profiles, the dashed one shows the case with a constant eddy diffusivity, and
the dotted one shows the effect of correcting for the dynamical losses.
Differences between the various cases is much smaller than the effect of
convection inhibition itself which can be seen by comparing with the no-
inhibition case (dash-dotted curve).

##  5 Observational markers of convection inhibition

As can be seen in Fig. 5, the presence or absence of a stable layer (i.e. of
convection inhibition) affects relatively mildly the temperature of the upper
troposphere. This is because the atmosphere above the stable layer is already
optically thick in most of the thermal part of the spectrum. As a result, the
temperature of the atmosphere below the stable layer affects only slightly the
outgoing flux and only a small change in the temperature of the photosphere is
needed to reach global radiative equilibrium.

To identify a more robust marker of convection inhibition in temperate sub-
Neptunes, we turn to the composition of the atmosphere. Indeed, we know that
the chemical composition of temperate atmospheres connected to a deep gaseous
envelope is mainly determined by the temperature of the level at which
chemical processes are quenched (Venot et al., 2018). Above the quench level,
chemical reactions are too slow compared to the mixing by the atmosphere.
Above that level, the composition is thus rather uniform up to the level where
photochemical rates start to dominate. For a given elemental abundance, the
molecular content of the atmosphere can therefore be a tracer of the deep
temperature.

The main implication of the presence of a stable layer is the higher
temperature in the deep atmosphere. This is particularly visible in the right
panel of Fig. 6 where we have extended our 1D model to 1000 bars. The increase
in temperature is due to two effects: i) the temperature jump in the stable
layer and ii) the fact that the dry troposphere, which has a much steeper
lapse rate, starts much higher. Hence, differences can be up to 1000 K at the
100 bar level.

To quantify the effect on the chemistry, we have computed the chemical
composition of the two model atmospheres in Fig. 6 using the chemistry module
of Exo-REM (Blain et al., 2021) assuming $K_{zz}=10^{4}\,$m2/s. We take this
value as characteristic of the deep convective region where the quenching
occurs. Our tests show a rather low sensitivity to the exact value of this
parameter. Exo-REM computes the quenching levels and resulting disequilibrium
composition by comparing the mixing timescale with the chemical timescale
using formulas from Zahnle & Marley (2014). This yields a much lower quenching
pressure of $\approx 20\,$bar for the baseline scenario compared to $\approx
600\,$bar for the no-inhibition case.

As can be seen in the left panel of Fig. 6, the two cases predict very
different chemical abundances in the stratosphere. Notably, the no-inhibition
case - with its low temperature interior - exhibits very low levels of CO and
CO2 in the stratosphere, most of the carbon forming CH4. On the contrary the
case where inhibition is taken into account, CO, CO2, and CH4 are all in
detectable quantities. This case also predicts a much lower abundance of NH3.
These conclusions are on par with the conclusions of Cavalié et al. (2017) for
Uranus and Neptune.

As water vapor is cold-trapped at the tropopause in all cases, carbon bearing
species dominate the transit spectrum. The two model atmospheres described
above would thus be easily distinguishable observationally with JWST.
Observing a Neptune-like planet receiving an insolation below the runaway
greenhouse threshold would thus enable us to infer the presence of convection
inhibition and a stable layer.

Figure 6: Vertical profiles of molecular volume mixing ratios in the
atmosphere of our prototypical temperate Neptune with (solid) and without
(dotted) convection inhibition. The metallicity is assumed to be
300$\times$solar. The right panel shows the thermal profiles used for the two
cases, where we see that the moist troposphere extends quite deep in the no-
inhibition case, resulting in a much lower quenching temperature.

##  6 Atmospheric constraints on the existence of liquid oceans on K2-18 b

It has been recently claimed that the absence of ammonia in the transit
spectrum of K2-18 b could be the sign of the existence of a liquid water ocean
below a relatively shallow atmosphere of less than a few bars (Madhusudhan et
al., 2023).

The main argument against such a scenario is that the irradiation level
received by the planet is above the critical runaway greenhouse threshold for
cloudless H2 rich atmospheres (Innes et al., 2023). But Madhusudhan et al.
(2023) argues that tropospheric clouds or hazes could increase the albedo of
the planet, stabilizing the liquid ocean.

In this section, we revisit these arguments using our newly developped 1D
atmospheric model and the improved knowledge of the atmospheric composition of
the planet provided by the recent JWST observations of the system (Madhusudhan
et al., 2023). In particular we show that

  * $\bullet$

When convection inhibition is accounted for, the planetary albedo required to
stabilize an ocean on K2-18 b is higher than previously estimated,

  * $\bullet$

Unlike on solar system planets, tropospheric clouds cannot provide such high
albedos because a significant part of the stellar flux does not reach the
troposphere to be reflected,

  * $\bullet$

When we add sufficient levels of stratospheric haze to our model atmospheres
to reach the albedos necessary to keep a liquid ocean, the methane features in
the transmission spectrum are too muted to be consistent with recent
observations.

###  6.1 Runaway greenhouse threshold and convection inhibition

Our first goal is to assess the maximum stellar irradiation K2-18 b can
receive while still sustaining a liquid ocean below the H2 dominated
atmosphere. This question has been investigated by Innes et al. (2023), who
showed that the greenhouse effect of an H2 dominated atmosphere is much
greater than that of an N2 dominated atmosphere, strongly lowering the runaway
greenhouse threshold. We revisit here these calculations for several reasons:
i) we want to use the exact planetary parameters of K2-18 b, ii) JWST data now
provide an estimate of the atmospheric composition of the planet, which allows
for more accurate opacity and mean molecular weight estimates, and iii) our 3D
simulations provided a better understanding of the turbulent transport
mechanisms in those atmospheres, whereas Innes et al. (2023) assumed that
energy was solely transported by radiation in the stable layers of the
atmosphere.

Because the irradiation received by the planet is rather well constrained, we
reframe the issue by asking what would the planetary Bond albedo need to be
for the absorbed insolation to remain below the runaway greenhouse threshold.
To answer this question, we perform the following experiment : we equilibrate
our 1D model for a given surface pressure with a stellar irradiation equal to
$F_{inc}(1-A)$ where $F_{inc}=342\,$W/m2 is the average insolation received by
K2-18 b and $A$ is an albedo that we choose arbitrarily. In this specific
experiment, we do not include any aerosols in the atmospheric model as their
effect is wholly incorporated in the $A$ parameter. Because K2-18 is a red
star and that significant amounts of methane have been detected in the
planet’s atmosphere, almost all the stellar light arriving at the atmospheric
top is absorbed in this setup so that $A$ is a good proxy for the bond albedo
of the model.

There are two differences with respect to 1D simulations shown in the previous
sections. First, we recomputed opacity tables with mixing ratios of methane
and carbon dioxide of 10-2ppmv, which seem to be the best match for the JWST
observations of Madhusudhan et al. (2023). Second, we change the surface
boundary condition to better mimic an ocean. Instead of fixing the mixing
ratio of vapor in the lowest layer equal to an expected value at depth, we now
treat the surface as an infinite source of water and vapor can freely
evaporate in the first layer until saturation is met. As in Innes et al.
(2023), we keep the total surface pressure constant during the evolution. We
then let the model evolve until it either reaches thermal equilibrium, in
which case we deem the ocean stable, or the mixing ratio of vapor reaches
unity at the surface, which we take as a proxy for the onset of runaway
greenhouse.

The results of this experiment are summarized in Fig. 7. As expected, the
higher the surface pressure, the higher the Bond albedo needs to be to keep a
stable surface ocean. These results are in rough agreement with the results of
Innes et al. (2023). First, we confirm that convection inhibition decreases
the threshold for the onset of runaway greenhouse. This is because near the
critical insolation, the moisture always becomes sufficient near the surface
to shut down moist convection, increasing the surface temperature for a given
outgoing flux, thus increasing the greenhouse effect of the atmosphere.
Second, Innes et al. (2023) find that the maximum stable insolation for a 1
bar atmosphere around an M star is around 110 W/m2 whereas our last stable
insolation is $\approx$ 130 W/m2 ($A=0.6$). The discrepancy could be partly
due to the presence of methane, whose anti-greenhouse effect can be rather
strong around late-type stars, but we think that the main difference is our
treatment of the turbulent transport in the stable layer where convection
inhibition operates. Because turbulence transports both sensible and latent
heat, the thermal gradient is much less steep in our model compared to a fully
radiative zone, which weakens the greenhouse effect of the stable layer. We
verify this by re-running this case with a much weaker turbulent transport and
find that it indeed enters a runaway greenhouse phase. Further comparison is
however difficult as the model currently relies on turbulence to transport
water vapor in non-convective zones so that removing it entirely causes the
stable layer to become unsaturated. But this shows that accounting correctly
for the dynamics of the atmosphere is important to derive quantitative limits.

However, the limits we find are much more stringent than the ones found by
Madhusudhan et al. (2021). For an M star like K2-18, they find that the
maximum equilibrium temperature (that is corrected for the Bond albedo) to
keep a liquid surface ocean is $\approx$ 410 K — which corresponds to a planet
averaged absorbed/thermally emitted flux of $\approx$ 1300 W/m2. This is to be
compared to our limit for the 1 bar case, which is estimated to be $\approx$
230 K (150 W/m2). To put this into context, the current equilibrium
temperature of the Earth is $\approx$ 255 K (240 W/m2), and recent estimates
of the runaway greenhouse limit for Earth-like planets yield estimates between
$\approx$ 260 and 270 K (270-300 W/m2), depending on the treatment of
continuum opacities, clouds, and atmospheric dynamics (Kopparapu et al., 2013;
Leconte et al., 2013; Yang et al., 2016). The fact that the limit for
hydrogen-dominated atmospheres occurs at lower fluxes is due to the increased
greenhouse effect of H2 compared to N2, and has been extensively studied (Koll
& Cronin, 2019; Chaverot et al., 2022; Innes et al., 2023).

The reason that Madhusudhan et al. (2021) find such high limits is less clear.
It seems to be due to their use of an ad hoc – and rather extreme –
approximation to treat aerosols: they assume that aerosols can be arbitrarily
efficient scatterers and model them by multiplying the Rayleigh scattering
coefficient of H2 by an arbitrarily large factor until the atmospheric Bond
albedo reaches the desired value. In addition to increasing the albedo, this
causes the stellar radiation to be scattered many times in the stratosphere,
which, counter-intuitively, enhances absorption there. Around redder stars,
this results in stratospheres that are as hot, if not hotter, than the
surface, which effectively suppresses the greenhouse effect of all the
atmospheric gases.

However, as we will see in the next section, the presence of such reflective
haze particles in the stratosphere is contradicted by observational data.

Figure 7: Constraints for the presence of a liquid surface ocean. Each marker
shows the outcome of a simulation for a given imposed albedo (the equivalent
effective flux is shown on the top axis) and atmospheric surface pressure.
Blue dots show cases where a liquid surface ocean is stable. Red crosses show
cases where a steam atmosphere forms. The top panel shows the results of
traditional models where convection inhibition is disregarded and the bottom
panel shows results with convection inhibition. The dashed (solid) line
roughly depicts the limit to the stability of an ocean in the case without
(with) inhibition. One can see that the inhibition limits the stability of
oceans to higher albedos, i.e. less irradiated planets.

###  6.2 Constraints on aerosols

In this section, we make an attempt at better quantifying the limits that can
be put on the albedo that aerosols (either clouds or hazes) can realistically
produce on a planet like K2-18 b.

A first hypothesis put forward by Madhusudhan et al. (2023) is that the
presence of deep, highly reflective tropospheric clouds could produce a
sufficient albedo to stabilize an ocean. Although such clouds are able to
produce high albedos for our solar system giant planets, we have to remember
that K2-18 is an M dwarf with an effective temperature around 3500 K and that
its radiation is easily absorbed in the stratosphere of the planet by the
multiple near-infrared methane bands. In a cloudless model of K2-18 b produced
with the fiducial methane and carbon dioxide mixing ratios of 10-2 ppmv found
by Madhusudhan et al. (2023), half the flux is absorbed above the 100 mb
level, which is still higher than the tropopause. So no scattering happening
below this level, however intense, could increase the albedo above 0.5 (and
that does not even account for the fact that scattered light would have to
cross the stratosphere a second time to escape). This is well illustrated by
the 3D global climate models from Charnay et al. (2021) who found that the
albedos of their models for K2-18 b barely exceed 0.12 even when thick dayside
tropospheric water clouds form.

Another hypothesis is the presence of highly reflective haze in the
stratosphere, although neither Madhusudhan et al. (2021) or Madhusudhan et al.
(2023) discuss which type of haze could meet the required constraints. This
solution works in principle because it is able to scatter incoming stellar
light high in the stratosphere, before it is efficiently absorbed. However, it
is easy to see that such a reflective haze should also strongly affect (e.g.
flatten) the transit spectrum of the planet in the visible and near-infrared,
whereas the recent JWST spectrum of Madhusudhan et al. (2023) shows relatively
deep methane absorption features with an amplitude in excess of 100 ppm
between 1 and 1.5 micron. To quantify this effect, we compute eclipse and
transmission spectra of the fiducial model of K2-18 b discussed above where we
add haze scattering following the parametrized approach of Madhusudhan et al.
(2021). In this approach, the amount of haze is encompassed in a so-called
haze factor ($n_{haze}$), which is used to multiply the cross section of
Rayleigh scattering of H2. Let us note that $n_{haze}=1$ corresponds to the
fiducial clear atmosphere model.

Figure 8: Eclipse (top) and transit (bottom) spectra of our models of K2-18 b
with hazes parametrized through the $n_{haze}$ factor (see text). One can see
that when the amount of haze increases, the amount of reflected light (hence
the albedo) increases but the amplitude of the methane bands in the transit
spectrum decreases.

The resulting spectra are shown in Fig. 8, where we see that the amount of
reflected light increases with $n_{haze}$, as expected. The corresponding
albedo can be seen in Fig. 9, with $A=0.03$ for the clear case and $A=0.72$
for $n_{haze}=10^{5}$. However, one can see in the bottom panel of Fig. 8 that
the increased scattering also mechanically suppresses the molecular methane
features in the near-infrared. We quantify the amplitude of these features by
taking the difference between the transit depth at two absorption peaks (1.
and 1.16$\mu$m) and the depth in the nearest windows (1.08 and 1.28$\mu$m
respectively). Those amplitudes are shown as a function of the model albedo in
Fig. 9. As could be expected, the amplitude of the molecular features
decreases when the albedo increases. But more importantly, we cannot find a
model where we have both a sufficient albedo and a transit amplitude that
matches the data.

Because the aforementioned parametrization of haze is rather ad-hoc, we have
tested various types of aerosols, including water ice particles and venusian
sulfuric acid cloud particles that are known for their high reflectivity.
Although we do not show all the results here, we always find a very similar
trend to what is shown in Figs. 8 and 9: models with high albedos produce very
flat spectra that do not match the observations. And it is difficult to
imagine a species that would reflect enough light in an unobserved part of the
spectrum because the NIRISS SOSS observation precisely cover the peak of the
stellar emission.

So we conclude that the observations of K2-18 b make the possibility of a
planet harboring a liquid ocean thanks to a haze-driven high-albedo very
unlikely. The only remaining possibility would be that hazes would form only
on the dayside to almost disappear at the limbs. However, the reader should
bear in mind that, unlike clouds, hazes cannot easily sublimate and thus
usually are much more uniformly distributed than clouds – as shown by solar
system examples like Titan or Giant Planets. We thus deem this possibility
rather unlikely as well.

Figure 9: Amplitude of two near-infrared methane molecular bands in the
transit spectrum of our models of K2-18 b with parametrized hazes as a
function of the Bond albedo. As already illustrated in Fig. 8, the amplitude
of the methane bands decreases when the albedo of the planet increases. The
numbers show the values of $n_{haze}$. The shaded area roughly depicts the
area of parameter space that could be consistent with both the observations
(transit amplitude greater than 100 ppm) and the albedo required to sustain a
surface ocean ($A\gtrsim 0.5$). The relatively high band amplitude found in
the observations rules out haze that are reflective enough to stabilize a
liquid surface ocean.

##  7 Conclusions

We have developed a cloud-resolving model able to simulate moist-convection in
vapor enriched atmosphere. Being integrated in the ecosystem of the Generic
PCM model, it is very flexible and can be easily adapted to a wide diversity
of planets. We then have investigated how moist convection behaves in H2
dominated atmospheres using K2-18 b as a prototypical temperate Neptune-like
planet.

Our main general findings are that

  * $\bullet$

Moist convection is effectively inhibited when the vapor abundance exceeds the
threshold abundance given by linear theory (Guillot, 1995; Leconte et al.,
2017)

  * $\bullet$

The atmospheric structure envisioned by Leconte et al. (2017) – i.e. the
formation of a stable layer between a moist troposphere above and a dry
troposphere below – is recovered in the 3D simulations, even though some of
the fundamental assumptions of the analytical theory are not verified. In
particular, almost no part of the atmosphere is fully saturated.

  * $\bullet$

The stable layer harbors some turbulence, the magnitude of which controls the
intensity of the vapor cycle in the atmosphere. Both the latent and sensible
heat transport that result contribute significantly to the energy flux through
the stable layer and need to be accounted for when determining the thermal
structure.

  * $\bullet$

The deep gaseous envelope of Neptune-like planets where condensation occurs
should be much hotter than envisioned by standard models. This impacts the
chemical composition of the atmosphere. This could itself be a way to
ascertain the presence of a stable layer in the atmosphere.

  * $\bullet$

The higher temperature at depth also decreases the limiting insolation at
which a surface liquid ocean can be stable under a H2 dominated atmosphere
(also see Innes et al. 2023).

Although our conclusions are based on simulations that focus on the water
condensation region, the principles are rather general and should readily
apply to the condensation level of any other condensable species that is
abundant enough: methane in Uranus and Neptune (Clément et al., 2023), iron or
silicates near the core of Neptune-like planets (Markham et al., 2022; Misener
& Schlichting, 2022), etc.

The fact that the higher temperatures at depth further reduce the insolation
at the so-called inner edge of the habitable zone for H2 dominated atmospheres
has direct implications for observed planets that are in this range of
insolations. In particular, the insolation received by K2-18 b seems too high
to find a configuration that could explain both the very high albedo necessary
to stabilize a surface ocean under its H2 atmosphere and the rather deep
methane features observed in transmission spectroscopy. Therefore, if the non-
detection of ammonia in this atmosphere is confirmed, it might be necessary to
invoke other mechanisms for the lack of this molecule than a shallow
atmosphere above a liquid ocean, like

###### Acknowledgements.

The authors acknowledge the support of the French Agence Nationale de la
Recherche (ANR), under grant ANR-20- CE49-0009 (project SOUND), the Programme
National de Planétologie (PNP), and CNES. ML acknowledges funding from the
European Union’s Horizon Europe research and innovation program under the
Marie Sklodowska-Curie grant agreement ”MuSICA-V”

## References

  * Allen et al. (1981) Allen, M., Yung, Y. L., & Waters, J. W. 1981, J. Geophys. Res., 86, 3617
  * Blain et al. (2021) Blain, D., Charnay, B., & Bézard, B. 2021, A&A, 646, A15
  * Bryan & Fritsch (2002) Bryan, G. H. & Fritsch, J. M. 2002, Monthly Weather Review, 130, 2917
  * Cavalié et al. (2020) Cavalié, T., Venot, O., Miguel, Y., et al. 2020, Space Sci. Rev., 216, 58
  * Cavalié et al. (2017) Cavalié, T., Venot, O., Selsis, F., et al. 2017, Icarus, 291, 1
  * Charnay et al. (2021) Charnay, B., Blain, D., Bézard, B., et al. 2021, A&A, 646, A171
  * Chaverot et al. (2022) Chaverot, G., Turbet, M., Bolmont, E., & Leconte, J. 2022, A&A, 658, A40
  * Clément et al. (2023) Clément, N., Leconte, J., Spiga, A., et al. 2023, Submitted to A&A
  * Cloutier et al. (2017) Cloutier, R., Astudillo-Defru, N., Doyon, R., et al. 2017, A&A, 608, A35
  * Daley-Yates et al. (2021) Daley-Yates, S., Padioleau, T., Tremblin, P., Kestener, P., & Mancip, M. 2021, A&A, 653, A54
  * Ding & Pierrehumbert (2016) Ding, F. & Pierrehumbert, R. T. 2016, ApJ, 822, 24
  * Forget et al. (2006) Forget, F., Montabone, L., & Lebonnois, S. 2006, in Mars Atmosphere Modelling and Observations, ed. F. Forget, M. A. Lopez-Valverde, M. C. Desjean, J. P. Huot, F. Lefevre, S. Lebonnois, S. R. Lewis, E. Millour, P. L. Read, & R. J. Wilson, 422
  * Garaud (2018) Garaud, P. 2018, Annual Review of Fluid Mechanics, 50, 275
  * Guillot (1995) Guillot, T. 1995, Science, 269, 1697
  * Guillot & Showman (2002) Guillot, T. & Showman, A. P. 2002, A&A, 385, 156
  * Habib & Pierrehumbert (2023) Habib, N. & Pierrehumbert, R. T. 2023, arXiv e-prints, arXiv:2310.08202
  * Hu et al. (2021) Hu, R., Damiano, M., Scheucher, M., et al. 2021, ApJ, 921, L8
  * Innes et al. (2023) Innes, H., Tsai, S.-M., & Pierrehumbert, R. T. 2023, arXiv e-prints, arXiv:2304.02698
  * Kasahara (1974) Kasahara, A. 1974, Monthly Weather Review, 102, 509
  * Koll & Cronin (2019) Koll, D. D. B. & Cronin, T. W. 2019, ApJ, 881, 120
  * Kopparapu et al. (2013) Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, ApJ, 765, 131
  * Lane et al. (2003) Lane, T. P., Sharman, R. D., Clark, T. L., & Hsu, H.-M. 2003, Journal of Atmospheric Sciences, 60, 1297
  * Laprise (1992) Laprise, R. 1992, Monthly Weather Review, 120, 197
  * Leconte (2021) Leconte, J. 2021, A&A, 645, A20
  * Leconte et al. (2013) Leconte, J., Forget, F., Charnay, B., et al. 2013, A&A, 554, A69
  * Leconte et al. (2017) Leconte, J., Selsis, F., Hersant, F., & Guillot, T. 2017, A&A, 598, A98
  * Ledoux (1947) Ledoux, P. 1947, ApJ, 105, 305
  * Lefèvre et al. (2018) Lefèvre, M., Lebonnois, S., & Spiga, A. 2018, Journal of Geophysical Research (Planets), 123, 2773
  * Lefèvre et al. (2021) Lefèvre, M., Turbet, M., & Pierrehumbert, R. 2021, ApJ, 913, 101
  * Li & Ingersoll (2015) Li, C. & Ingersoll, A. P. 2015, Nature Geoscience, 8, 398
  * Madhusudhan et al. (2021) Madhusudhan, N., Piette, A. A. A., & Constantinou, S. 2021, ApJ, 918, 1
  * Madhusudhan et al. (2023) Madhusudhan, N., Sarkar, S., Constantinou, S., et al. 2023, ApJ, 956, L13
  * Manabe & Wetherald (1967) Manabe, S. & Wetherald, R. T. 1967, Journal of Atmospheric Sciences, 24, 241
  * Markham et al. (2022) Markham, S., Guillot, T., & Stevenson, D. 2022, A&A, 665, A12
  * Markham & Stevenson (2021) Markham, S. & Stevenson, D. 2021, PSJ, 2, 146
  * Misener & Schlichting (2022) Misener, W. & Schlichting, H. E. 2022, MNRAS, 514, 6025
  * Mlawer et al. (2012) Mlawer, E. J., Payne, V. H., Moncet, J. L., et al. 2012, Philosophical Transactions of the Royal Society of London Series A, 370, 2520
  * Nakajima et al. (2000) Nakajima, K., Takehiro, S.-i., Ishiwatari, M., & Hayashi, Y.-Y. 2000, Geochim. Res. Lett., 27, 3129
  * Parmentier et al. (2013) Parmentier, V., Showman, A. P., & Lian, Y. 2013, A&A, 558, A91
  * Richard et al. (2011) Richard, C., Gordon, I., Rothman, L., et al. 2011, Journal of Quantitative Spectroscopy and Radiative Transfer,
  * Robinson & Catling (2014) Robinson, T. D. & Catling, D. C. 2014, Nature Geoscience, 7, 12
  * Rosenblum et al. (2011) Rosenblum, E., Garaud, P., Traxler, A., & Stellmach, S. 2011, ApJ, 731, 66
  * Schwarzschild & Härm (1958) Schwarzschild, M. & Härm, R. 1958, ApJ, 128, 348
  * Selsis et al. (2023) Selsis, F., Leconte, J., Turbet, M., Chaverot, G., & Bolmont, É. 2023, Nature, 620, 287
  * Skamarock (2006) Skamarock, W. C. 2006, Monthly Weather Review, 134, 2241
  * Skamarock et al. (2019) Skamarock, W. C., Klemp, J. B., Dudhia, J., et al. 2019, NCAR Technical Note NCAR. January 2019. Mesoscale and Microscale Meteorology Division. National Center for Atmospheric Research. Boulder, 1
  * Spiga & Forget (2009) Spiga, A. & Forget, F. 2009, Journal of Geophysical Research (Planets), 114, E02009
  * Stern (1960) Stern, M. E. 1960, Tellus, 12, 172
  * Venot et al. (2020) Venot, O., Cavalié, T., Bounaceur, R., et al. 2020, A&A, 634, A78
  * Venot et al. (2018) Venot, O., Drummond, B., Miguel, Y., et al. 2018, Experimental Astronomy, 46, 101
  * Wordsworth et al. (2011) Wordsworth, R. D., Forget, F., Selsis, F., et al. 2011, ApJ, 733, L48
  * Yang et al. (2016) Yang, J., Leconte, J., Wolf, E. T., et al. 2016, ApJ, 826, 222
  * Zahnle & Marley (2014) Zahnle, K. J. & Marley, M. S. 2014, ApJ, 797, 41

##  Appendix A Calculation of the energy budget

To help visualize the flow of energy in our simulations, we compute the net
upward energy fluxes as follows. In the 3D model, the physical
parametrizations give us at each timestep the specific heating rate due any
diabatic heating process in each cell ($Q_{\mathrm{i}}$). The net upward flux
for process $\mathrm{i}$ at any level defined by the hydrostatic pressure
$\pi$ is defined as the integral of the heating between that level and the
surface

$\displaystyle
F_{\mathrm{i}}=-\int_{\pi}^{\pi^{s}}Q_{\mathrm{i}}\frac{\mathrm{d}\pi}{g},$
(25)

where the minus sign comes from the ordering of the integral boundaries.

The sensible heat transport is computed from the velocity and potential
temperature fields using $F_{\mathrm{dyn}}\equiv c_{p}\langle\rho
w^{\prime}\theta^{\prime}\rangle$, where the primes denote departures from the
average values.