Convection modeling of pure-steam atmospheres

Condensable species 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. 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 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 model parameters, 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.


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 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 . Extremely closein 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 xianyu.tan@physics.ox.ac.uk the major condensable constituents being CO 2 and CH 4 , respectively.
Several studies have targeted the climate dynamics with non-dilute condensable vapor.  examined the energy budget associated with non-dilute condensation and precipitation, and proposed a scheme for convection parameterization applicable in general conditions.  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 nondilute conditions. Yamashita et al. (2016) performed two-dimensional convection modeling for a pure-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 . 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.
We use a plane-parallel, twostream 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 m 2 kg −1 , same as that used in . 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.

Condensation, rainout and evaporation
In  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, 1 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 ) 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: where ρ 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 1 ρ in Eq.
(1) drops out, and there is a conserved quantity ρk − p before and after condensation. In a pure-steam atmosphere, this quantity (denoted by A) is A ≡ ρ(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 cool-ing. Then where ρ 1 = p1 RT1 . 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 in which we made use of conservation of mass (ρ 1 = ρ 2 ). With 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 ρ 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 S to overcome surface tension (Houze 2014), in which 1 + S = p 1 /p sat . We allow 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 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 A. Pre-evaporation temperature and pressure in the evaporation regions are denoted as T 3 and p 3 , and the total density is ρ 3 = p3 RT3 + ρ cond , where ρ cond is the density of the soon-to-be-evaporated precipitation.

The conserved quantity is
After evaporation, temperature T 4 and pressure p 4 are finally obtained via 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.

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 δz ranges between 600 and 680 m depending on the surface temperature. Our canonical models have 320 × 320 × 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 m s −2 , the heat capacity at constant pressure for water vapor c p = 1870 J kg −1 K −1 and for condensates c pl = 4190 J kg −1 K −1 , the specific gas constant R = 461.9 J kg −1 K −1 , and latent heat L = 2.5 × 10 6 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  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.

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 S = 10 −7 and three surface temperatures of 600, 700 and 800 K. In both condensing and dry regions, the horizontal rootmean-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 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 ( 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 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 ( 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 .
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 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 puresteam 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). 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 downwardmoving 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 , 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 down- Radiative-cooling driven condensation Overshoot driven condensation Evaporation Fig. 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 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 ward 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.

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 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 S = 0.01, 0.1 and 0.2. In the pure-steam condition, the characteristic permitted temperature variation before and after condensation is δT ≈ ST 2 R/L when S 1. Higher 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 ( 1 mbar), higher S leads to larger RMS vertical velocity (left panel) and larger RMS temperature variation ( 0.2 K for S ≥ 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 S. In the lower dry region, the RMS vertical velocity slightly decreases with increasing S, although the decrement appears to saturate when S 0.1. This results in weaker overshoot-driven condensation when S is larger as shown on the right. At a given time, condensation in the radiative cooling region is widespread in the case with S = 10 −7 but is highly sparse when 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 S. Given the same radiative cooling rates, area fraction of condensation decreases with increasing S to balance radiative cooling. Interestingly, the case with S = 10 −7 shows a larger domain-mean condensation rate at the photosphere that is almost balanced by radiative cooling, while others with S 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 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.
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 radiativeconvective equilibrium modeling applied in Earth's tropical conditions suggest that rotation triggers hurricanelike 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 . 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 × 10 −4 s −1 (the value at Earth's North pole) and another one 10 times of that f = 1.45 × 10 −3 s −1 . The domain-mean properties are shown in the bottom row of Figure 3. All three models have a supersaturation ratio 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 vor-tices. Figure 4 shows snapshots of vertical velocity at different pressure levels from the model with f = 1.45 × 10 −3 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).

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  and its application in  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.  Fig. 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 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 × 10 −4 and 1.45 × 10 −3 s −1 , and a supersaturation ratio S = 0.1. 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 highresolution, 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 ( 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 im-plies weak transport of trace species from the deep atmosphere, which may have interesting implications for the chemistry of such atmospheres.