# Convection modeling of pure-steam atmospheres

Xianyu Tan, Maxence Lefèvre and Raymond T. Pierrehumbert  Atmospheric Oceanic
and Planetary Physics, Department of Physics, University of Oxford, OX1 3PU,
United Kingdom<EMAIL_ADDRESS>
###### Abstract

Condensable species in either vapor or condensed form are crucial in shaping
planetary climate. A wide range of planetary climate systems involve
understanding non-dilute condensable substances and their influence on climate
dynamics. There has been progress on large-scale dynamical effects and on 1D
convection parameterization, but resolved 3D moist convection remains
unexplored in non-dilute conditions, though it can have a profound impact on
temperature/humidity profiles and cloud structure. In this work, we tackle
this problem for pure-steam atmospheres using three-dimensional, high-
resolution numerical simulations of convection in post-runaway atmospheres
where the water reservoir at the surface has been exhausted. We show that the
atmosphere is comprised of two characteristic regions, an upper condensing
region dominated by gravity waves and a lower noncondensing region
characterized by convective overturning cells. Velocities in the condensing
region are much smaller than those in the lower noncondensing region, and the
horizontal temperature variation is small ($\lesssim 1$ K) overall.
Condensation in the thermal photosphere is largely driven by radiative cooling
and tends to be statistically homogeneous. Some condensation also happens
deeper, near the boundary of the condensing region, due to triggering by
gravity waves and convective penetrations and exhibit random patchiness. This
qualitative structure is insensitive to varying parameters in the model, but
quantitative details may differ. Our results confirm theoretical expectations
that atmospheres close to the pure-steam limit do not have organized deep
convective plumes in the condensing region. The generalized convective
parameterization scheme discussed in Ding & Pierrehumbert (2016) is
appropriate to handle the basic structure of atmospheres near the pure-steam
limit but is difficult to capture gravity waves and their mixing that appear
in 3D convection-resolving models.

###### Subject headings:

Exoplanet atmospheres; Planetary atmospheres

††slugcomment: Accepted to AAS Journals

##  1\. introduction

Characterization of the smaller exoplanets is beginning to come into view,
(e.g., Kreidberg et al. 2014; Knutson et al. 2014; Benneke et al. 2019;
Tsiaras et al. 2019), and their atmospheres can present greater compositional
diversity than the $\mathrm{H_{2}/He}$ dominated atmospheres of hot Jupiters.
Condensable species, either in vapor or condensed form, are important in
determining planetary climate (Pierrehumbert 2010). The need to better
understand the planetary climate dynamics with non-dilute condensable
constituents is demanding. Non-diluteness refers to a situation wherein the
mass of condensable substances can be comparable to that of the noncondensing
components. Exoplanets with sizes in between Neptune and super-Earths may have
water-rich atmospheres (e.g., Zeng et al. 2019; Mousis et al. 2020; Otegi et
al. 2020; Harman et al. 2021). Terrestrial planets inward of the habitable
zone may experience a runaway greenhouse wherein the ocean evaporates and
water vapor dominates the atmosphere (Kasting et al. 1993). Steam atmospheres
are relevant to terrestrial planets during the magma ocean phase immediately
following accretion (e.g., Zahnle et al. 1988; Hamano et al. 2013), and for
young planets whose water vapor has not yet condensed into an ocean (Turbet et
al. 2021). Extremely close-in rocky planets may have rock-vapor atmospheres on
the dayside which are condensable during the transport to the nightside
(Castan & Menou 2011). Finally, Mars and (marginally) Titan are in this regime
as well with the major condensable constituents being ${\rm CO_{2}}$ and ${\rm
CH_{4}}$, respectively.

Several studies have targeted the climate dynamics with non-dilute condensable
vapor. Ding & Pierrehumbert (2016) examined the energy budget associated with
non-dilute condensation and precipitation, and proposed a scheme for
convection parameterization applicable in general conditions. Pierrehumbert &
Ding (2016) demonstrated the novelty on large-scale dynamics arisen from mass
transport by precipitation and constraints from the condensation
thermodynamics, then presented general circulation models (GCMs) in non-dilute
conditions. Yamashita et al. (2016) performed two-dimensional convection
modeling for a pure-${\rm CO_{2}}$ atmosphere in Martian conditions. Ding &
Pierrehumbert (2018) further demonstrated the importance of horizontal heat
transport in determining global surface temperature variation of pure-steam
atmospheres. Ding & Pierrehumbert (2020) showed that phase-curve information
may be used to distinguish richness of water vapor in atmospheres of slowly-
rotating, tidally-locked terrestrial planets. Turbet et al. (2021) performed
long-term GCM simulations for early hot steam atmospheres of Earth and Venus
and suggested the important role of cloud radiation effects.

Understanding the nature of convection is vital as it is an important form of
energy and mass transport in the atmospheres and can profoundly impact
climates. Three dimensional resolved convection simulations are a valuable
tool for advancing the understanding of convection, and for testing the
parameterizations that are essential for representing convection in most
global circulation models. Convection modeling has been carried out in the
context of exoplanets (Zhang et al. 2017; Sergeev et al. 2020; Lefèvre et al.
2021; Song et al. 2021), although mostly in Earth-like conditions. Convection
in non-dilute conditions is in a novel regime that is little explored, though
work on 1D hydrostatic parameterizations provides some hypotheses as to the
expected behavior (Ding & Pierrehumbert 2016). 1D radiative-convective models
of pure-steam atmospheres typically assume a temperature structure of a dry
adiabat attached to the surface temperature and a dew-point adiabat of water
vapor within the saturated region (Goldblatt et al. 2013; Hamano et al. 2013;
Boukrouche et al. 2021). This assumption needs to be validated using
convection models that self-consistently capture relevant dynamical and
condensation processes.

In this work, we investigate the 3D nature of convection for the limiting case
of a pure-steam atmosphere. This is easier to understand than general
nondilute conditions and sets a baseline for future work that incorporates
varying fractions of non-condensable gases, but there are also numerous
planetary phenomena for which the pure-steam limit is relevant. In Section 2,
we introduce our numerical model; in Section 3 we present our basic results
and sensitivity exploration; and finally we discuss and conclude in Section 4.

##  2\. Model

###  2.1. CM1

We utilize the Cloud Model 1 (CM1), a 3D, non-hydrostatic model that has been
widely used for convection studies (Bryan & Fritsch 2002, see also
https://www2.mmm.ucar.edu/people/bryan/cm1/), to perform high-resolution,
idealized experiments of convection over a regional domain. We use a plane-
parallel, two-stream radiative transfer scheme with a grey approximation in
the thermal emission. The numerical tool TWOSTR (Kylling et al. 1995) is
employed to solve the radiative transfer equations and is coupled to the
dynamics of CM1. The gas opacity of water vapor is assumed to be $0.1\;{\rm
m^{2}kg^{-1}}$, same as that used in Ding & Pierrehumbert (2016). We apply a
fixed surface temperature as a lower boundary condition and treat it as a free
parameter. Similar lower boundary conditions have been widely used in modeling
studies of convective aggregation (see a recent review by Wing et al. 2017) as
well as 1D radiative-convective models for runaway climate (Boukrouche et al.
2021). We assume a solid surface, and there is no evaporative surface flux to
the atmosphere. This is a post-runaway situation, after the ocean has
completely evaporated into the atmosphere leaving a subsaturated layer near
the ground. The atmosphere is assumed to be transparent to instellation. Note
that water vapour has strong near-IR absorption, which allows the incoming
stellar radiation to heat the atmosphere. In this paper, the effect is
neglected so as to focus on the essentials of the problem. We utilize the
fully-compressible equation set with only water vapor in our models, but with
a scheme to deal with condensation, rainout and evaporation as described
below. This study is the first extraterrestrial use of the CM1.

###  2.2. Condensation, rainout and evaporation

In Ding & Pierrehumbert (2016) the First Law of Thermodynamics was used to
adjust an entire atmospheric column to an energetically consistent state
neutrally stable to convection, after buoyancy-generated kinetic energy has
been dissipated as heat. Here, because we use a nonhydrostatic model that
explicitly resolves the conversion of potential to kinetic energy by buoyancy,
the mixing it causes, its subsequent dissipation as heat,111 In our 3D
convection resolving model, the sources of kinetic energy dissipation include
surface stress from the lower boundary layer, parameterized subgrid turbulent
dissipation, and numerical dissipation. The dissipation of kinetic energy
turns into a heat source in the thermodynamics equation, and this is a built-
in option in the CM1 model.  and the pressure adjustment triggered by the
removal of precipitation, our only use of the First Law is to determine the
amount of condensate produced or evaporated within individual grid cells, and
the effects on local pressure and temperature. An important difference between
our 3D modeling and the 1D parameterization scheme (Ding & Pierrehumbert 2016)
is that the adjustment of the column is handled explicitly by the dynamical
core of our 3D model. The adjustment occurring in our model – which is not
required to be complete – is a natural emergent property of the dynamics.

A two step process is modeled, in which precipitation is first produced in
situ and then removed by rainout. Let’s suppose that we start with a mass of
atmosphere with no condensate present, but that radiative cooling or adiabatic
ascent creates some supersaturation; a small amount of condensate will form,
the pressure will adjust to the phase boundary, and the latent heat release
will slightly increase the temperature of the parcel. Without energy and mass
exchange with the environment during condensation, the air parcel follows the
following relation based upon the first law of thermodynamics and mass
conservation (Emanuel 1994; Ding & Pierrehumbert 2016):

$d\left(k-\frac{p}{\rho}\right)+pd\frac{1}{\rho}=0,$ (1)

where $\rho$ is the atmospheric density which includes all gases and the
condensates, $k=c_{pl}q_{t}T+Lq_{c}$, $c_{pl}$ is that for the condensates,
$L$ is the latent heat, $q_{t}$ is the mass concentration of the total
condensable in both phases, and $q_{c}$ is the mass concentration of the
condensable vapor. $k$ is the moist enthalpy in the pure steam limit, in which
we also have $q_{t}=1$. Because CM1 is a nonhydrostatic model using $z$ rather
than $p$ as vertical coordinate, we have written the First Law in constant-
volume rather than constant-pressure form.

For condensation occurring at a grid point in the model without loss of mass,
the term $pd\frac{1}{\rho}$ in Eq.(1) drops out, and there is a conserved
quantity $\rho k-p$ before and after condensation. In a pure-steam atmosphere,
this quantity (denoted by $\mathcal{A}$) is
$\mathcal{A}\equiv\rho(c_{pl}T+Lq_{c})-p$. We make use of this quantity to
determine post-condensation temperature and pressure and the amount of
condensates. We start with a pre-condensation situation wherein there is only
gas with a temperature $T_{1}$ and pressure $p_{1}$, and the gas is
supersaturated due to radiative or adiabatic cooling. Then

$\mathcal{A}=\rho_{1}(c_{pl}T_{1}+L)-p_{1}=\rho_{1}c_{pl}T_{1}+(\frac{L}{RT_{1}}-1)p_{1},$
(2)

where $\rho_{1}=\frac{p_{1}}{RT_{1}}$. After condensation, condensates form
and temperature and pressure are adjusted, which are denoted as $T_{2}$ and
$p_{2}$, respectively. At this point, the condensates are retained in the
atmosphere. The temperature of the condensates is also $T_{2}$. $T_{2}$ and
$p_{2}$ are restricted to the phase boundary $p_{sat}(T)$, given by the
Clausius-Clayeyron relation, because of the coexistence of condensed and
vapour phases. The conserved quantity in the post-condensation is then

$\begin{split}\mathcal{A}&\equiv\rho_{2}(c_{pl}T_{2}+Lq_{c})-p_{2}\\\
&=\rho_{1}c_{pl}T_{2}+L\rho_{c}-p_{2}\\\
&=\rho_{1}c_{pl}T_{2}+\left(\frac{L}{RT_{2}}-1\right)p_{sat}(T_{2}),\end{split}$
(3)

in which we made use of conservation of mass $(\rho_{1}=\rho_{2})$. With
$\mathcal{A}$ given by the initial values using Eq. (2), we solve for $T_{2}$
in Eq. (3) using Newton’s method. Because the range of temperatures
encountered within the condensing layer is not large, we found it sufficient
to use the analytic constant-$L$ form of $p_{sat}$, though it is an assumption
that would be easy to relax. The condensate density $\rho_{\rm cond}$ is
obtained from the difference of the gas density before and after the
condensation.

With suitable cloud condensation nuclei, condensates form with a small
supersaturation ratio $\mathcal{S}$ to overcome surface tension (Houze 2014),
in which $1+\mathcal{S}=p_{1}/p_{\rm sat}$. We allow $\mathcal{S}$ to be
nonzero but still much smaller than 1. It is assumed to be $10^{-7}$ in most
simulations to represent the no-surface-barrier limit but is varied up to 0.2
in some cases for sensitivity exploration. Note that the post-condensation
pressure is adjusted to $p_{\rm sat}$.

We do not model the detailed processes of cloud condensational growth and
coagulation. Instead, the condensates are assumed to instantaneously fall out
after condensation occurs and evaporate in the subsaturated regions.
Evaporation is likewise treated using the conserved quantity $\mathcal{A}$.
Pre-evaporation temperature and pressure in the evaporation regions are
denoted as $T_{3}$ and $p_{3}$, and the total density is
$\rho_{3}=\frac{p_{3}}{RT_{3}}+\rho_{\rm cond}$, where $\rho_{\rm cond}$ is
the density of the soon-to-be-evaporated precipitation. The conserved quantity
is $\mathcal{A}=\rho_{3}(c_{pl}T_{3}+Lq_{c})-p_{3}$, where
$q_{c}=\frac{p_{3}}{RT_{3}}/\rho_{3}$. After evaporation, temperature $T_{4}$
and pressure $p_{4}$ are finally obtained via

$\begin{split}\mathcal{A}&=\rho_{3}(c_{pl}T_{4}+L)-p_{4}\\\
&=\rho_{3}T_{4}(c_{pl}-R)+\rho_{3}L.\end{split}$ (4)

We adopt a highly idealized approach to determine the location of evaporation.
All precipitation evaporates within the atmosphere. Evaporation occurs only
when the relative humidity is below 90%. The evaporative mass is assumed to
evenly distribute over a fixed accumulative depth (this depth needs not be
continuous), and this depth is treated as a free parameter. Our fiducial
models assume an evaporative depth of 10 km, but we will investigate the
sensitivity of results to the evaporative depth.

There are two additional considerations to conserve energy associated with
rainout. The first is the dissipation of gravitational energy of the falling
condensates. We assumed that this is dissipated via friction, and the
frictional heating goes into the gas instantaneously. The second is to
consider the temperature difference between the condensates and the air along
the falling path. This heat exchange is assumed to occur via conduction and
reaches equilibrium instantaneously.

###  2.3. Numerical setup and convergence

We use the large-eddy-simulation setup that integrates the filtered Navier-
Stokes equations. Stresses and fluxes on the resolved flow by sub-grid
turbulence are parameterized using a prognostic turbulence kinetic energy
scheme similar to that used in Deardorff (1980). Acoustic waves are treated
explicitly in both horizontal and vertical directions using a time-splitting
technique. The horizontal boundaries are periodic and the vertical boundaries
are impermeable. A Rayleigh damping is applied to winds in the top 10 km
layers, and surface stress and heat flux are calculated using the original CM1
formulation at the bottom boundary. We applied a 6th-order numerical diffusion
to maintain numerical stability. The model domain is in Cartesian geometry
with a typical grid space of 600 m in the horizontal directions. The vertical
grid space $\delta z$ ranges between 600 and 680 m depending on the surface
temperature. Our canonical models have $320\times 320\times 250$ (in $x$, $y$
and $z$) grid points. The resolution is chosen to better resolve convection
and waves while remaining computationally feasible. Following the default
values in CM1, we adopt the following constants: surface gravity $g=9.81\;{\rm
m\;s^{-2}}$, the heat capacity at constant pressure for water vapor
$c_{p}=1870\;{\rm J~{}kg^{-1}~{}K^{-1}}$ and for condensates
$c_{pl}=4190\;{\rm J~{}kg^{-1}~{}K^{-1}}$, the specific gas constant
$R=461.9\;{\rm J~{}kg^{-1}~{}K^{-1}}$, and latent heat $L=2.5\times
10^{6}\;{\rm J\;kg^{-1}}$. The models are initialized with a dry adiabat
attached to the surface and a moist adiabat when gas is saturated. Most
simulations assume a 1-bar surface pressure which allows exploration of the
essential dynamical features of the problem without the computational expense
of a thicker atmosphere. The setup corresponds to the behaviour of a planet
which started with a low mass ocean; we can speculate that it would
qualitatively mimic the behaviour of the upper bar of a much deeper atmosphere
with a thicker noncondensing region. This setup treats the post-runaway
climate that is interesting to a branch of 1D climate studies (e.g., Kasting
et al. 1993; Goldblatt et al. 2013; Hamano et al. 2013; Boukrouche et al.
2021), whereas the parametereization in Ding & Pierrehumbert (2016) did not
address this geometry.

We ran a model with a lower resolution of 1 km in all directions up to 88
simulation days. The total kinetic energy reaches a statistical equilibrium
after only a few model days. The total internal and potential energy reach a
statistical equilibrium after about 25 days and then slightly oscillate around
the mean. However, quantitative properties of convective flows and gravity
waves are insensitive to the long-term convergence of the model. For our
canonical runs with higher resolution, we integrate the system only up to 5 to
11 simulation days as they are computationally more costly, and statistical
results are obtained by outputs of the last 6 hours.

Figure 1.— Results from CM1 simulations with three surface temperature of 600
K, 700 K and 800 K, surface pressure of 1 bar and supersaturation ratio
$\mathcal{S}=10^{-7}$, showing basic structure that is characterized by an
upper condensing region and a lower dry convecting region. Statistical results
are averaged over the last 6 hours of integration. _Panel (a):_ randomly
selected instantaneous temperature-pressure profiles (black lines). Note that
there are multiple black lines per case but they have very small horizontal
temperature differences and so they appear to be merged. Thick grey lines are
dry adiabats corresponding to different surface temperatures; the thick green
line is the moist adiabat for pure-water atmosphere. _Panel (b):_ horizontal-
and time-averaged root mean square (RMS) temperature variations as a function
of mean pressure. _Panel (c):_ horizontal- and time-averaged RMS vertical
velocities as a function of mean pressure. _Panel (d):_ horizontal- and time-
mean condensation mass rate (positive meaning condensation and negative
meaning evaporation) as a function of mean pressure.

##  3\. Results

###  3.1. Basic structure and dynamics

Our simulations show that the atmospheric domain is generally comprised of two
characteristic regions — a relatively quiescent, stratified upper condensing
region and a vigorously convecting lower dry region, separated by the first
condensation level. This is similar to those found in 2D convection
simulations of Yamashita et al. (2016). We start with describing horizontal-
mean properties of the simulations. Figure 1 shows results as a function of
mean pressure from three models with a surface pressure of 1 bar,
supersaturation ratio $\mathcal{S}=10^{-7}$ and three surface temperatures of
600, 700 and 800 K. In both condensing and dry regions, the horizontal root-
mean-square (RMS) temperature variations are typically much smaller than 1 K
(panel b) in the domain except near the first condensation level where
convective overshoots occur. The instantaneous temperature-pressure profiles
in panel (a) of Figure 1 appear to be merged to a single dry adiabat in the
lower region and to the moist adiabat in the upper region. The RMS vertical
velocities in panel (c) reach several to more than 10 $\;{\rm m~{}s^{-1}}$ in
the lower dry region but rapidly decrease above the first condensation level.
The overall vertical speed increases with increasing surface temperature. The
domain-mean condensation/evaporation mass rates in panel (d) show evaporation
below the first condensation level but a sharp transition to condensation
right above the first condensation level. The latter is driven by the
convective overshoots and gravity waves. Low pressures ($\lesssim 1$ mbar)
show a vertically broad and smooth condensing region that is primarily driven
by radiative cooling.

The upper condensing zone further exhibits two subdivisions. One is above the
first condensation level but below the radiative cooling zone. This region is
optically thick and has small radiative cooling rates (on the order of
$10^{-5}\;{\rm Ks^{-1}}$). While the ascending branch of wave motions
supersaturates and is adjusted to saturation, the descending branch is free of
thermal damping. Waves can vertically propagate through this region and exert
certain vertical velocity and temperature fluctuations. The other subdivision
at the thermal photosphere ($\lesssim 1$ mbar) show negligible RMS temperature
and vertical velocity. Radiatively driven condensation acts uniformly in the
horizontal direction and sets stringent constraints on the local condition.
The collapse of pure-steam gas into a unique property determined by the
Clausius-Clapeyron relation eliminates dynamical perturbations (Pierrehumbert
& Ding 2016).

Figure 2.— Instantaneous snapshots of the model with surface temperature of
800 K, surface pressure of 1 bar and supersaturation ratio
$\mathcal{S}=10^{-7}$. _Panel (a):_ Vertical velocity at a slice as a function
of $y$ and pressure. _Panel (b) to (d)_ : Vertical velocities at different
pressure levels, with (b) representing those at the lower dry convective zone,
(c) at the overshooting zone slightly above the first condensation level, and
(d) at the radiatively cooling zone. _Panel (e) and (f):_ condensation mass
rate at different pressure levels with (e) in the overshooting zone and (f) in
the radiative cooling zone.

Note that since the condensation scheme is applied at each time step, and
condensate is not allowed to accumulate, at most a small portion of the vapour
is converted to condensate and subsequently removed by precipitation. There is
never a large proportion of atmospheric mass removed (see the small
condensation mass rate in panel (d) of Figure 1); the interesting aspect of
diluteness is that the small mass loss at each condensation step can add up
over time to transport a significant proportion of atmospheric mass. The
condensed phase occupies very little volume, so its removal does not
significantly change the $(p,T)$ of the gas phase. In a non-hydrostatic model,
removal of precipitation pushes the atmospheric column out of hydrostatic
balance, but the dynamics handles the resulting pressure adjustment
explicitly. If the expansion caused by unburdening lower air parcels causes
cooling and supersaturation, that will be handled by condensation at
subsequent time steps.

Now we present the 3D structure of convection and condensation. Dynamics in
the lower noncondensing region is characterized by a major convective
overturning cell that spans across the whole horizontal domain and with peak
velocities reaching several tens of $\;{\rm m~{}s^{-1}}$. Negative buoyancy is
generated in the upper parts of the noncondensing region by evaporative
cooling, resulting in denser currents that penetrate deep down. Figure 2
displaces snapshots of vertical velocity in a vertical cross section and
several horizontal cross sections, as well as condensation mass rates in two
cross sections from the model with a surface temperature of 800 K. Panel (a)
and (b) show a broad and coherent upwelling cell and narrower ridges of
downwelling. As a consequence of the area asymmetry and continuity, the
downdrafts have larger speeds than the updrafts.

In the condensing region, there is no deep penetrating plume even though
latent heat is released, which is in stark contrast to the deep moist
convection on Earth. The lack of buoyancy generation in the pure-steam
atmosphere is because once the atmosphere is saturated, the temperature-
pressure profiles collapse into a single adiabat that is uniquely determined
by the Clausius–Clapeyron relation. The condensing air parcel is neutrally
stable with regard to the environment because the air parcel shares the same
density-pressure trajectory as the environment (Colaprete & Toon 2003; Ding &
Pierrehumbert 2016). On the other hand, the condensing region is strongly
stratified against downward motions that are not associated with condensation.
While the saturated air may still have the possibility to move freely upward
as it is neutrally buoyant, its motions will be limited by the stratification
experienced by downward-moving air via mass continuity. This excludes
convective instability in the saturated region even with vigorous
perturbations from the lower noncondensing layer. Behaviour of this sort was
hypothesized in the 1D treatment of (Ding & Pierrehumbert 2016), but the
nonhydrostatic model is free to deviate from it; the confirmation of the
supposition in a fully dynamic model is thus significant.

In our simulations all the condensation and cloud formation is generated by
small vertical displacements by gravity waves or by radiatively driven
condensation in the thermal photosphere. Internal gravity waves are triggered
by the dry convective penetration and show much smaller spatial structure and
velocity magnitudes than those in the dry convective region. Panel (c) and (d)
show characteristic wave patterns. Gravity waves generated from convective
perturbation have long been identified in Earth and planetary atmospheric
modeling (e.g., Fovell et al. 1992; Baker et al. 2000; Lefèvre et al. 2018).
Time evolution of these fields show a good correlation between locations of
upwelling convective plume and sources of gravity waves. Gravity-wave induced
temperature variations can either trigger condensation in the ascending
regions or permit evaporation in the descending regions. Panel (e) illustrates
that precipitation fall from higher altitudes can evaporate at locations with
downward velocity (and thereby are hotter and subsaturated), whereas regions
with upward velocity can generate condensation. In the radiative zone shown in
panel (f), only condensation can occur and it is more widespread. Its spatial
pattern generally correlates well with the gravity waves.

###  3.2. Sensitivity exploration

We perform experiments with varying parameters to evaluate the qualitative and
quantitative differences of models to those shown in Section 3.1. These
experiments have the same parameters and setup as the one with a surface
pressure of 1 bar, a surface temperature of 800 K, a supersaturation ratio
$\mathcal{S}=10^{-7}$, evaporation depth of 10 km and no rotation, except
those specified in the following text. We show that the basic structure and
dynamics remain qualitatively the same, but the quantitative dynamical details
can differ in some cases.

We first conduct additional experiments with different supersaturation ratio
$\mathcal{S}=0.01$, 0.1 and 0.2. In the pure-steam condition, the
characteristic permitted temperature variation before and after condensation
is $\delta T\approx\mathcal{S}T^{2}R/L$ when $\mathcal{S}\ll 1$. Higher
$\mathcal{S}$ permits more flexibility of temperature variation in the
condensing region, and potentially stronger waves there. The first row in
Figure 3 shows results of the experiments. Indeed, in the upper radiative
cooling zone ($\lesssim 1$ mbar), higher $\mathcal{S}$ leads to larger RMS
vertical velocity (left panel) and larger RMS temperature variation ($\gtrsim
0.2$ K for $\mathcal{S}\geq 0.1$, not shown). Somewhat surprisingly, the RMS
vertical velocity at the region below the radiative cooling zone and above the
dry convective zone remain invariant with different $\mathcal{S}$. In the
lower dry region, the RMS vertical velocity slightly decreases with increasing
$\mathcal{S}$, although the decrement appears to saturate when
$\mathcal{S}\gtrapprox 0.1$. This results in weaker overshoot-driven
condensation when $\mathcal{S}$ is larger as shown on the right.

At a given time, condensation in the radiative cooling region is widespread in
the case with $\mathcal{S}=10^{-7}$ but is highly sparse when $\mathcal{S}$ is
no longer tiny (not shown). Dynamically, when a certain region condenses and
adjusts back to a saturated state, its perturbation rapidly propagates. This
maintains other regions “warm” such that condensation does not easily occur.
Energetically, heat released from a single condensing event increases with
increasing $\mathcal{S}$. Given the same radiative cooling rates, area
fraction of condensation decreases with increasing $\mathcal{S}$ to balance
radiative cooling. Interestingly, the case with $\mathcal{S}=10^{-7}$ shows a
larger domain-mean condensation rate at the photosphere that is almost
balanced by radiative cooling, while others with $\mathcal{S}\gtrapprox 0.01$
show smaller condensation rates. This implies that gas dynamics becomes
increasingly important in energy transport with a small but non-negligible
supersaturation ratio.

Next, we present experiments with varying location and depth of the
evaporation shown in the second row of Figure 3. These affect the the
generation of negative buoyancy in the dry convective region because
evaporation is the major cooling mechanism to drive the dry convection. In the
first experiment, we extend the evaporation depth to 45 km as appose to the
original 10 km; and in the second case, evaporation is only permitted in the
lower 10 km above the ground. These settings can be visualized by the
evaporative mass rate profiles shown on the right column. In both cases,
although the basic structure and characteristic dynamics remain the same, the
magnitude of the velocities is smaller than the canonical case. This is not
surprising as negative buoyancy are expected to be smaller in both cases. The
case with an evaporation depth of 45 km shows the weakest activity, with RMS
vertical velocity $<3\;{\rm m~{}s^{-1}}$. The case with a lower 10 km
evaporation depth is in between. The RMS vertical velocities above the first
condensation level of both experiments are negligible compared to the nominal
case, which is related to the weaker penetrative plume at the top of the dry
convective zone and therefore a weaker generation of gravity waves. As a
result, both cases show little overshoot-driven condensation. In terms of the
spatial organization, the dry convective cell spans the whole horizontal
domain in all cases (not shown). But the updrafts and downdrafts in the case
with an evaporation depth of 45 km appear to be more symmetric than the other
two cases.

Figure 3.—  Horizontal- and time-averaged RMS vertical velocity (left column)
and condensation mass rate (right column) as a function of pressure for
various models. _First row:_ models with surface temperature of 800 K and
surface pressure of 1 bar, but with four supersaturation ratio
$\mathcal{S}=10^{-7},~{}0.01,~{}0.1$ and 0.2. _Second row:_ models with
different evaporation height and location indicated by the thicknesses and
locations of the evaporative mass rates on the right. _Third row:_ models with
different surface pressures but with the same dry adiabat as that with a
surface temperature of 800 K and surface pressure of 1 bar. _Bottom row:_
models with zero and two different Coriolis parameter $f=1.45\times 10^{-4}$
and $1.45\times 10^{-3}\;{\rm s^{-1}}$, and a supersaturation ratio
$\mathcal{S}=0.1$.

We then carry out experiments with a varying surface pressure of 0.5 bar and 2
bars adjusting the surface temperature so as to keep the same dry adiabatic
profile as that of the nominal case. The similarity of both the RMS vertical
velocities and condensation mass rates shown in the third row of Figure 3
demonstrates that this level of variation on the surface pressure changes the
quantitative dynamics little. This is a bit surprising as we would have
expected that the higher the surface pressure, the more vigorous the dry
convection because more vertical length would facilitate larger convective
available potential energy for a nondiluted convective plume. Perhaps the dry
convection is sufficiently turbulent (see Figure 2) that such undiluted
convective plumes are not present. This experiment also suggests that the
behavior near and in the condensing layer is not too sensitive to processes
occurring near the surface or in the noncondensing deeper layers.

Lastly, we examine the role of rotation in shaping the convection and whether
coherent vortex would form in small-scale nonhydrostatic models in pure-steam
conditions. One might first expect vigorous hurricanes to form because the
atmosphere already contains enormous quantities of latent heat, whereas for
Earth-like hurricanes the latent heat has to be imported by evaporation from
the ocean surface (Emanuel 2018). 3D radiative-convective equilibrium modeling
applied in Earth’s tropical conditions suggest that rotation triggers
hurricane-like vortices (e.g., Bretherton et al. 2005; Held & Zhao 2008;
Merlis et al. 2016). However, GCMs with rich condensable vapor and Earth’s
rotation rate did not generate hurricane-like features (Pierrehumbert & Ding
2016). It is unclear whether that is due to the hydrostatic nature of GCMs, or
low horizontal resolution, or an intrinsic dynamical property of condensable-
rich atmospheres. Here, we explore two Coriolis parameters, $f=1.45\times
10^{-4}\;{\rm s^{-1}}$ (the value at Earth’s North pole) and another one 10
times of that $f=1.45\times 10^{-3}\;{\rm s^{-1}}$. The domain-mean properties
are shown in the bottom row of Figure 3. All three models have a
supersaturation ratio $\mathcal{S}=0.1$ because we expect that models with
higher flexibility of temperature have a higher chance of being sculpted by
the rotation. The inclusion of rotation results in small changes in the RMS
vertical velocity and condensation rate, and a slight change in the vertical
profile of the RMS vertical velocity in the dry convective zone. Similarly,
differences in the condensation mass rates are small.

We do not find existence of long-lasting, coherent vortices. Figure 4 shows
snapshots of vertical velocity at different pressure levels from the model
with $f=1.45\times 10^{-3}\;{\rm s^{-1}}$. The strong rotation only tends to
limit the horizontal size of dry convective structures. Rotation could result
in swirling structures in the dry convection (for instance, the feature at
around x=100 km and y=50 km in the upper panel). Wave properties in the upper
condensing zone are also influenced by rotation. Perhaps the lack of vortex
formation is not surprising because a systematic baroclinic dynamical
structure is likely needed to maintain the available potential energy against
dissipation due to the surface drag and related vortex dissipation mechanism
(such as the Ekman pumping). Such a baroclinic structure is difficult to
maintain in the pure-steam atmospheres due to the strong constraint of the
Clausius-Clapeyron relation (Colaprete & Toon 2003; Ding & Pierrehumbert 2016;
Pierrehumbert & Ding 2016).

Figure 4.— Snapshots of vertical velocities at the dry convective zone (upper
panel), the overshooting zone (middle panel) and the radiative cooling zone
(bottom panel) from the model with surface temperature 800 K, surface pressure
of 1 bar and supersaturation ratio $\mathcal{S}=0.1$, and including rotation
with $f=1.45\times 10^{-3}\;{\rm s^{-1}}$. Pressures are indicated above each
panel.

##  4\. Discussion

The lack of convective instability and organized convective structure in the
condensing region of pure-steam atmospheres suggests that convection
parameterization in both relevant GCMs and 1D models may safely neglect the
role of deep penetrating plumes and mainly consider condensation,
precipitation and evaporation. The basic parameterization scheme proposed by
Ding & Pierrehumbert (2016) and its application in Pierrehumbert & Ding (2016)
should be sufficient for atmospheres closed to the pure-steam limit. The
domain-mean temperature structure of our full 3D simulations maintains the dry
and moist adiabats as those typically assumed for 1D models, which is
encouraging for 1D radiative-convective models for post-runaway climate
calculations using sophisticated radiative transfer (e.g., Goldblatt et al.
2013; Kopparapu et al. 2013; Boukrouche et al. 2021). GCM simulations of steam
atmospheres in the early stages of Earth and Venus before water vapor has
condensed into an ocean illustrates the importance of understanding climate
dynamics of such atmospheres Turbet et al. (2021). Our configuration is
similar to that in these GCM calculations except for the lack of the
suppression of condensation and dry convection by shortwave heating in the
intermediate layer. This intermediate stratified layer serves to suppress the
strong convective perturbations on the bottom of the saturated layer and
eliminates additional thick cloud formation near there due to these
overshooting, which supports results in Turbet et al. (2021).

Although our models do not include cloud radiative forcing, our results
provide implications for cloud parameterization and observatonal implications
for atmospheres closed to the pure-steam limit. Condensation is driven by
radiative cooling in the thermal photosphere, though some does happen deeper,
near the boundary of the condensing region, due to triggering by gravity
waves. In the photosphere, radiative cooling acts uniformly in the horizontal
direction, which implies that clouds tend to be statistically homogenized at
the photosphere of steam atmospheres. We do not see evidence for aggregation
of condensing patches in the thermal photosphere. The lower condensing region
below the photosphere is characterized by convective penetration and gravity
waves. The resulting cloud structure would be spatially correlated to the
convective plumes and waves. Although these cloud structures cloud would be
highly time variable and patchy, the horizontal mean cloud mixing ratio can be
comparable to those driven by radiative cooling (see the condensation rate in
panel d of Figure 1). Such possible cloud configuration —most clouds in a
steam atmosphere would form in the thermal photosphere, though they can
nonetheless exhibit random patchiness (see Figure 2 d and f)— should have
observational implications for observations of steam atmospheres of both sub
Neptunes and terrestrial exoplanets in the runaway stage. Of course, the
amount of small cloud droplets that are able to keep aloft is sensitive to
microphysics, which in turn is affected by the mode of vertical transport.
This aspect can be quantified only with fully coupled dynamics-microphysics
models. Finally, the suppression of deep convection in the condensing layer
implies weak transport of trace species from the deep atmosphere, which would
have interesting implications for the chemistry of such atmospheres.

To summarize, in this study, we have performed high-resolution, non-
hydrostatic simulations to simulate convection in pure-steam atmospheres. We
find that the atmosphere is characterized by an upper condensing region and a
lower dry convecting region. The condensing region is stratified and
characterized by gravity waves that are triggered by convective overshoots.
The lower dry region is characterized by convective overturning cells that
span over the model domain. Magnitude of velocities in the condensing region
is much smaller than that in the lower dry region. Horizontal temperature
variation is small ($\lesssim 1$ K) overall. Condensation in the thermal
photosphere is driven by radiative cooling and tends to be statistically
homogeneous. Near the boundary of the condensing region, condensation is
triggered by gravity waves and convective penetrations and exhibit random
patchiness. Models with different parameters show a similar qualitative
picture but can be quantitatively different in some cases. Our results should
also be applicable in exoplanets with thick gaseous envelopes filled with
mostly condensable species.

Future extensions would be to include various fractions of non-condensable
gases and to explore convection from dilute to non-dilute conditions with
potentially deeper atmospheres. Effects of clouds, including mass loading and
radiative effects in both visible and thermal bands, may yield more complexity
and feedbacks in controlling the convective systems. Realistic cloud
microphysics schemes might yield a time delay of the full relaxation back to
the phase equilibrium. This could affect the amount of condensates in the air
as well as the heating or cooling rates and therefore influence the generation
of buoyancy on the convective system. If water vapor absorption of
instellation is included, it would reduce the net radiative cooling rate in
the condensing layer, and the convective mass transport in the subsaturated
noncondensing layer (e.g., Turbet et al. 2021). This aspect would be
especially interesting for planets around M dwarfs whose peak spectral energy
is close to near-IR. Gravity waves in our current model are obtained without
large-scale vertical wind shear. Its existence could impact the generation and
vertical propagation of the waves as well as cloud formation (Lefèvre et al.
2020), which is worth examining in future studies. Finally, the suppression of
convection in the condensing layer implies weak transport of trace species
from the deep atmosphere, which may have interesting implications for the
chemistry of such atmospheres.

This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation program (grant
agreement No. 740963/EXOCONDENSE).

## References

  * Baker et al. (2000) Baker, R. D., Schubert, G., & Jones, P. W. 2000, Journal of the atmospheric sciences, 57, 184
  * Benneke et al. (2019) Benneke, B., Knutson, H. A., Lothringer, J., Crossfield, I. J., Moses, J. I., Morley, C., Kreidberg, L., Fulton, B. J., Dragomir, D., Howard, A. W., et al. 2019, Nature Astronomy, 3, 813
  * Boukrouche et al. (2021) Boukrouche, R., Lichtenberg, T., & Pierrehumbert, R. T. 2021, arXiv e-prints, arXiv:2107.14150
  * Bretherton et al. (2005) Bretherton, C. S., Blossey, P. N., & Khairoutdinov, M. 2005, Journal of the atmospheric sciences, 62, 4273
  * Bryan & Fritsch (2002) Bryan, G. H. & Fritsch, J. M. 2002, Monthly Weather Review, 130, 2917
  * Castan & Menou (2011) Castan, T. & Menou, K. 2011, The Astrophysical Journal Letters, 743, L36
  * Colaprete & Toon (2003) Colaprete, A. & Toon, O. B. 2003, Journal of Geophysical Research: Planets, 108
  * Deardorff (1980) Deardorff, J. W. 1980, Boundary-Layer Meteorology, 18, 495
  * Ding & Pierrehumbert (2016) Ding, F. & Pierrehumbert, R. T. 2016, ApJ, 822, 24
  * Ding & Pierrehumbert (2018) Ding, F. & Pierrehumbert, R. T. 2018, The Astrophysical Journal, 867, 54
  * Ding & Pierrehumbert (2020) —. 2020, The Astrophysical Journal Letters, 901, L33
  * Emanuel (2018) Emanuel, K. 2018, Meteorological Monographs, 59, 15
  * Emanuel (1994) Emanuel, K. A. 1994, Atmospheric convection (Oxford University Press)
  * Fovell et al. (1992) Fovell, R., Durran, D., & Holton, J. 1992, Journal of Atmospheric Sciences, 49, 1427
  * Goldblatt et al. (2013) Goldblatt, C., Robinson, T. D., Zahnle, K. J., & Crisp, D. 2013, Nature Geoscience, 6, 661
  * Hamano et al. (2013) Hamano, K., Abe, Y., & Genda, H. 2013, Nature, 497, 607
  * Harman et al. (2021) Harman, C., Kopparapu, R. K., Stefánsson, G., Lin, A. S., Mahadevan, S., Hedges, C., & Batalha, N. E. 2021, arXiv preprint arXiv:2109.10838
  * Held & Zhao (2008) Held, I. M. & Zhao, M. 2008, Journal of the atmospheric sciences, 65, 2003
  * Houze (2014) Houze, R. A. 2014, Cloud dynamics, Vol. 104 (Academic press)
  * Kasting et al. (1993) Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108
  * Knutson et al. (2014) Knutson, H. A., Dragomir, D., Kreidberg, L., Kempton, E. M.-R., McCullough, P., Fortney, J. J., Bean, J. L., Gillon, M., Homeier, D., & Howard, A. W. 2014, The Astrophysical Journal, 794, 155
  * Kopparapu et al. (2013) Kopparapu, R. K., Ramirez, R., Kasting, J. F., Eymet, V., Robinson, T. D., Mahadevan, S., Terrien, R. C., Domagal-Goldman, S., Meadows, V., & Deshpande, R. 2013, The Astrophysical Journal, 765, 131
  * Kreidberg et al. (2014) Kreidberg, L., Bean, J. L., Désert, J.-M., Benneke, B., Deming, D., Stevenson, K. B., Seager, S., Berta-Thompson, Z., Seifahrt, A., & Homeier, D. 2014, Nature, 505, 69
  * Kylling et al. (1995) Kylling, A., Stamnes, K., & Tsay, S.-C. 1995, Journal of Atmospheric Chemistry, 21, 115
  * Lefèvre et al. (2018) Lefèvre, M., Lebonnois, S., & Spiga, A. 2018, Journal of Geophysical Research: Planets, 123, 2773
  * Lefèvre et al. (2020) Lefèvre, M., Spiga, A., & Lebonnois, S. 2020, Icarus, 335, 113376
  * Lefèvre et al. (2021) Lefèvre, M., Turbet, M., & Pierrehumbert, R. 2021, The Astrophysical Journal, 913, 101
  * Merlis et al. (2016) Merlis, T. M., Zhou, W., Held, I. M., & Zhao, M. 2016, Geophysical Research Letters, 43, 2859
  * Mousis et al. (2020) Mousis, O., Deleuil, M., Aguichine, A., Marcq, E., Naar, J., Aguirre, L. A., Brugger, B., & Gonçalves, T. 2020, The Astrophysical journal letters, 896, L22
  * Otegi et al. (2020) Otegi, J., Bouchy, F., & Helled, R. 2020, Astronomy & Astrophysics, 634, A43
  * Pierrehumbert (2010) Pierrehumbert, R. T. 2010, Principles of planetary climate (Cambridge University Press)
  * Pierrehumbert & Ding (2016) Pierrehumbert, R. T. & Ding, F. 2016, Proceedings of the Royal Society of London Series A, 472, 20160107
  * Sergeev et al. (2020) Sergeev, D. E., Lambert, F. H., Mayne, N. J., Boutle, I. A., Manners, J., & Kohary, K. 2020, The Astrophysical Journal, 894, 84
  * Song et al. (2021) Song, Q., Yang, J., Luo, H., Li, C., & Fu, S. 2021, arXiv e-prints, arXiv:2108.04143
  * Tsiaras et al. (2019) Tsiaras, A., Waldmann, I. P., Tinetti, G., Tennyson, J., & Yurchenko, S. N. 2019, Nature Astronomy, 3, 1086
  * Turbet et al. (2021) Turbet, M., Bolmont, E., Chaverot, G., Ehrenreich, D., Leconte, J., & Marcq, E. 2021, Nature, 598, 276
  * Wing et al. (2017) Wing, A. A., Emanuel, K., Holloway, C. E., & Muller, C. 2017, Shallow clouds, water vapor, circulation, and climate sensitivity, 1
  * Yamashita et al. (2016) Yamashita, T., Odaka, M., Sugiyama, K.-i., Nakajima, K., Ishiwatari, M., Nishizawa, S., Takahashi, Y. O., & Hayashi, Y.-Y. 2016, Journal of the atmospheric sciences, 73, 4151
  * Zahnle et al. (1988) Zahnle, K. J., Kasting, J. F., & Pollack, J. B. 1988, Icarus, 74, 62
  * Zeng et al. (2019) Zeng, L., Jacobsen, S. B., Sasselov, D. D., Petaev, M. I., Vanderburg, A., Lopez-Morales, M., Perez-Mercader, J., Mattsson, T. R., Li, G., Heising, M. Z., Bonomo, A. S., Damasso, M., Berger, T. A., Cao, H., Levi, A., & Wordsworth, R. D. 2019, Proceedings of the National Academy of Science, 116, 9723
  * Zhang et al. (2017) Zhang, X., Tian, F., Wang, Y., Dudhia, J., & Chen, M. 2017, The Astrophysical Journal Letters, 837, L27