© The Authors 2024

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 via 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 below 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, which can be hindered by other dynamical processes, needs to be able to develop. For example, when there is a compositional gradient in the atmosphere – whether it be 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 2024). 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,

(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_υ 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^⋆ 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 is higher than the critical concentration, q_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 temperature is the same in 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 mostly been studied analytically with unidimensional arguments, and many caveats remain:

• The analytical theory of Leconte et al. (2017) assumes a fully saturated medium, whereas we know that on Earth, moist convective regions are all but fully saturated. We do not know whether a stable layer would really form in a realistic, sub-saturated atmosphere.

• 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 atmospheres of ice giants.

• We do not know how convection inhibition affects the atmosphere of Neptune-like exoplanets or if this is observable with our current methods and observatories.

To address these issues, we developed 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 atmospheres of Solar System ice giants have radiative timescales on the order of decades to centuries, we focused on simulating the atmosphere of a prototypical temperate sub-Neptune: K2-18b (Cloutier et al. 2017). This both enables the equilibration of the simulation with reasonable computing resources and helps us shed light on the atmospheric structure of such temperate- and 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 show that our more realistic 3D simulations largely support the thermal structure envisioned in the earlier study by 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 (Sec. 4). We then discuss how the chemical composition of temperate Neptune-like planets is affected by convection inhibition, potentially acting as a tracer of this process (Sect. 5).

Interestingly, during the writing of this manuscript, new near- and mid-infrared transit observations of K2-18 b with the James Webb Space telescope (JWST) were released (Madhusudhan et al. 2023). One of the most intriguing features of these observations is the potential non-detection of NH₃, which has been interpreted as a 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 tried to quantify in more detail 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 Sec. 6 we use our updated 1D model to place 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, which are relatively hard to achieve around a late-type star. We further demonstrate that the kinds of aerosol properties needed to create such a high albedo are incompatible with the current transit observations, which 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 fourth version of the Weather Research and Forecast model (hereafter WRF V4) described in Skamarock et al. (2019), which we coupled with the physical parametrizations of the LMD Generic Planetary Climate Model (Generic PCM, i.e., 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 simulated the compressible, non-hydrostatic Euler equations using the philosophy described in Laprise (1992), which uses a terrain-following hydrostatic-pressure vertical coordinate defined as

(2)

where π^t and 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 .

The hydrostatic pressure, π, 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, φ = gz, and the specific volume of the fluid, α, one can write the equivalent of the hydrostatic equation:

(3)

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

(4)

such that

(5)

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

(6)

and the geopotential φ that is not written in flux form as it is not a conserved quantity. The usual velocities are denoted v = (u, v, w), and we define the horizontal velocity for sake of compacity as is the potential temperature, p₀ is an arbitrary reference pressure, R the specific constant of the gas, and r_i is the mass mixing ratio of the various tracer species with respect to dry air (in particular water vapor with index υ and condensed water with index c). For the sake of compactness of notation, the specific concentration of the tracers, q_i = r_i/(1 + r_υ), 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:

(7)

(8)

(9)

(10)

(11)

(12)

(13)

Here the various differential operators are defined as follows:

(14)

(15)

(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_i the specific source/sink rates for tracer i (in our case, the vapor and condensed phases).

Our set of equations differ from those of Daley-Yates et al. (2021) or Habib & Pierrehumbert (2024), 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:

(17)

where R_υ and R_d are the specific gas constant of the vapor and dry air, and γ = 1/(1 – R/c_p) is kept constant. We notice that the ratio R/cp 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_dR_d + q_υR_υ will change with the vapor concentration, the c_p used in the right-hand side of Eq. (12) must be varied accordingly, following c_p = q_dc_p,d + q_υc_p,υ. 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 α_d/α = 1 + ∑_i r_i terms that account for the mass loading effect of both phases.

2.2 Dynamical core and numerical implementation

To solve these equations, we used WRF V4, as 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 64 × 64 × 256 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 Courant–Friedrichs–Lewy (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 turned off all the parametrizations of sub-grid-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 fifth-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 parameterizations

The various physical source terms in the equations of motion – in particular the diabatic heating and the tracer source terms, Q and S_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 WRF V4 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 modeled 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_i) and the physical parametrizations in PCM that use tracer specific concentration with respect to the whole gas (q_i = r_i/(1 + r_υ)) – 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 (π) 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, including dry and moist convective adjustment, and sub-grid-scale humidity distribution. 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 10⁷ 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 collision-induced absorptions (CIA) for H₂-H₂, H₂-He are taken from the HITRAN database (Richard et al. 2012). The H₂O-H₂O and H₂O-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 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, which we kept arbitrarily small in this first study to avoid the radiative feedback of clouds. Precipitations fall instantaneously. Because the planets we modeled 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 micro-physics, we wanted to keep a physically motivated model without any free parameter. We thus assumed 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 this is clearly an upper limit on the pressure of the re-evaporation 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 assumed that falling precipitations reach their equilibrium velocity instantaneously such 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_int, at each physical timestep. We kept the bottom of the model below the re-evaporation level and verified 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, here we try to identify what the possible sources of such energy losses are and to identify possible areas of improvement for the model in the future.

We 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

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

• 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.

• 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 the variation of c_p into account in Eq. (12), this variation entails a slight non-conservation, especially when the mixing ratio of vapor is not negligible.

3 Atmospheric dynamics and thermal structure of temperate Neptunes

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 Sec. 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× 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, and hence on the composition through the mean molecular weight. We thus decided to use a simple approach to convert metallicity into molecular abundances: we assumed that the quenching in our temperate planet will occur near 10 bar and 1000 K and used the chemical code from Venot et al. (2020) to compute the abundances. With this approach, the volumic concentration of the main absorbers are and . For water, this corresponds to a deep specific concentration q_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 × 10⁻³ kg mol⁻¹.

Because K2-18b receives an insolation that is very close to the runaway greenhouse threshold (see the 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 used an average insolation of 175 Wm⁻² to run our 3d simulations. This would correspond to an effective bond albedo of ≈0.5, which allows us to be in the right regime to test the conclusions of Madhusudhan et al. (2023) in Sec. 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_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, penetrates this layer, so 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 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, we 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⁻⁵ kg s⁻¹. As could be expected, the re-evaporation of these rainfalls occurs below the base of the convective plumes. Indeed, a first requirement for re-evaporation is to reach an unsaturated region below the updraft.

An interesting finding is that the level of re-evaporation seems to set the bottom of the stable layer. Indeed, below this level, condensation is impossible, so there are no vapor sources or sinks. As a result, any vapor gradient in this lower region – called the dry troposphere in Fig. 1 – would be short-lived and stellar energy deposited at depth is sufficient to trigger standard, dry convection and efficiently mix the atmosphere. But this raises the question as to why the dry troposphere cannot extend above the re-evaporation level. We should 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 re-evaporation 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 Leconte et al. (2017) model: they did not consider that the precipitations formed in the moist convective layer would necessarily re-evaporate below the stable layer, meaning 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 and thus 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.

Fig. 1 Atmospheric structure of the baseline simulation. From left to right: temperature, water vapor specific concentration, vertical velocity (in m s−1), and relative humidity. The first two panels show horizontal and temporal averages. The dash-dotted black line in the first panel shows the standard moist-adiabat profile for the same conditions. In the second panel, the dotted and dashed black 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 facilitate the comparison of the altitudes of the various features. The rising moist plume (with maximal velocities of around 8 m s−1) is mirrored by a descending cold plume in the dry region (-15 m s−1) caused by the re-evaporation of rains at the bottom of the stable layer.

Fig. 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 the calculation method.

3.4 Turbulent mixing

Even at the (500 m) scale resolved in our simulations, turbulence spontaneously appears¹. 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 estimated the equivalent mixing coefficient (the so-called eddy diffusivity or K_zz) in the simulations using a passive tracer whose concentration (r_tra) is fixed at the surface and undergoes advection by the flow and local exponential decay with a timescale τ_tra. Then we used two different methods. In the first one, the eddy flux approach, we just computed the average turbulent flux of tracer in the simulations and assumed that

(18)

where ⟨⟩ denotes temporal and horizontal averaging. In the second method, or integral the approach, we used the fact that if the mixing exhibits a diffusive-like behavior, the steady-state tracer profile should obey the following law:

(19)

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

(20)

We note in passing that because the transport is not perfectly diffusive (especially in convective regions), this approach can yield negative K_zz, and so 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.

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⁻², which is comparable to values inferred for the lower stratosphere on Earth (Allen et al. 1981). As discussed in more detail in Sec. 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 re-evaporation in that region, so, in a steady state, the upward vapor flux (∝ K_zz(dq/dz)) 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.

Fig. 3 Average vertical profiles of simulated eddy diffusion coefficients and velocities. 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 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.

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⁻¹ m s⁻¹). The only exception occurs in the thin boundary layer just above the stable layer below, around 2 × 10⁴ 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 that 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 re-evaporation 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 – a factor of 2–3 higher in velocity – because both the vapor loading and the cooling due to re-evaporation (which can reach several kelvins) accelerate the plume downward.

Finally, we 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 Stabilization and sub-saturation

An intriguing feature of our simulations is that the atmosphere is sub-saturated 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 is a perfect example of this (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 sub-saturated 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 sub-saturated 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)² 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 focus mainly on the new features of the model that pertain to the inhibition of convection.

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 (θ). 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,

(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 (and hence a 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 (2024) 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_υ > q_cri. In the convective regions, the excess vapor condensed during the adjustment process is assumed to instantaneously rain down to the re-evaporation 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 used two tracers that are mixed by convection as described above. For the thermodynamics and the microphysics of clouds and precipitations, we used exactly the same parametrization as in our 3D model, which is described in Sec. 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(d/dz), 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 m²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 implemented an alternative formulation, referred to as the baseline scenario, where the eddy diffusivity assumes a larger value () just at the top boundary of the dry convective region (also called the radiative-convective boundary, whose pressure is p_<sc>rcb</sc>) and drops off very rapidly as some high power α of the pressure before settling to a lower constant value higher up:

(22)

The value of α is set to reproduce the sharp decrease in the turbulence above the convective region seen in Fig. 3. This yields α ≈ 13. Then, and are tuned to reproduce the vapor flux and the thermal gradient just above the radiative-convective boundary, respectively. This yields ≈ 0.08 m² s⁻¹ and ≈ 3 m² s⁻¹, which is also consistent with the numerical values from the 3D simulation.

Fig. 4 Atmospheric structure of the 1D model for the baseline case. From left to right: temperature, water vapor specific concentration, dynamical fluxes, and radiative flux 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 the convective zone, where the value of the diffusivity is irrelevant because mixing is not diffusive. 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 panel shows that most of the radiative cooling to space occurs in the stable layer and above.

4.3 1D-3D comparison

In Fig. 4, we present the atmospheric structure for the baseline scenario. We 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 in 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 × 10⁴ 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 sub-saturated on average.

Figure 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. This also causes an increase in 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 ~200 K.

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 1000K at the 100 bar level.

To quantify the effect on the chemistry, we computed the chemical composition of the two model atmospheres shown in Fig. 6 using the chemistry module of Exo-REM (Blain et al. 2021) assuming K_zz = 10⁴ m² s⁻¹. We took 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 ≈20 bar for the baseline scenario compared to ≈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 CO₂ in the stratosphere, most of the carbon forming CH₄. On the contrary the case where inhibition is taken into account, CO, CO₂, and CH₄ are all in detectable quantities. This case also predicts a much lower abundance of NH₃. 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.

Fig. 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 are much smaller than the effect of convection inhibition itself, which can be seen by comparing with the no-inhibition case (dash-dotted curve).

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 developed 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

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

• 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,

• 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 H₂-dominated atmosphere. This question has been investigated by Innes et al. (2023), who showed that the greenhouse effect of an H₂-dominated atmosphere is much greater than that of an N₂ 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 reframed the issue by asking what the planetary Bond albedo would need to be for the absorbed insolation to remain below the runaway greenhouse threshold. To answer this question, we performed the following experiment : we equilibrated our 1D model for a given surface pressure with a stellar irradiation equal to F_inc(l – A), where F_inc = 342 W m⁻² is the average insolation received by K2-18 b and A is an albedo that we chose arbitrarily. In this specific experiment, we did 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; therefore, 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⁻² ppmv, which seem to be the best match for the JWST observations of Madhusudhan et al. (2023). Second, we changed 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 treated 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 kept the total surface pressure constant during the evolution. We then let the model evolve either until it reached thermal equilibrium, in which case we deemed the ocean stable, or the mixing ratio of vapor reached unity at the surface, which we took 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 m⁻², whereas our last stable insolation is ≈130 W m⁻² (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 verified this by rerunning 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; removing it would cause the stable layer to become entirely 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 ≈410 K – which corresponds to a planet averaged absorbed/thermally emitted flux of ≈1300 W m⁻². This is to be compared to our limit for the 1 bar case, which is estimated to be ≈230 K (150 W m⁻²). To put this into context, the current equilibrium temperature of the Earth is ≈255 K (240 W m⁻²), and recent estimates of the runaway greenhouse limit for Earth-like planets yield estimates between ≈260 and 270 K (270–300 W m⁻²), 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 H₂ compared to N₂, 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 H₂ 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, counterintuitively, 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.

Fig. 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×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.

Fig. 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 in which 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⁻² 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) nor 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 computed eclipse and transmission spectra of the fiducial model of K2-18 b discussed above, adding 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 H₂. We note that n_haze = 1 corresponds to the fiducial clear atmosphere model.

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⁵. However, we see in the bottom panel of Fig. 8 that the increased scattering also mechanically suppresses the molecular methane features in the near-infrared. We quantified the amplitude of these features by taking the difference between the transit depth at two absorption peaks (1. and 1.16 μm) and the depth in the nearest windows (1.08 and 1.28 μm, respectively). These 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 gobservation 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 such as Titan or the giant planets. We thus deem this possibility rather unlikely as well.

Fig. 8 Eclipse (top) and transit (bottom) spectra of our models of K2-18 b with hazes parametrized through the nhaze factor (see the main 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.

Fig. 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 nhaze. 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 ≳ 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 atmospheres. Being integrated into the ecosystem of the Generic PCM model, it is very flexible and can be easily adapted to a wide diversity of planets. We then investigated how moist convection behaves in H₂-dominated atmospheres using K2-18 b as a prototypical temperate Neptune-like planet.

Our main general findings are:

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

• The atmospheric structure envisioned by Leconte et al. (2017) – that is, 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.

• 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.

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

• The higher temperature at depth also decreases the limiting insolation at which a surface liquid ocean can be stable under an H₂-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, for example the methane in Uranus and Neptune (Clément et al. 2023) and the iron or silicates near the core of Neptune-like planets (Markham et al. 2022; Misener & Schlichting 2022).

The fact that the higher temperatures at depth further reduce the insolation at the so-called inner edge of the habitable zone for H₂-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 for there to be a configuration that explains 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 mechanisms to explain the lack of this molecule other than a shallow atmosphere above a liquid ocean.

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. M.L. acknowledges funding from the European Union’s Horizon Europe research and innovation program under the Marie Sklodowska-Curie grant agreement “MuSICA-V”.

Appendix A Calculation of the energy budget

To help visualize the flow of energy in our simulations, we computed 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_i). The net upward flux for process i at any level defined by the hydrostatic pressure, π, is defined as the integral of the heating between that level and the surface:

(A.1)

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_dyn ≡ c_p⟨ρw′θ′⟩, where the primes denote departures from the average values.

References

All Tables

All Figures

Fig. 1 Atmospheric structure of the baseline simulation. From left to right: temperature, water vapor specific concentration, vertical velocity (in m s−1), and relative humidity. The first two panels show horizontal and temporal averages. The dash-dotted black line in the first panel shows the standard moist-adiabat profile for the same conditions. In the second panel, the dotted and dashed black 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 facilitate the comparison of the altitudes of the various features. The rising moist plume (with maximal velocities of around 8 m s−1) is mirrored by a descending cold plume in the dry region (-15 m s−1) caused by the re-evaporation of rains at the bottom of the stable layer.

In the text

	Fig. 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 the calculation method.
In the text

Fig. 3 Average vertical profiles of simulated eddy diffusion coefficients and velocities. 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 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.

In the text

Fig. 4 Atmospheric structure of the 1D model for the baseline case. From left to right: temperature, water vapor specific concentration, dynamical fluxes, and radiative flux 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 the convective zone, where the value of the diffusivity is irrelevant because mixing is not diffusive. 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 panel shows that most of the radiative cooling to space occurs in the stable layer and above.

In the text

Fig. 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 are much smaller than the effect of convection inhibition itself, which can be seen by comparing with the no-inhibition case (dash-dotted curve).

In the text

Fig. 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×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.

In the text

Fig. 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 in which 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.

In the text

	Fig. 8 Eclipse (top) and transit (bottom) spectra of our models of K2-18 b with hazes parametrized through the nhaze factor (see the main 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.
In the text

Fig. 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 nhaze. 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 ≳ 0.5). The relatively high band amplitude found in the observations rules out haze that are reflective enough to stabilize a liquid surface ocean.

In the text