© The Authors 2023

1 Introduction

Understanding the set of conditions that allow rocky planets to have liquid water on their surface in the form of lakes, seas, or oceans is a major scientific step in determining the fraction of planets potentially suitable for the emergence and development of life as we know it on Earth. Based on our experience on Earth, liquid water as a solvent is thought to be a necessary (but not sufficient) condition for the existence of life, along with the presence of an energy source and nutrients C, H, N, O, P, and S (Forget et al. 2013; Shields et al. 2016). While liquid water is now known to exist in the subsurface of a large sample of planetary bodies (icy moons, Pluto, Triton; see Lunine 2017 and references therein), planets that have liquid water on their surface have two additional favorable properties (Turbet & Selsis 2023). First, liquid water is exposed to stellar radiation and second, liquid water is in interaction with the atmosphere, which facilitates the remote search for tracers of life (using the exoplanetary atmosphere to search for biomarkers).

This has motivated numerous theoretical studies to determine the range of parameters required for a planet to sustain surface liquid water reservoirs (e.g., type of star, incoming stellar radiation) and is the main driver for the definition of the habitable zone (HZ; Kasting et al. 1993; Kopparapu et al. 2013, 2014). The HZ is now a widely used first-order criterion to determine a zone to search for life-compatible planets, and to evaluate the potential habitability of detected planets (see, e.g., Borucki et al. 2011).

To date, most numerical climate studies in the field have focused on evaluating the conditions needed for a planet to stabilize surface liquid water, whether using one-dimensional (1D) climate models (Kasting et al. 1993; Kopparapu et al. 2013, 2014; Ramirez & Kaltenegger 2014, 2016, 2017, 2018; Yang et al. 2016; Koll & Cronin 2019; Chaverot et al. 2022) or three-dimensional (3D) global climate models (GCMs; Wordsworth et al. 2011; Yang et al. 2013, 2014, 2019a,b; Shields et al. 2013; Bolmont et al. 2016; Kopparapu et al. 2016, 2017; Turbet et al. 2016, 2017; Boutle et al. 2017, 2020; Wolf 2017; Wolf et al. 2017, 2020; Turbet 2018; Way et al. 2018; Del Genio et al. 2019; Kodama et al. 2019; Colose et al. 2021). However, this approach implicitly assumes surface liquid water was available in the first place. This assumption was recently challenged in Turbet et al. (2021), where the conditions of surface ocean formation by condensation of the primordial water vapor reservoir were explored. Terrestrial planets are indeed believed to form hot due to their initial accretion energy, and thus to cross a magma ocean stage (Zahnle et al. 1988; Abe 1997; Hamano et al. 2013; Lebrun et al. 2013; Salvador et al. 2017; Lichtenberg et al. 2023) where superficial water is present only in the form of vapor before evolving toward their final state. For liquid water oceans to ever appear on a planetary surface, water initially present in the young and warm planetary atmosphere must be able to condense on the surface in the first place (Hamano et al. 2013; Lebrun et al. 2013; Turbet et al. 2021).

The conditions leading to the condensation of a water ocean after the magma ocean phase were first studied with 1D numerical climate models (Abe & Matsui 1988; Hamano et al. 2013; Lebrun et al. 2013; Salvador et al. 2017), which neglected the effects of atmospheric dynamics and clouds. Recently, Turbet et al. (2021) simulated water condensation for the first time in a GCM, specifically for the cases of Earth and Venus, but with potentially far-reaching consequences for a broader range of planets (Kasting & Harman 2021). They showed that in water-rich atmospheres water clouds tend to preferentially form on the nightside and at the poles, owing to the strong subsolar water vapor absorption. Clouds thus have a strong net warming effect that inhibits surface water condensation even at modest insolations, with strong implications for the early habitability of Earth and Venus.

Here we seek to extend the work of Turbet et al. (2021) by extending GCM simulations to the cases of exoplanets, and especially to those orbiting stars cooler than the Sun. Terrestrial planets orbiting low-mass, cool, red stars (also known as M stars) are indeed the best targets we currently have for remote sensing of their atmospheres, and thus for searching for tracers of habitability and life (Morley et al. 2017; Kaltenegger 2017; Lustig-Yaeger et al. 2019; Wunderlich et al. 2019; Fauchez et al. 2019). These planets are thus potentially unique natural laboratories where we can confront our theories on planetary habitability with observations, in the event that they were able to preserve their atmospheres (Zahnle & Catling 2017).

2 Method

All the numerical experiments presented in this study are based on the Generic Planetary Climate Model (PCM), which is historically known as the Laboratoire de Météorologie Dynamique (LMD) Generic GCM.

2.1 The Generic PCM

The Generic PCM is a sophisticated 3D GCM that was previously developed and used to simulate a wide range of exoplanetary atmospheres ranging from that of temperate rocky planets to warm mini-Neptunes (Wordsworth et al. 2011; Leconte et al. 2013b; Charnay et al. 2015a,b, 2021; Turbet et al. 2016; Turbet 2018; Fauchez et al. 2019). Specifically, the model was adapted to simulate water-rich planetary atmospheres (Leconte et al. 2013a) with numerous applications: studying the runaway greenhouse effect on Earth (Leconte et al. 2013a), post-impact atmospheres on Mars (Turbet et al. 2020c), water cloud formation on temperate mini-Neptunes (Charnay et al. 2021), and water condensation on early Earth and Venus (Turbet et al. 2021).

We summarize the main features of the model below. The GCM includes a complete radiative transfer that takes into account the absorption and scattering by the atmosphere, the clouds, and the surface. The radiative transfer calculations, which are based on the correlated-k approach (Fu & Liou 1992), use 55 spectral bands in the thermal infrared (from 0.65 to 100 μm, designed to capture the thermal emission of the planet) and 45 bands in the visible domain (from 0.3 to 6.5 μm, designed to capture the incoming stellar radiation). Opacity tables were computed as in Turbet et al. (2021) using the HITRAN 2008 database for wavelengths above 1 μm and the HITEMP 2010 database for those below. We made this choice because most of the differences between HITRAN and HITEMP are below 1 μm (see comparisons done in Goldblatt et al. 2013), and that this significantly reduces the computation time of the opacity tables. The opacity tables also account for the continuum absorptions (from collision-induced absorptions, dimer, and far-wings) of N₂–N₂, H₂O–H₂O, and H₂O–N₂, using the HITRAN collision-induced absorption (CIA) database (Karman et al. 2019) and the MT_CKD v3.5 database (Mlawer et al. 2012). Subgrid-scale dynamical processes (turbulent mixing and convection) were parameterized as described in Leconte et al. (2013a). Moist convection was taken into account following a moist convective adjustment scheme originally derived from Manabe & Wetherald (1967) and later generalized for water-dominated atmospheres in Leconte et al. (2013a). Relative humidity evolves freely and is limited to 100% (no supersaturation). In practice, when an atmospheric grid cell reaches 100% saturation and the corresponding atmospheric column has an unstable vertical temperature profile, the moist convective adjustment scheme is performed to obtain a stable moist adiabatic lapse rate. When condensing, water vapor forms liquid water droplets and/or water ice particles, depending on the atmospheric temperature and pressure, and forms clouds. We used a fixed number of activated cloud condensation nuclei (CCNs) per unit mass of air to determine, based on the amount of condensed water, the local, effective radii of H₂O cloud particles (Turbet et al. 2021). The effective radius is then used to compute the radiative properties and the sedimentation velocity of cloud particles. Water precipitation is divided into rainfall and snowfall, depending on the nature of the cloud particles, determined solely by the atmospheric temperature. Rainfall is parameterized to account for the conversion of cloud liquid droplets to raindrops by coalescence with other droplets (Boucher et al. 1995). Rainfall is considered to be instantaneous (i.e., it goes directly to the surface), but it can evaporate while falling through subsaturated layers (Gregory 1995). In hot and steamy simulations, re-evaporation of precipitation is always complete, which means that precipitation always fully evaporates in the dry lower atmosphere before it reaches the surface. The snowfall rate is calculated using the sedimentation velocity of particles, assumed to be equal to the terminal velocity that we approximate by a Stokes law modified with a “slip-flow” correction factor (Rossow 1978).

The simulations presented in this paper were performed at a spatial resolution of 64 × 48 in longitude × latitude. In the vertical direction, the model is composed of 40 distinct atmospheric layers for a 10 bar surface pressure atmosphere, ranging from the surface up to a few pascals. The exact number of atmospheric layers and model top pressure vary slightly depending on the surface pressure considered. The dynamical time step of the simulations is ~15 s, but can vary slightly from one simulation to another. The radiative transfer and the physical parameterizations (such as condensation, convection) are calculated every 100 and 5 dynamical time steps, respectively.

The main input parameters of our simulations are summarized in Table 1. Additional details on the model and the parameterizations used in this study can be found in Turbet et al. (2021), and all the references therein. Compared to Turbet et al. (2021), we made several changes: We varied the type of host star; we varied the rotation period and rotation mode into synchronous rotation, and fixed the obliquity to 0°; we varied the values of the masses and radii of planets, depending on the numerical experiments, either to match real planets (in particular the TRAPPIST-1 inner planets) or to make comparisons with previous studies (in particular Yang et al. 2013); we varied the internal heat flux by simply adding a constant (through time, and at all latitudes and longitudes) heat flux at the bottom of the model. These changes are detailed in Table A.1.

2.2 Strategy of the numerical experiments

Using the Generic PCM, we set up a series of 3D GCM simulations, assuming a N₂+H₂O-rich atmospheric composition, designed to identify the conditions required to trigger water condensation from the atmosphere onto the surface.

For this we followed a strategy similar to that used in Turbet et al. (2021). Given a host star (T_eff between 2600 and 5780 K), we first start with a GCM simulation with partial pressures $P_{{\text{N}}_2}=1$ bar and $P_{{\text{H}}_2\text{O}}=10$ bar (in our baseline scenario), with a high incoming stellar radiation (500 W m⁻², corresponding to 1.47 of the insolation on Earth). A value of $P_{{\text{H}}_2\text{O}}=10$ bar is far below the possible content of a primordial steam atmosphere (Earth oceans would correspond to ~300 bar). Using 10 bar ensures the feasibility of the GCM simulations, and should be sufficient to reach conclusions valid for higher contents, as demonstrated in Sect. 3.3. For this baseline simulation we assume an isothermal cloud-free atmospheric temperature profile at 10³ K for the initial state. While the upper atmospheric layers evolve very rapidly toward equilibrium, we used a strategy proposed in Turbet et al. (2021) to converge the deepest layers of the atmosphere, which is detailed below. We ran the GCM for typically 10 Earth years starting from the initial state described earlier, and then proceeded as follows. In step 1 we ran the GCM for 2 Earth years. In step 2 we extrapolated the evolution of the temperature field over n_{Earth years} (the number of Earth years on which we decided to extrapolate the temperature field) using $$T_{i,j,k}\left(t+n_{\text{Earth\hspace{0.17em}years}}\right)=T_{i,j,k}\left(t\right)+n_{\text{Earth\hspace{0.17em}years}}\times \text{Δ}{T}_{mean,k},$$(1)

with T_i,j,k the temperature at the i, j,k spatial coordinates (corresponding respectively to longitude, latitude, and altitude coordinates) and ∆T_mean,k the change in the mean horizontal temperature field (averaged over all longitudes and latitudes) at the altitude layer k calculated over the second Earth year in step 1. The value of n_{Earth years} was first set arbitrarily to 100, then 50, 20, 10, and finally 0, when the planetary atmosphere is close to convergence. The transition from one n_{Earth years} value to another was chosen manually to minimize the convergence time, while avoiding numerical instabilities. In step 3 we repeated steps 1 and 2 until the top of the atmosphere (TOA) radiative imbalance was lower than ~1 W m⁻². We also checked that the surface temperatures were stable.

We repeated the same procedure for simulations at lower incoming stellar radiation (ISR) until runaway water condensation was reached, except that we used the final equilibrium state of simulations at (previous) higher ISR values as the initial state. In comparison to Turbet et al. (2021), the convergence time of the simulations is much longer, in particular for simulations around the coolest stars (TRAPPIST-1 and Proxima Cen). This is due to the fact that for cooler stars, a very small fraction of the stellar flux reaches the low atmosphere and the surface, due to the near-IR absorption by water vapor. For instance, for a GCM simulation of 10 bar of water vapor around TRAPPIST-1 near the water condensation limit, it takes about 150 Earth years of GCM simulations to converge the simulation. This is equivalent to around 2000 Earth years of evolution, taking extrapolation into account. We discuss this aspect in more detail in Sect. 3 as it also has an impact on the thermal profiles (Selsis et al. 2023). We note that Turbet et al. (2021) performed some sensitivity tests to ensure that the accelerated convergence scheme was working effectively. We thus did not perform any additional tests here.

On the whole, and taking advantage of the convergence scheme presented above, we managed to converge all the simulations presented in this work using a total of 1300 kh CPU on the French supercomputer OCCIGEN. This corresponds to a total of 70 simulations, which are listed in this work, in addition to numerous preliminary and sensitivity simulations (included in the total CPU used and in the carbon budget provided in the Acknowledgments).

Finally, we note that for each combination of ISR and star type, we calculated the rotation period T self-consistently with the assumption of synchronous rotation, using the formula $$T=2\pi \frac{{\left(d_{\odot -\oplus }^2\left(\frac{L_{\star }}{L_{\odot }}\right)\left(\frac{F_{\oplus }}{F_p}\right)\right)}^{\frac{3}{4}}}{\sqrt{G{M}_{\star }}},$$(2)

where $d_{\odot -\oplus }^2$ is the Sun-Earth distance (in m); L_★ and L_⊙ the stellar and solar luminosity, respectively; F_p and F_⊕ the stellar flux received at the TOA of the planet and on Earth, respectively; G the gravitational constant (in m³ kg⁻¹ s⁻²); and M_★ the stellar mass (in kg). The formula is simply a reformulation of Kepler’s third law of planetary motion and assumes that the rotation period is equal to the orbital period (i.e., it assumes synchronous rotation). The luminosity of the star is taken as equal to that at the zero age main sequence (ZAMS), which is at the beginning of the main sequence phase (which can be much lower than during the pre-main sequence (PMS), in particular for M stars. This is done because here we are interested in determining the minimum ISR at which a planet has condensed its water reservoir at the start of the main sequence. This roughly corresponds to the minimum luminosity during the lifetime of stars. For each star type, we used the following mass-luminosity properties:

• For K5 stars: L_K5 = 0.17 L_⊙, M_K5 = 0.7 M_⊙;

• For M3 stars: L_M3 = 1.2 × 10⁻² L_⊙, M_M3 = 0.31 M_⊙ (the stellar properties of L98-59; Kostov et al. 2019);

• For M5.5 stars: L_M5.5 = 1.5 × 10⁻³ L_⊙, M_M5.5 = 0.12 M_⊙ (the stellar properties of Proxima Centauri; Faria et al. 2022);

• For M8 stars: L_M8 = 5.5 × 10⁻⁴ L_⊙, M_M8 = 0.09 M_⊙ (the stellar properties of TRAPPIST-1; Agol et al. 2021).

2.3 Ensemble of simulations

We set up a grid of GCM simulations (summarized in Table A.1) designed to explore the conditions of primordial ocean condensation for various types of terrestrial exoplanets. The simulations are listed below:

• SUN (from Turbet et al. 2021), K5, M3, Pcen, and T1: Designed to explore surface water condensation on an Earth-like planet (1M_⊕, 1R_⊕) synchronously rotating (except for the Sun) around its host star. The simulations assume respectively a G2V star (Sun spectrum), a K5 star (4400 K blackbody spectrum), an M3 star (3400 K blackbody spectrum), an M5.5 star (Proxima Cen PHOENIX BT-Settl synthetic spectrum), and an M8 star (TRAPPIST-1 PHOENIX BT-Settl synthetic spectrum). All the simulations assume a H₂O partial pressure of 10 bar, and a N₂ partial pressure of 1 bar.

• M3-YANG: Designed to reproduce the same planetary parameters (2 R_⊕, 60 Earth day rotation period) as used in the standard simulations of Yang et al. (2013), except the simulations assume a warm start (1 bar of N₂, 10 bar of H₂O).

• M3-3bar, M3-1bar, and M3-0.3bar: Designed to explore the impact of water content on the surface water condensation.

• T1b, T1c, T1d, T1e, and Pcen-b: Designed to explore surface water condensation (and observability) on the known exoplanets TRAPPIST-1b, TRAPPIST-1c, TRAPPIST-1d, TRAPPIST-1e, and Proxima b, respectively.

• T1b-3bar and T1b-1bar: Designed to explore the role of water content on the climate (and observability) of TRAPPIST-1b.

• T1b-Fgeo1, T1b-Fgeo5, and T1b-Fgeo25: Designed to explore the role of internal heat flux on the climate (and observability) of TRAPPIST-1b. The internal heat flux is possibly high on TRAPPIST-1b due to tidal heating (Turbet et al. 2018).

2.4 Computation of observables

Transmission and emission spectra were simulated using the Planetary Spectrum Generator (PSG;¹ Villanueva et al. 2018, 2022). The PSG is an online radiative transfer model integrating up-to-date radiative transfer methods and spectroscopic parameterizations, including a realistic treatment of multiple scattering in layer-by-layer spherical geometry. Multiple scattering from aerosols in the atmosphere is computed using the discrete ordinates method, which employs the plane-parallel atmosphere in which the radiation field is approximated by a discrete number of streams angularly distributed. The angular dependence of the aerosol scattering phase function is represented with Legendre polynomials whose number of expansion terms is equal to the number of stream pairs. For each available aerosol type (e.g., Massie & Hervig 2013), look-up-tables of Legendre expansion coefficients are pre-computed using an assumed particle size distribution. The PSG contains billions of spectral lines of thousands of species from several constantly updated spectroscopic repositories, such as HITRAN, ExoMol, JPL, CDMS, and GSFC-Fluor. In this work we used the latest HITRAN database (Gordon et al. 2022). The PSG noise model for the NIRISS, NIRSpec, and MIRI instruments were benchmarked with the JWST Exposure-Time-Calculator (ETC)². The PSG noise model was benchmarked with current post-JWST launch sensitivities by using JWST ETC version 1.7. We used the NASA JWST website³ to establish the optics temperature post commissioning (~40 K). For each instrument, we computed the unsaturated signal-to-noise ratio (S/N) from the ETC by using the appropriate filter, subarray, readout pattern, group per integration, and integrations per exposure. We adjusted the number of exposures, exposure time, readout noise, and dark current in the PSG to match these numbers. The total background level for the observations (e- s⁻¹) was retrieved from the ETC, and we used a consistent background level in the PSG. To estimate the effective instrument throughput efficiencies needed by the PSG, we compared the extracted flux from the source (e- s⁻¹) in the ETC and we adjusted the throughput parameter in the PSG to match these levels. Finally, we checked that the simulated S/N from the ETC matched those we computed with the PSG for a set of representative cases.

To ingest the 3D atmospheric fields (temperature, pressure, volume mixing ratios of molecular species, and mass mixing ratios and particle sizes for liquid and water ice clouds) produced by the Generic PCM and create synthetic spectra, we used the PSG module Global Exoplanetary Spectra⁴ (GlobES; Villanueva et al. 2022). For transmission spectroscopy, GlobES computes the radiative transfer at each GCM latitude chunk across the terminator to create a local spectrum. All the local spectrum are then averaged to create an average terminator spectrum for the planet. For emission spectra, GLobES performs radiative transfer simulations across the whole observable disk, and the individual spectra are integrated considering the projected area of each bin. Because radiative transfer calculations across the full observable disk are very expensive, we horizontally binned GCM data by a factor of 3, as typically done to speed up the calculation, leading to 3 × 3 (latitude × longitude = 9) binning, and a speed up of nine times, while sacrificing a little accuracy. Thermal phase curves were computed at all phase angles (with a step of 5 degrees) along the orbit. Broadband phase curves were integrated over the MIRI 5–25 μm range (see Figs. B.3–B.9). For comparison with MIRI observations in photometric bands, we integrated the high-resolution emission spectra using the spectral transmission function of each MIRI band (for F1280W and F1500W bands, in this work).

3 Results

3.1 Preferential nightside cloud formation and climate multi-stability

Our most striking result is that in all the simulations we performed (listed in Table A.1), water clouds form preferentially on the nightside and at the poles. This is illustrated in Fig. 1 (top row) that shows the column-integrated water cloud distribution for a tidally locked (TL) planet (2 R_⊕, 1.4 g_⊕) around an M3 star, for a fixed rotation period of 60 Earth days, and for various levels of ISR (the M3-YANG cases in Table A.1). The parameters used here for the star and the planet closely match those used in the standard simulations of Yang et al. (2013), except that it is assumed here that the entire water reservoir is initially vaporized in the atmosphere. While “cold start” GCM simulations in Yang et al. (2013), with a global surface liquid water ocean, predict the formation of low-altitude convective dayside clouds, our “hot start” GCM simulations, with a steam atmosphere and a dry surface, predict the formation of high-altitude stratospheric nightside clouds.

Figure 2b shows the horizontal distribution of clouds and winds in the TL coordinates (Koll & Abbot 2015). In the TL coordinate system, the TL latitude is equal to 90° at the substellar point and −90° at the antisubstellar point. Figure 2b was derived by converting Fig. 2a (which is similar to Fig. 1, top row, right column) using the Python package of Koll (2020). The TL coordinates highlight well the strong dichotomy in cloud distribution between the dayside and the nightside of the planet, with two maxima located near the two wind gyre structures (near −120° in longitude, ±60° in latitude).

The cloud formation mechanism at play has already been identified and discussed for the cases of early Earth and Venus in Turbet et al. (2021), and we find that the same mechanism also applies here to TL exoplanets. In short, water vapor which is the main component of the atmosphere efficiently absorbs the incident stellar flux. This shortwave heating breaks moist convection on the dayside and thus prevents the formation of convective clouds. This shortwave heating also drives strong winds in the stratosphere that transport air parcels from the dayside to the nightside, which is qualitatively similar to the results of Fujii et al. (2017). As the air parcels move to the nightside, they are cooled by longwave emission (i.e., thermal emission to space), which reduces their temperature to the saturation temperature of water vapor. This leads to large-scale condensation, and thus to the formation of stratospheric clouds preferentially located on the nightside and at the poles. The 3D distribution of water clouds and winds is depicted in Figs. 2b,c using the TL coordinates. We note that we observe in Figs. 2a–c a slight decrease in the cloud cover at the antisubstellar point. This is particularly visible in the TL coordinates, with a decrease in clouds between −60° and −90° TL latitudes. This reduction in cloud content is associated with a downward movement of air parcels near the antisubstellar point (see the wind field in Fig. 2c). This subsidence produces adiabatic warming, thus reducing the amount of clouds in this region, which is qualitatively similar to the results of Charnay et al. (2021). We also note that the asymmetry between the northern and southern hemispheres in the cloud distribution (see, e.g., Fig. 2a) is due to temporal variability. The asymmetry in fact vanishes when averaging over a longer period of time (e.g., 5 Earth years, instead of 250 Earth days, as illustrated in Fig. 4).

The preferential accumulation of clouds on the nightside and lack of clouds on the dayside has two consequences: (1) lowering the bond albedo (see Fig. 1, middle row); and (2) producing a strong greenhouse effect on the nightside (see Fig. 1, bottom row). Figure 1 (bottom row) illustrates that the planetary thermal radiation displays a strong day-to-night amplitude. First, this is due to the greenhouse effect of nightside clouds and second, it is due to the fact that a significant fraction of the ISR is absorbed high in the atmosphere, which is qualitatively similar to the results of Wolf et al. (2019). This shortwave heating heats up the atmosphere on the dayside, which thermally re-radiates to space, further increasing the day-night outgoing thermal radiation (OTR) contrast. This 3D effect is the direct consequence of the distribution of the stellar flux, combined with the strong direct absorption of this flux by water vapor.

As a consequence, we find the same result as in Turbet et al. (2021) that, depending on the initial state (all water initially in the form of vapor in the atmosphere versus all water initially condensed at the surface), there is a wide range of ISR values for which climate end states are different. This hysteresis is illustrated in Fig. 3 for the surface temperature (Fig. 3a), planetary albedo (Fig. 3c), upper atmosphere water vapor mixing ratio (Fig. 3b) and greenhouse effect (Fig. 3d). The greenhouse effect G is defined as the difference between the surface temperature and the equilibrium temperature of the planet $$G=T_{\text{surf}}-{\left(\frac{\text{OTR}}{\sigma }\right)}^{\frac{1}{4}},$$(3)

with T_surf the global mean surface temperature, OTR the global mean thermal radiation to space, and σ the Stefan-Boltzmann constant. Our new simulations (Fig. 3, red branches) are compared to the standard simulations of Yang et al. (2013; Fig. 3, blue branches) in the standard case of a 2 R_⊕,1.4 g_⊕ planet synchronously rotating (fixed rotation period of 60 Earth days) around an M3 star.

Figure 3a shows the evolution of the surface temperature as a function of ISR, for two types of initial states. The arrows indicate an estimate of the ISR threshold limit at which the radiative TOA budget becomes unstable, leading to a runaway greenhouse effect (for the cold start sequence) and runaway water condensation (for the hot start sequence). For the hot start sequence, the arrows are positioned halfway between the last simulation (or lowest ISR) that does not lead to surface water condensation and the first simulation (or highest ISR) that does lead to surface water condensation (i.e., the simulation in which rainwater hits the surface). For the cold start sequence, the arrows are positioned halfway between the last simulation (or highest ISR) that avoids the runaway greenhouse effect and the first simulation (or lowest ISR) that goes into a runaway greenhouse effect (Yang et al. 2013). As already identified in Turbet et al. (2021), near the water condensation limit the number of dayside clouds grows when reducing the ISR, leading to an increase in bond albedo (see Fig. 3c, red branch). As soon as the water condensation limit is attained, the dayside becomes largely covered by clouds, which makes the bond albedo jump. This produces a strong TOA radiative desequilibrium, which inevitably leads to the condensation of water vapor on the surface. This imbalance is maintained until almost all the water reservoir has condensed onto the surface. This cloud-driven positive feedback is what triggers the runaway water condensation.

Fig. 1 Water cloud column (top), albedo (middle), and thermal emission to space (bottom) horizontal maps, for an initially hot and steamy 2 R⊕ planet (with a 10 bar H2O + 1 bar N2 atmosphere) synchronously rotating (fixed rotation period of 60 Earth days) around an M3 star. The simulations shown here were forced to ISR (defined here and throughout the manuscript as the global mean average value, i.e., a factor of 4 lower than the stellar flux at the substellar point) ranging from 500 W m−2 or 1.47 F⊕ (right) down to 300 W m−2 or 0.88 F⊕ (left), and the outputs were averaged over 250 Earth days. The associated simulations names are M3-YANG-1, 2, 3, 4, and 5. For reference, insolation is 340.5 W m−2 on present-day Earth.

Fig. 2 Water cloud horizontal and vertical distributions. Water cloud Col. (a) in latitude vs. longitude coordinates and (b) tidally locked (TL) coordinates (b). Water cloud content vertical distribution (c) as a function of TL latitude (averaged over all TL longitudes). The plots were computed for an initially hot and steamy 2 R⊕ planet (with a 10 bar H2O + 1 bar N2 atmosphere) synchronously rotating (fixed rotation period of 60 Earth days) around an M3 star, and forced to an ISR of 500 W m−2 (simulation M3-YANG-1). The outputs were averaged over 250 Earth days. For reference, insolation is 340.5 W m−2 on present-day Earth. Wind fields are superimposed on the plots. Winds are plotted at a pressure of ~50 Pa in panels a and b, and averaged over all TL longitudes in panel c.

3.2 Water condensation around different types of stars

We conducted a grid of simulations to explore the conditions for water condensation depending on the type of host star (from Sun-like stars to late M stars; see cases SUN, K5, M3, Pcen, and T1 in Table A.1) for an Earth-like planet (1 M_⊕, 1 R_⊕) synchronously rotating (except in the case of the Sun) around its host star. Again, we find that in all the simulations, water clouds preferentially form on the nightside and at the poles. This is illustrated in Fig. 4 (top row), which shows the column-integrated water cloud distribution for the planet depending on the type of host star. Here again, the albedo map (Fig. 4, middle row) and the OTR map (Fig. 4, bottom row) are correlated and anti-correlated with the distribution of clouds, respectively.

As discussed in the previous section, we observe again a slight decrease in the cloud cover at the antisubstellar point (Fig. 4, top row), associated with subsidence warming. We note that the cloud patterns are broadly similar between all types of host star, in particular all those assuming a synchronous rotation. In all the simulations the cloud content is highest near the two wind gyre structures, the coldest parts of the atmospheres, although their exact location and extension is affected by the parameters of the simulations (in particular, the rotation rate; see Carone et al. 2014, 2015, 2016).

While the thermal emission maps (Fig. 4, bottom row) are qualitatively similar between all the simulations, we observe a much higher dayside emission for the lowest mass stars (TRAPPIST-1, Proxima Cen) than the other cases. This is first due to the bond albedo being lower (see Fig. 4, middle row; see also discussions below), but also, and above all, to the fact that low-mass stars emit more near-IR flux that gets absorbed higher up in the atmosphere (Fujii et al. 2017) and directly re-emitted to space on the dayside (Wolf et al. 2019).

For each type of star we used the GCM simulations to build their water condensation sequences, illustrated in Fig. 5. Figure 5a shows the evolution of the surface temperature as a function of ISR, for all the types of host stars. The dashed arrows indicate the ISR threshold limit at which the radiative TOA budget becomes imbalanced, leading to runaway water condensation.

There are several noticeable points in Fig. 5. First, we find that the surface temperatures are lower for cooler stars even though more ISR is absorbed by the planet (as illustrated in Fig. 5b showing the bond albedos). This is somewhat counterintuitive as cooler stars have been shown to produce consistently warmer climates on planets with Earth-like atmospheres (Kopparapu et al. 2013; Shields et al. 2013; Godolt et al. 2015; Wolf et al. 2019). We discuss this effect in more detail below.

Second, the bond albedo decreases for cooler stars, which was expected (see, e.g., Kopparapu et al. 2013; Shields et al. 2013). Moreover, the bond albedo increases when the ISR decreases. This effect, which is discussed in Sect. 3.1, is what drives the sharp radiative budget deficit at the TOA, and thus what defines the runaway water condensation for all types of stars. As the ISR reaches the highest values we explored, Fig. 5b shows that the bond albedos seem to converge asymptotically toward the results from 1D cloud-free calculations (here taken from Kopparapu et al. 2013). This is due to the fact that the dayside of the planet is getting warmer and therefore has fewer clouds.

Third, the water condensation limit ISR threshold decreases monotonically with decreasing stellar effective temperature. In other words, this indicates that for a given ISR, water condensation is more difficult to attain for planets around M stars than around Sun-like stars. The implications of this are discussed in Sect. 3.4.

To understand in more detail why surface temperatures decrease (in thick H₂O-rich atmospheres) for cooler stars, we show in Figs. 6a–c the vertical profiles of the atmospheric temperatures, cloud mixing ratios, and absorbed stellar flux. While atmospheric temperatures are higher for cooler stars in the upper atmosphere (in the stratosphere, typically above 10⁴ Pa), this is reversed below in the troposphere where atmospheric temperatures become lower for cooler stars. Following Selsis et al. (2023), we found that the absorption of the ISR by water bands, combined with Rayleigh scattering of H₂O molecules, is so efficient that not enough shortwave flux reaches the lower atmosphere to sustain convection. This is illustrated in Fig. 6c, which shows that the absorbed stellar radiation drops in the low atmosphere due to most of the flux being absorbed higher in the atmosphere. Selsis et al. (2023) explored and discussed this effect in detail, and we verify here that the same effect appears in all our GCM simulations of thick H₂O-rich atmospheres. The clouds have little effect on the temperature structure in the low atmosphere, given that the dayside is almost entirely depleted in clouds (Fig. 4, top row). In other words, from the point of view of the shortwave radiative transfer calculations, the atmosphere behaves similarly to the cloud-free case, thus validating the approach of Selsis et al. (2023).

Fig. 3 Hysteresis loops and conditions of ocean formation for an initially hot and steamy 2 R⊕ planet (with a 10 bar H2O + 1 bar N2 atmosphere) synchronously rotating (fixed rotation period of 60 Earth days) around an M3 star (simulation M3-YANG-1). The figure depicts the evolution of surface temperature (top left), stratospheric water mixing ratio (top right), bond albedo (bottom left), and greenhouse effect (bottom right) as a function of the incoming stellar flux. The blue branches correspond to simulation results (Yang et al. 2013) that assume water is initially condensed on the surface. The red branches (this study) correspond to initially hot and steamy simulation results that assume water is initially in vapor form in the atmosphere.

Fig. 4 Water cloud column (top), albedo (middle), and thermal emission to space (bottom) horizontal maps, for an initially hot and steamy 1 R⊕ planet (with a 10 bar H2O + 1 bar N2 atmosphere) orbiting around different types of stars: TRAPPIST-1 (T★ ~ 2600 K), Proxima Centauri (T★ ~ 3050 K), an M3 star (T★ ~ 3400 K), a K5 star (T★ ~ 4400 K), and the Sun (T★ ~ 5780 K). The associated simulations names are T1-1, Pcen-1, M3-1, K5-1, and SUN-1, respectively. The simulations shown here were forced to an ISR of 500 W m−2 (i.e., 1.47 F⊕), and the outputs were averaged over 100 Earth days. A synchronous rotation for all the simulations was assumed, except for those of the Sun. In this case the results were taken directly from Turbet et al. (2021): rotation period of 1 Earth day, orbital period of 365 Earth days, and obliquity of 23.5°. For this simulation only, the maps were calculated in the heliocentric frame (i.e., keeping the subsolar point at 0° longitude and 0° latitude), and using an average of five Earth years as in Turbet et al. (2021).

3.3 Impact of the water content

We conducted a series of simulations to explore the effect of the water mass fraction in our simulations (M3, M3-3bar, M3-1bar, and M3-0.3bar from Table A.1), and how it potentially affects the water condensation limit. Figure 7 shows the water condensation sequences for H₂O atmospheric pressures from 0.3 to 10 bar (all simulations assume 1 bar of N₂ as a background gas).

We first find that the surface temperatures increase with water content, as already seen in Turbet et al. (2021), and as expected from the increasing greenhouse effect of water vapor. We note however that for planets orbiting very cool stars, the surface temperature increases very little with H₂O partial pressure because the stellar flux in the lower atmosphere is too low to maintain convection, and thus the thermal profile reaches an isotherm (see Sect. 3.5 where this is discussed in the context of TRAPPIST-1b).

The bond albedo sequence looks very similar for all H₂O partial pressures explored, except that at $P_{{\text{H}}_2\text{O}}=0.3$ bar. For this case the water vapor content is low enough to reduce the fraction of ISR absorbed by the atmosphere, which leads to an increase in bond albedo (as a reminder, the surface albedo is fixed at 0.2 in these simulations). A similar behavior is seen in the 1D cloud-free simulations of Kopparapu et al. (2013, see, e.g., their Fig. 6a) where the bond albedo for an Earth-like planet around an M-type star decreases as surface temperature increases, due to the increasing water vapor content. Moreover, given that less ISR is absorbed on the dayside by water vapor, the dayside cloud positive feedback (producing the runaway water condensation introduced in previous subsections) is even more efficient. This is visible in Fig. 7b revealing that the bond albedo increases more rapidly for the $P_{{\text{H}}_2\text{O}}=0.3$ bar simulation when decreasing ISR. As a result, the water condensation limit occurs at slightly higher ISR for the $P_{{\text{H}}_2\text{O}}=0.3$ bar simulation (see Fig. 7a) than for the other H₂O partial pressures explored. It is not clear if this trend would remain for even lower H₂O atmospheric mixing ratios and/or H₂O partial pressures, given that the saturation temperature also decreases with partial pressure. This deserves to be explored in more detail, possibly with some interesting implications for the case of early Venus, and more generally for the case of land exoplanets (Kodama et al. 2015; Ding & Wordsworth 2020).

The fact that the evolutions of the bond albedo as a function of ISR (see Fig. 7b) look very similar for all GCM simulations with $P_{{\text{H}}_2\text{O}}\ge 1$ bar is a good indication that we are capturing well the climate regime in which H₂O is a dominant gas (starting from $P_{{\text{H}}_2\text{O}}=1$ bar, up to pressures as high as desired). Although the amount of water has an impact on the surface temperature (see Fig. 7a), for planets more irradiated than the water condensation limit the surface temperatures are always higher than that of the saturation temperature of the assumed H₂O partial pressure. The fact that the $P_{{\text{H}}_2\text{O}}=10$ bar simulation produces surface temperatures even higher than the critical temperature of water (T_crit = 647 K) indicates that for $P_{{\text{H}}_2\text{O}}\ge 10$ bar, the water condensation limit would remain unchanged. More specifically, H₂O-rich atmospheric simulations at different total pressures have thermal profiles that overlap significantly (see Fig. 13c, in the context of TRAPPIST-1b). Increasing the surface pressure to higher values would likely produce an extension of the thermal profile, while keeping a similar profile at low atmospheric pressures.

For planets orbiting ultra-cool stars, Selsis et al. (2023) showed that the thermal profile in the lower atmosphere can converge to an isotherm at a temperature that is lower than that of the critical point of water, potentially permitting the co-existence of a steam atmosphere and liquid water oceans for planets more irradiated than the runaway greenhouse limit. Although we did not explore this possibility in detail, our simulations for TRAPPIST-1 inner planets (see Figs. 12 and 13) support this result. However, this result might be restricted to a purely theoretical experiment since a tiny internal heat flux is enough to increase the surface temperatures above the critical temperature of water (Selsis et al. 2023; see also Fig. 13b in the context of TRAPPIST-1b).

To summarize, the sensitivity experiment described in this subsection shows that the choice of $P_{{\text{H}}_2\text{O}}=10$ bar for our GCM calculations is a relevant choice to represent all the cases where water is a major component of the atmosphere. Additional radiative active gases (e.g., CO₂) could change the TOA radiative budget, likely producing more warming (see Turbet et al. 2021), but this deserves to be explored quantitatively (see discussions in Sect. 4).

Fig. 5 Conditions of ocean formation for initially hot and steamy 1 R⊕ planets (with a 10 bar H2O + 1 bar N2 atmosphere) orbiting around different types of stars: TRAPPIST-1 (T★ ~ 2600 K), Proxima Centauri (T★ ~ 3050 K), an M3 star (T★ ~ 3400 K), a K5 star (T★ ~ 4400 K) and the Sun (T★ ~ 5780 K). The associated simulations names are T1-1 to 8, Pcen-1 to 7, M3-1 to 7, K5-1 to 6 and SUN-1 to 6, respectively. The figure depicts the evolution of surface temperature (top) and bond albedo (bottom) as a function of the incoming stellar flux. We made the assumption of a synchronous rotation for all the simulations, except for those of the Sun (where we took the results directly from Turbet et al. (2021) in which rotation period is equal to 1 Earth day, orbital period is equal to 365 Earth days, and obliquity is equal to 23.5°). Horizontal dashed lines indicate the bond albedos of planets in a runaway greenhouse effect, obtained from 1D cloud-free calculations (Kopparapu et al. 2013). The dashed arrows indicate the ISR threshold limit at which the radiative TOA budget becomes imbalanced, leading to a runaway water condensation.

3.4 Consequences for the habitable zone

We synthesized our results to define a water condensation zone (WCZ), which we propose to use to redefine the inner boundary of the HZ. We recall that the ISR at which the surface water condensation starts is necessarily less than or equal to the ISR required to trigger the runaway greenhouse effect. It is therefore a more restrictive criterion. Figure 8 highlights in bold red the water condensation limit calculated with GCM simulations for five types of host stars (Sun-like, K5, M3, M5.5, M8 stars). Planets receiving lower ISR than this limit are within the WCZ of their host star, meaning ocean formation by condensation is permitted. To put this result in context, we added various inner edge boundaries of the HZ, based on the results of 1D calculations (Kopparapu et al. 2013) and 3D calculations for slow (Kopparapu et al. 2017; Way & Del Genio 2020) and fast rotators (Leconte et al. 2013a; Wolf & Toon 2015). We note that the inner edge of the HZ is at much higher ISR for slowly rotating planets because GCM simulations predict for these planets the formation of highly reflective dayside convective clouds (Yang et al. 2013, 2014) that significantly expand the HZ toward higher ISR.

To derive the water condensation limit, we followed two steps. In the first step, based on the results presented in Sect. 3.2, we estimated for the five types of host stars (Sun-like, K5, M3, M5.5, M8 stars) five estimates of the water condensation limit (Fig. 8, red brackets) bracketed with the ISR of two GCM simulations: the ISR of the first GCM simulation reaching runaway water condensation, and the ISR of the last GCM simulation not reaching runaway water condensation. In the second step we fitted the five brackets with the empirical HZ inner edge of Kopparapu et al. (2013) shifted by 0.065 F_⊕ (i.e., 0.065 the insolation on Earth). This corresponds roughly to a difference of 22 W m⁻² between the inner edge of the HZ calculated in Kopparapu et al. (2013) and our water condensation limit. We did not attempt to fit the points with a more sophisticated curve since the Kopparapu et al. (2013) fit with a shift is sufficient to satisfactorily match the results of our GCMs simulations. The condensation curve looks broadly similar to the HZ inner edge of Kopparapu et al. (2013) because low-mass stars emit more flux in the near-IR, which produces more shortwave warming that delays water condensation. This is the same effect that makes the runaway greenhouse effect occur at lower ISR for planets around low-mass stars (Kopparapu et al. 2013). Compared to 1D cloud-free calculations, nightside clouds add an additional greenhouse effect, which is mainly responsible for the shift of the HZ inner edge toward lower ISR. The net effect of clouds is roughly similar across the various types of host stars (see Fig. 9).

It is important to note that the ISR that matters for determining whether a planet can condense its primordial water vapor reservoir into oceans is the ISR the planet received at the beginning of the main sequence phase of its host star. This is commonly referred to as the ZAMS. For the types of stars we are studying here, the ZAMS corresponds to the moment when the stellar luminosity is minimal, while the planets are already fully formed. In Fig. 8, we indicated the ISR of all the currently known planets and exoplanets (with a mass ≤ 5 M_⊕ and/or a radius ≤ 1.6 R_⊕) in or near the inner boundary of the HZ are shown (black squares), along with the ISR they received at the ZAMS (gray brackets). The planets are based on the following publications: HD85512b (Pepe et al. 2011), GJ667Cc (Bonfils et al. 2013; Anglada-Escudé et al. 2013), τ Ceti e (Tuomi et al. 2013), Kepler-138e (Rowe et al. 2014; Piaulet et al. 2023), Kepler-442b and Kepler-438b (Torres et al. 2015), K2-3d (Crossfield et al. 2015), Kepler-445d (Muirhead et al. 2015), Proxima b (Anglada-Escudé et al. 2016), K2-72c and e (Dressing et al. 2017), TRAPPIST-1c, d and e (Gillon et al. 2017), Kepler-1652b (Torres et al. 2017), Wolf 1061c, Luyten b (Astudillo-Defru et al. 2017), Kepler-1649b and c (Angelo et al. 2017; Vanderburg et al. 2020), Ross 128b (Bonfils et al. 2018b), GJ 1132c (Bonfils et al. 2018a), Teegarden b (Zechmeister et al. 2019), L98-59e (Kostov et al. 2019), TOI-700d and e (Gilbert et al. 2020, 2023), TOI-1266c (Demory et al. 2020), GJ 1061c and d (Dreizler et al. 2020), LP-890-9c (Delrez et al. 2022), Ross 508b (Harakawa et al. 2022), GJ 1002b (Suárez Mascareño et al. 2023), Wolf 1069b (Kossakowski et al. 2023), TOI-715b (Dransfield et al. 2023). While the previous limits of the inner edge of the HZ (Kopparapu et al. 2013, 2017; Leconte et al. 2013a; Wolf & Toon 2015; Way & Del Genio 2020) should be compared with the ISR that planets receive today, the water condensation limit (Fig. 8, bold red line) should be compared with the ISR that planets received at the ZAMS.

The ISR of planets at the ZAMS were calculated by first using stellar properties, including age estimates and associated uncertainties based on up-to-date publications, and then by calculating the variation in the luminosity of the star between the ZAMS and its current measured age. The gray brackets in Fig. 8 correspond to the ±1σ age estimates. For planets whose age constraints are either very loose or non-existent (mostly for low-mass stars), we arbitrarily used age estimates between 1 and 10 Gyr, which is the maximum value in the grid of stellar models of Baraffe et al. 2015). While ISR changes significantly with age (between now and the ZAMS) for planets around solartype and K-type stars, they change very little around low-mass (M-type) stars. Sun-like stars for instance have a luminosity that increases by about 40% during the first 4.5 Gyr of evolution in their main sequence. Ultra-cool stars like TRAPPIST-1 have an extended PMS phase during which their luminosity can decrease by more than two orders of magnitude within the first gigayear of evolution (Bolmont et al. 2017). However, once they reach their minimum luminosity (at the ZAMS), then their luminosity increases very little on a timescale less than or equal to the age of the universe.

Figure 8 clearly illustrates the findings of Turbet et al. (2021), by showing that Venus never received a low enough insolation to condense its oceans, and also that the Earth was able to condense its oceans thanks to the fact that the Sun’s luminosity was fainter in the past than it is today (known as the faint young sun opportunity). It also shows that the standard water condensation curve (bold red line in Fig. 8) can be compared directly to the ISR received today by planets orbiting around M stars because their luminosity evolves very little on the main sequence. We predict that several planets (e.g., K2-3d, TOI-700d, TOI-700e, K2-72e, Luyten b, LP890-9c) that are located well inside the HZ calculated using 3D models for tidally locked planets (Kopparapu et al. 2017), may in fact never have been able to form surface oceans. Interestingly, for TOI-700d this prediction depends on the age of the planet, which controls the minimum ISR it received at the ZAMS. This result highlights the importance of accurately measuring the age of planetary systems (at least, for planets around early M and more massive stars) to assess their ability to form oceans in the first place.

In parallel, and to facilitate the usability of our results by the community, we show in Fig. 10 a series of water condensation curves at different ages. The curves indicate, for a planet of age X (here with X = 0.5, 1, 4.5, and 10 Gyr), that if that planet is to the left of the age X curve, then it has never been able, or will never be able, to condense its water at the surface. The numerical values of the curves are provided in Table 2. These curves were first calculated from the standard water condensation curve (black dashed line in Fig. 10, which is the same as the bold red line in Fig. 8) and then multiplied by the ratio $\frac{L_{\star ,\text{at\hspace{0.17em}age\hspace{0.17em}=\hspace{0.17em}X}}}{L_{\star ,\text{at\hspace{0.17em}the\hspace{0.17em}ZAMS}}}$. This ratio was determined using the grid of stellar models of Baraffe et al. (2015). We note that the most massive stars in Fig. 10 are too short-lived for the 4.5 and 10 Gyr limits to be defined. We also note that to evaluate the potential habitability of a planet, the curves in Fig. 10 should be used in tandem with the traditional HZ (Kopparapu et al. 2013, 2017). In the case of planets orbiting highly evolved stars (e.g., the Earth in more than 1 billion years), the inner edge of the HZ may be more restrictive than the water condensation limit. In other words, the planet may have formed its first oceans by condensation, but then lost them several billion years later when a runaway greenhouse effect was triggered by the increase in the star’s luminosity (Leconte et al. 2013a; Wolf & Toon 2015).

In summary, using a grid of GCM simulations, we have defined for the first time a water condensation zone (WCZ), that we propose to use as a new standard to define the inner edge of the HZ. It is fortunate that this limit lies very close, albeit at lower ISR, to the 1D cloud-free calculations of Kopparapu et al. (2013), which are the most widely used definition in the exo-planet community. This is mainly because water vapor profiles are similar between 1D and 3D climate models, and because the absence of dayside clouds makes the atmosphere behave as if it were cloud-free from the point of view of shortwave heating. There is also a significant impact of the greenhouse effect of nightside clouds, and of the temporal evolution of stellar luminosity, which both quantitatively affect the water condensation limit (e.g., compared to Kopparapu et al. 2013).

An important point to note is that the question of the ability of planets to condense their primordial atmospheric water vapor reservoir is in fact even more constraining for planets orbiting M stars (than Sun-like stars) since the level of ISR they receive is much higher during their PMS, which can last up to 1 Gyr (Baraffe et al. 2015) for the lowest mass stars. To evaluate the ability of a planet to condense its primordial water reservoir, the water condensation limit should be compared to the minimum ISR received by the planet over the course of its evolution (in general, at the beginning of the main sequence phase, also known as the ZAMS), and not to its currently observed ISR. As a consequence, unlike planets orbiting solar stars, planets around M stars cannot benefit from a lower ISR period at the beginning of their evolution (Turbet et al. 2021) favorable for water condensation. This is clearly seen in Fig. 10 by looking at the difference between the water condensation limit at the ZAMS (dashed black lines) and the water condensation limit at several ages (solid colored lines).

Fig. 6 Vertical profiles (global mean, envelope defined by global maximum and minimum values) of atmospheric temperatures (left), water cloud mixing ratio (middle), and absorbed stellar fluxes (right). The profiles were computed for an initially hot and steamy 1 R⊕ planet (with a 10 bar H2O + 1 bar N2 atmosphere) orbiting around different types of stars: TRAPPIST-1 (T★ ~ 2600 K), Proxima Centauri (T★ ~ 3050 K), an M3 star (T★ ~ 3400 K), a K5 star (T★ ~ 4400 K), and the Sun (T★ ~ 5780 K). The curves were drawn from the same GCM simulations as in Fig. 4. The associated simulations names are T1-1, Pcen-1, M3-1, K5-1, and SUN-1, respectively.

Fig. 7 Conditions of ocean formation for initially hot and steamy 1 R⊕ planets synchronously rotating around an M3 star (T★ ~ 3400 K). The figure depicts the evolution of surface temperature (top) and bond albedo (bottom) as a function of the incoming stellar flux, for four different bulk atmospheric compositions: 10 bar H2O + 1 bar N2 (M3), 3 bar H2O + 1 bar N2 (M3-3bar), 1 bar H2O + 1 bar N2 (M3-1bar), and 0.3 bar H2O + 1 bar N2 (M3-0.3bar).

Fig. 8 Various inner edge boundaries of the habitable zone (HZ) for different types of stars, based on the results of 1D calculations (Kopparapu et al. 2013) and 3D calculations for slow (Kopparapu et al. 2017; Way & Del Genio 2020) and fast rotators (Leconte et al. 2013a; Wolf & Toon 2015). The black solid and dashed lines were calculated based on 1D climate calculations, for the runaway greenhouse effect and the moist greenhouse effect, respectively (Kopparapu et al. 2013). The gray solid and dashed lines were calculated based on 3D climate calculations of tidally locked planets, for the runaway greenhouse effect and the moist greenhouse effect, respectively (Kopparapu et al. 2017). The red bold line (and red brackets) indicates the newly calculated water condensation limit. All the currently known planets and exoplanets (with a mass ≤ 5 M⊕ and/or a radius ≤ 1.6 R⊕) in or near the inner boundary of the HZ are shown (black squares). Also shown (gray brackets) are the incoming stellar radiation (ISR) they received at the ZAMS based on their age estimates, using the grid of stellar models of Baraffe et al. (2015). The water condensation limit (red bold line) should be compared with the ISR that planets received at the ZAMS (gray brackets).

3.5 Application to TRAPPIST-1 planets and their observability with JWST

3.5.1 Water condensation and clouds

We extended our GCM simulations of H₂O-rich atmospheres to known exoplanets, in particular to TRAPPIST-1 planets (see Table A.1 for a detailed list). Figure 11 shows the cloud distribution, bond albedo, and OTR for TRAPPIST-1 inner planets (based on the stellar and planet properties provided in Agol et al. 2021, following Mann et al. 2019 and Ducrot et al. 2020). Again, clouds are located mostly on the nightside, high in the atmosphere. It is clearly visible that the cloud horizontal extension reduces with increasing ISR. The extension of clouds has impacts on the observables, in particular transit spectra, which we discuss in the following subsections.

For the three TRAPPIST-1 inner planets, in none of the cases that we have explored is the water vapor able to condense on the surface, and thus to form oceans. TRAPPIST-1d, which lies very close to the runaway greenhouse limit (and could be either in or out, depending on subtle cloud feedbacks; see Wolf 2017; Turbet et al. 2018), is unable to reach runaway water condensation and should thus lie outside the HZ. Simulations of TRAPPIST-1e (and Proxima b, for which we also performed a simulation as it is a prime target for the JWST and the ELTs) systematically lead to the quasi-complete condensation of the water vapor reservoir onto the surface. This confirms what had been calculated in Sect. 3.2 in the case of an Earth-like planet around TRAPPIST-1, that the water condensation limit is located somewhere between the orbit of TRAPPIST-1d and e. Given that TRAPPIST-1b, c, and d are outside the WCZ of their star, they are unlikely to be habitable planets. More generally, there is a significant number of small exoplanets that had been previously put inside or close to the HZ calculated with 3D models (Kopparapu et al. 2017) and that may never have condensed their primordial water reservoir. This includes TOI-700d (Gilbert et al. 2020), TOI-700e (Gilbert et al. 2023), LP-890-9c (Delrez et al. 2022), Teegarden b (Zechmeister et al. 2019), K2-72e (Dressing et al. 2017), and K2-3d (Crossfield et al. 2015).

We performed several sensitivity simulations to explore the thermal structure of the TRAPPIST-1 planetary atmospheres, which are summarized in Figs. 12 and 13. Again, our results confirm the 1D accelerated time-marching model results of Selsis et al. (2023, see a comparison in Extended Data Fig. 3 in Selsis et al. 2023) that the thermal profiles strongly diverge from that calculated from 1D inverse radiative-convective models. In particular, due to the lack of shortwave flux reaching the lower atmosphere and surface, the low atmosphere is much colder than predicted from the adiabatic profile. The water vapor content also has a limited impact on the thermal structure (see Fig. 13c). The internal heat flux, which is possibly high on TRAPPIST-1b due to tidal heating (Turbet et al. 2018), has a strong impact on the lower part of the atmosphere (see Fig. 13b) by transporting energy from the surface to the upper part of the atmosphere. This brings an extra heating term warming the atmospheric layers near the surface. This confirms with a 3D GCM the results of Selsis et al. (2023). Both the water content and internal heat fluxes (in the range of simulations explored here) have almost no impact on the temperature structure above 0.1 bar, and on the location and amount of water clouds.

The fact that the thermal profiles calculated with the GCM are very different from the 1D inverse climate calculations has several important consequences (Selsis et al. 2023), including the reduction of the lifetime of a magma ocean and modification of the mass–radius relationships. All the vertical profiles discussed in this section are used as input for radiative transfer calculations to compute various types of observables (transit spectra, eclipses, phases curves) described in the following subsections.

Fig. 9 Cloud forcings for initially hot and steamy 1 R⊕ planets (with a 10 bar H2O + 1 bar N2 atmosphere) orbiting around different types of stars: TRAPPIST-1 (T★ ~ 2600 K), Proxima Centauri (T★ ~ 3050 K), an M3 star (T★ ~ 3400 K), a K5 star (T★ ~ 4400 K), and the Sun (T★ ~ 5780 K). The associated simulations names are T1-1 to 8, Pcen-1 to 7, M3-1 to 7, K5-1 to 6 and SUN-1 to 6, respectively. The figure depicts the greenhouse effect (top) of clouds, the amount of incident radiation (middle) reflected back by the clouds (the more negative the value, the greater the reflected flux), and the net radiative effect (bottom) of clouds (positive values mean warming), as a function of the incoming stellar flux. In all these simulations, clouds produce a strong atmospheric warming.

Fig. 10 Various positions of the water condensation limit, depending on the age of the planetary system considered. The standard water condensation limit (black dashed line) is valid at the ZAMS. It can be used as is to evaluate the condensation of water on a planet, provided that the insolation of the planet is taken at the ZAMS (see, e.g., gray brackets in Fig. 8). The other water condensation limits (solid colored lines) indicate, for a planet of age X (here with X = 0.5, 1, 4.5, and 10 Gyr), that if that planet is to the left of the age X curve, then it has never been able, or never will be able, to condense its water at the surface. These curves were calculated from the standard water condensation curve (black dashed line), and shifted using the grid of stellar models of Baraffe et al. (2015). Numerical values of all the curves are provided in Table 2.

3.5.2 Transit spectra

We used the various GCM simulations of TRAPPIST-1b, c, and d to compute their expected signature in transit spectroscopy. Figure 14 summarizes our results for H₂O-dominated atmospheres that we compare, for context, to spectra for CO₂-dominated atmospheres, also calculated from dedicated GCM simulations (not shown).

First, we observe that the amplitude of water bands in the transit spectra is quite strong, up to 150 ppm for TRAPPIST-1b. We compared our spectra with those of Koll et al. (2019, calculated in the range 1.75–3 μm) for cloud-free H₂O-dominated atmospheres. Our cloud-free spectra (see Fig. 14b, dashed lines) look qualitatively similar. When correcting the amplitude of the transit depth features by the scale height ratio (about 35% change) using the latest mass measurement of Agol et al. (2021) (compared to that of Grimm et al. 2018), the cloud-free spectra of Koll et al. (2019) are also quantitatively similar to ours. The amplitude of water bands is significantly stronger than that of CO₂ bands (up to 75 ppm for TRAPPIST-1b), due mainly to a larger scale height (with a lower mean molecular weight, and a higher temperature for H₂O than CO₂). We note that the 10 bar CO₂ spectra look qualitatively similar to those of Lustig-Yaeger et al. (2019) and Koll et al. (2019, calculated in the range 1.75–3 μm), although the amplitudes of CO₂ absorption bands are lower. This is again most likely because we used the planet properties of Agol et al. (2021), while Lustig-Yaeger et al. (2019) and Koll et al. (2019) used old values (with a higher mass for TRAPPIST-1b) from Grimm et al. (2018). Figure 14a illustrates that the relative transit depths (i.e., the transit depth corrected from an arbitrary offset) are very similar for the range of atmospheric pressures simulated (from 1 to 10 bar for H₂O, which confirms the results of Koll et al. 2019; from 0.1 to 10 bar for CO₂, which confirms the results of Morley et al. 2017). This is simply due to the fact that the atmospheric structure in the upper atmosphere (i.e., typically for pressures lower than 0.1 bar, which is also the radiative part of the atmosphere, to which transit spectroscopy is most sensitive) is weakly impacted by the total pressure (see, e.g., Fig. 13c). This result differs for CO₂ atmospheres from the predictions of Koll et al. (2019), which used a parameterization of the heat redistribution that produces day-to-night temperature variations that are significantly larger than what our GCM simulations predict. This impacts the temperature at the terminator, the atmospheric height scale, and thus the amplitude of the molecular bands in the transit spectra as a function of surface pressure.

Second, we observe that the impact of clouds varies strongly from one planet to another. This is illustrated in Fig. 14b, which compares the transit spectra of TRAPPIST-1b, c, and d with and without the effect of clouds (solid vs dashed lines). The horizontal extension of the clouds (toward the dayside, and thus the terminator regions) increases with decreasing ISR (see Fig. 11, top row). This is simply because the planet is warmer, which reduces the area of cloud stability. The contrast is particularly striking between TRAPPIST-1d, which has a significant amount of water clouds reaching the terminator, and the two inner planets, which lack clouds. For the TRAPPIST-1 system, we thus find that the ISR required for the water clouds not to significantly blur the transit spectra is somewhere between that of TRAPPIST-1c and d. More generally, we find, using the grid of simulations performed for this study, that removing most water clouds of the terminator typically requires an ISR at least twice that of the Earth.

TRAPPIST-1b, c, and d are transit spectroscopy targets of the JWST space-based observatory, starting from Cycle 1, with the NIRSpec and the NIRISS instruments. To prepare and anticipate the outcome of these observations, we computed error bars on the synthetic spectra by cumulating all transit observations of JWST Cycle 1 for the three TRAPPIST-1 inner planets (see Figs. B.1 and B.2). This includes a total of two transits of TRAPPIST-1b with NIRISS (Lim et al. 2021), four transits of TRAPPIST-1c with NIRSpec (Rathcke et al. 2021) and two with NIRISS (Lim et al. 2021), and two transits of TRAPPIST-1d with NIRSpec (Lafreniere 2017). Our spectra reveal that water vapor dominated atmospheres (i.e., without H₂ and He as the dominant gases) are theoretically within reach of JWST Cycle 1 observations using transit spectroscopy, but only for TRAPPIST-1b and c. For TRAPPIST-1d the H₂O bands are much weaker due to the effect of clouds forming at the terminator, and a detection would require a larger number of transits. Detection thresholds are provided in Table 3 for each of the atmospheric scenarios (H₂O- and CO₂-dominated atmospheres) and observations planned for the first cycle of JWST. We note that we did not account for the transit stellar contamination effect here (Rackham et al. 2018, 2023), which will be an issue (see review in Turbet et al. 2020a and references therein).

Fig. 11 Water cloud column (top), albedo (middle), and thermal emission to space (bottom) horizontal maps, for TRAPPIST-1b (left), TRAPPIST-1c (middle), and TRAPPIST-1d (right). The associated simulations names are T1b, T1c, and T1d, respectively. The simulations assume initially hot and steamy tidally locked planets (with a 10 bar H2O + 1 bar N2 atmosphere), and the outputs were averaged over 100 Earth days. The dashed white lines indicate the location of the terminator.

Fig. 12 Temperature vertical profiles for TRAPPIST-1b (left), TRAPPIST-1c (middle), and TRAPPIST-1d (right). The associated simulation names are T1b, T1c, and T1d, respectively. The solid lines are global mean, temporal average (over 100 Earth days) temperature profiles calculated from GCM simulations. The shaded lines are snapshots taken at different latitude and longitude points to show the range of profiles within a simulation. The dashed lines are calculated with a 1D inverse radiative-convective model (Turbet et al. 2019, 2020b) similar to those used in Kopparapu et al. (2013) to calculate the inner edge of the HZ.

Fig. 13 Temperature vertical profiles (a) comparing TRAPPIST-1b, c, and d (simulations T1b, T1c, and T1d), (b) several internal heat fluxes for TRAPPIST-1b (simulations T1b, T1b-Fgeo1, T1b-Fgeo5, and T1b-Fgeo25), and (c) several water partial pressures for TRAPPIST-1b (simulations T1b, T1b-3bar, T1b-1bar). The solid lines are global mean, temporal average (over 100 Earth days) temperature profiles calculated from GCM simulations. The dashed lines are calculated with a 1D inverse radiative-convective model (Turbet et al. 2019, 2020b) similar to those used in Kopparapu et al. (2013) to calculate the inner edge of the HZ.

3.5.3 Secondary eclipses

Using the same GCM simulations, we computed the expected signals in thermal emission of TRAPPIST-1b and c at the secondary eclipse. Figure 15 summarizes our results by showing the emission spectra of H₂O-dominated atmospheres compared, for context, with CO₂-dominated atmospheres and with airless planets. TRAPPIST-1b and c are secondary eclipse targets of the JWST, starting from Cycle 1, with the MIRI instrument. This includes a total of ten eclipses for TRAPPIST-1b (five in the F1280W filter in Lagage & Bouwman 2017, five in the F1500W filter in Greene et al. 2017), and four eclipses for TRAPPIST-1c in the F1500W filter (Kreidberg et al. 2021). We include in Fig. 15 the first results of this observational campaign in the F1500W filter for TRAPPIST-1b (Greene et al. 2023) and TRAPPIST-1c (Zieba et al. 2023). In Fig. 15 we note that we artificially reduced the observed secondary eclipses of TRAPPIST-1b and c by 9.8% to correct for the difference between the observed stellar flux in the MIRI F1500W filter (equal to 2.595 mJy, based on Greene et al. 2023 and Zieba et al. 2023, which provide an absolute stellar flux in very good agreement) and that predicted by the stellar spectra used for TRAPPIST-1 (using Kurucz 2005 stellar templates). We computed error bars near 15 μm (F1500W MIRI filter) on the synthetic spectra by using the observed values from Greene et al. (2023) for TRAPPIST-1b (± 99 ppm) and from Zieba et al. (2023) for TRAPPIST-1c (± 94 ppm), divided by the square root of the ratio between the synthetic Kurucz stellar spectrum of TRAPPIST-1 and the observed stellar flux in the MIRI F1500W filter (Greene et al. 2023; Zieba et al. 2023). For TRAPPIST-1b, we evaluated the error bars near 12.8 μm (F1280W MIRI filter) by using the error bars near 15 μm (F1500W MIRI filter) rescaled by the square root of the ratio of synthetic Kurucz stellar fluxes between the F1280W and F1500W filters.

We observe that the eclipse depth is strong for H₂O-rich atmospheres (up to 400 ppm and 250 ppm for TRAPPIST-1b and c, respectively, near 15 μm). We note that the spectra are very similar for the explored atmospheric pressure range (from 1 to 10 bar of H₂O, which confirms the results of Koll et al. 2019 for TRAPPIST-1b, calculated in the range 5-12 μm). This is because most of the planetary thermal emission comes from atmospheric layers at or above 1 bar (Boukrouche et al. 2021). The spectra are not sensitive here to the contribution of the surface thermal emission (Hamano et al. 2015) because the surface is not hot enough to emit in the visible. This result is likely to persist even for water content well above 10 bar (Selsis et al. 2023).

By comparison, emission spectra of CO₂-dominated atmospheres exhibit stronger variations when changing surface pressure, which confirms the results of Koll et al. (2019; calculated in the range 5–12 μm) and Ih et al. (2023; for TRAPPIST-1b). This is due to the fact that CO₂-dominated atmospheres have many more spectral windows (compared to H₂O), where the planetary thermal emission is thus sensitive to the surface temperatures (which depends on the greenhouse effect of CO₂, and thus on the total pressure of CO₂). It is clearly visible in Fig. 15 that CO₂ spectra look similar in the CO₂ absorption bands (e.g., near 15 μm), but look different in the CO₂ spectral windows.

Our emission spectra reveal for the two planets that CO₂, H₂O, and airless atmospheres leave distinct signatures, which are detectable by the MIRI instrument on board JWST (see Table 4). The differences between our simulated cases are particularly strong near 15 μm, where JWST Cycle 1 observations have been obtained for the two planets (Greene et al. 2023; Zieba et al. 2023). Our GCM-based calculations confirm that thick CO₂-dominated atmospheres are unlikely for TRAPPIST-1b and c (rejected at >6 σ for TRAPPIST-1b; rejected at >3 σ for TRAPPIST-1c).While a thick H₂O-dominated atmosphere is also unlikely for TRAPPIST-1b (rejected at >3 σ), it is compatible within ~1 σ with the F1500W MIRI observations of TRAPPIST-1c. This makes TRAPPIST-1c a promising candidate to test for the presence of water vapor and clouds using transit spectroscopy as well as secondary eclipses in other MIRI filters.

To isolate the contribution of the planet (and thus get rid of corrections on the stellar flux), and facilitate comparison with previous works (Lustig-Yaeger et al. 2019), we converted the eclipse depths into brightness temperatures. This was done by first computing a grid of secondary eclipse spectra at fixed brightness temperatures (to find the correspondence between eclipse depth and brightness temperature), and then by interpolating our thermal emission spectra (based on GCM simulations) using this grid. Figure 16 summarizes our results for TRAPPIST-1b and c, along with published observations (Greene et al. 2023; Zieba et al. 2023). We note that we recalculated the observed brightness temperatures and found values a few degrees lower than those in Greene et al. (2023) and Zieba et al. (2023) (${493}_{-26}^{+26}\text{\hspace{0.17em}K}$ instead of ${503}_{-27}^{+26}\text{\hspace{0.17em}K}$ for TRAPPIST-1b; ${376}_{-33}^{+31}\text{\hspace{0.17em}K}$ instead of ${380}_{-31}^{+31}\text{\hspace{0.17em}K}$ for TRAPPIST-1c). In our brightness temperature calculations our aim was to treat observations and models using the same procedure in order to make comparisons as fair as possible. We thus computed first, for a grid of brightness temperatures, a grid of eclipse depth spectra. We then converted the grid of spectra into a grid of F1500W MIRI eclipse depths by integrating the fluxes in the F1500W MIRI filter. Finally, we interpolated in this grid using the measured F1500W eclipse depth values (from Greene et al. 2023 for the TRAPPIST-1b grid and from Zieba et al. 2023 for the TRAPPIST-1c grid) to derive the observed brightness temperatures. For CO₂-dominated atmospheres (at least the 10 bar simulation), we find results that are qualitatively similar to those of Lustig-Yaeger et al. (2019). Most of the differences are likely due to the different atmospheric compositions used (e.g., Lustig-Yaeger et al. 2019 have lower brightness temperature near 8.5 μm likely due to the SO₂ absorption included), which makes it difficult to compare other effects (3D vs 1D for example). We note that although the measurements at 15 μm are incompatible for TRAPPIST-1b with almost all our models of CO₂ and H₂O atmospheres, this does not demonstrate the absence of an atmosphere that could have a composition very different from that explored in the present manuscript.

Fig. 14 Various transit spectra computed with PSG for the three innermost planets of the TRAPPIST-1 system calculated from 3D GCM simulations (assuming H2O-dominated composition, and also CO2-dominated for comparison). In the upper panel a, transit spectra of TRAPPIST-1b for various compositions and pressures are compared; in the lower panel b, transit spectra of TRAPPIST-1b, c, and d for two compositions only (10 bar H2O + 1 bar N2; 10 bar CO2) are compared. The dashed lines indicate H2O cloud-free atmospheres, while solid lines assume H2O cloudy atmospheres. The dotted lines indicate CO2 atmospheres (including the effect of CO2 ice clouds, although extremely limited here; instead, CO2 condensation was cut off for TRAPPIST-1b and c simulations with 100 mbar of CO2 as this would cause the atmosphere to collapse on the nightside.). These are relative transit depths; the curves are shifted with arbitrary offsets to compare them more easily.

3.5.4 Phase curves

Last but not least, using again the same GCM simulations, we computed their expected signals in thermal phase curves. Similar observations are planned for Cycle 2 of JWST (Program 3077, 55 h of observations, i.e., about 1 orbit of TRAPPIST-1c or 1.5 orbit of TRAPPIST-1b) with the MIRI instrument using the F1500W filter. This observation mode is a promising avenue to characterize the atmospheres of TRAPPIST-1b and c. Figure 17 summarizes our results by showing the thermal emission phase curves of H₂O-dominated atmospheres compared, for context, with CO₂-dominated atmospheres and with airless planets. For this figure, we integrated the thermal emission on the F1500W filter (centered around 15 μm) as it is the selected observing mode for the phase curve measurements to be obtained within JWST Cycle 2.

There are significant differences between the different types of atmospheres. The phase curves of H₂O-dominated atmospheres have relatively low day-to-night amplitudes. For TRAPPIST-1b, c, and d there is a factor of 1.7 to 2 (in the F1500W filter) between dayside and nightside emissions. We interpret this behavior, which is similar to what is observed in Wolf et al. (2019) for the incipient runaway case, as follows. At first order, the signal of the phase curve can be decomposed into two parts, a continuous signal (which persists at all phase angles) corresponding to the thermal emission of the Simpson-Nakajima limit (Nakajima et al. 1992; Goldblatt & Watson 2012), and a second source of emission localized on the dayside, which is produced by the re-emission of the ISR absorbed high in the atmosphere (Wolf et al. 2019). The first component is lower than the Simpson-Nakajima limit due to part of the emission being absorbed by nightside clouds. This is clearly visible by looking at the OTR maps in Figs. 1, 4, and 11 (bottom rows). The presence of nightside water clouds reduces the planet’s thermal emission on the nightside, leaving a visible signature on the phase curve. The effect is qualitatively comparable to that of silicate clouds on the phase curves of hot giant exoplanets (Parmentier et al. 2021). Clouds are more prominent, and thus have a stronger impact on the phase curves for planets receiving insolation close to the water condensation limit.

Again, we note that the phase curves look very similar for all the surface pressures explored (from 1 to 10 bar of H₂O). Interestingly, we observe a significant hotspot phase shift in the H₂O-rich simulations (which is less pronounced in the CO₂ phase curves). We interpret this hotspot shift as a consequence of the stellar flux being absorbed high in the atmosphere (combined with the effect of the mean molecular weight; see Zhang & Showman 2017) first because the shortwave heating combined with rotation rate produces strong westerly winds transporting heat from the dayside to the nightside. This is clearly illustrated in Fig. 2, which shows the strong horizontal winds in the stratosphere. This is also because the thermal emission comes from low-pressure atmospheric layers, where horizontal winds are strong. Interestingly, CO₂-dominated cases have a much less pronounced hotspot offset (near 15 μm, which is the strongest CO₂ absorption band), mostly because atmospheric winds are much weaker due to less shortwave flux being absorbed high in the atmosphere (where most of the OTR is emitted at 15 μm). The amplitude of the phase curve and hotspot offset depend on the wavelength considered. This is clearly illustrated in Figs. B.4, B.7, and B.9. We note that the simulated emission spectra of TRAPPIST-1b, c, and d are provided at different phase angles (including at secondary eclipse, discussed in the previous Section) in Figs. B.3–B.9.

Fig. 15 Emission spectra of TRAPPIST-1b (upper panel) and c (lower panel), calculated at the secondary eclipse. The spectra are calculated using GCM simulations of H2O-rich (blue to gray) and CO2-rich atmospheres (red to orange). Emission spectra are also computed for airless planets (in bold black). The 1σ error bars are calculated using the JWST Cycle 1 observation program (for more details, see Sect. 3.5.3).

4 Discussion

4.1 Comparisons with other GCMs

Three-dimensional simulations in Turbet et al. (2021) and in this study were performed with the same GCM, the Generic PCM. This raises the question of whether nightside cloud formation in H₂O-dominated atmospheres is also obtained in other 3D climate models. The multi-model validation (Fauchez et al. 2020, 2021) is indeed a crucial step to validate this cloud effect, which has potentially strong impacts on the habitability of planets (Turbet et al. 2021; Kasting & Harman 2021). The best illustration of this multi-model approach is the demonstration that substellar convective clouds appear in all models that have simulated a synchronously rotating aquaplanet, giving great confidence in this cloud feedback (see GCM intercomparisons in Yang et al. 2019b; Sergeev et al. 2022).

We thus surveyed the literature to look for evidence (or the absence of evidence) of a similar cloud feedback regarding the formation of nightside water clouds in other models. First of all, and to the best of our knowledge, there are no published results of 3D GCM simulations of post-runaway atmospheres (i.e., water-dominated atmospheres with an exhausted surface water reservoir), with the exception of Turbet et al. (2021). Boukrouche (2022, PhD manuscript) recently performed cloud-free GCM simulations of post-runaway water-dominated atmospheres and confirmed, using the ExoFMS-SOCRATES GCM, the preferential formation of clouds on the nightside, based on distance to saturation. Apart from this work, the best analogs we have can be found in numerical studies that simulated the climate of planets in incipient runaway greenhouse conditions (Kopparapu et al. 2017; Wolf 2017; Wolf et al. 2019; Kane 2022; Chaverot et al. 2023). Using the ExoCAM GCM, Kopparapu et al. (2017) first studied the climate and clouds of tidally locked rocky exoplanets transitioning into runaway greenhouse conditions. They found that as soon as a planet reaches the runaway greenhouse transition, its substellar convective cloud layer dissipates, drastically reducing the albedo and accelerating the runaway greenhouse effect. The same effect is seen with ExoCAM in Wolf (2017) with GCM simulations of TRAPPIST-1d. A snapshot of cloud distribution on TRAPPIST-1d as it evolves into the runaway greenhouse transition (see Wolf 2017, Fig. 4) not only reveals the absence of clouds on the dayside, but also a preferential accumulation of clouds on the nightside, with two local maxima near the wind gyre structures, similar to our findings. Similar cloud distributions are obtained again with ExoCAM (see Wolf et al. 2019, Figs. 4 and 5), for a broader range of host star types (from T_eff = 2600 K to 4500 K). Using ROCKE-3D GCM, Kane et al. (2018) also performed GCM simulations of the exoplanet Kepler-1649b in incipient runaway greenhouse conditions. The runaway greenhouse effect produces a sharp decrease in the bond albedo. We analyzed the netCDF simulation file, provided in Kane (2022); it reveals that water clouds are also preferentially located on the nightside in the upper part of the atmosphere.

In summary, our literature review reveals that the preferential formation of nightside clouds is also obtained in other GCMs (ExoCAM, ROCKE-3D, ExoFMS) as soon as they reach conditions where water vapor becomes a dominant gas in the atmosphere (called incipient runaway greenhouse conditions). This cloud feedback deserves to be studied and validated in the context of post-runaway water-dominated atmospheres with other GCMs (including the effects of water clouds).

Fig. 16 Same as Fig. 15, but with the emission converted into brightness temperature units.

Fig. 17 Thermal phase curves of TRAPPIST-1b (upper panel) and c (lower panel), integrated in the F1500W filter (centered around 15 μm), for various atmospheric scenarios: H2O-rich (blue to gray), CO2-rich atmospheres (red to orange), airless planet (bold black). The H2O-rich atmospheres also include 1 bar of N2.

4.2 Detecting the cloud feedback

Our simulations of synthetic spectra (see Sects. 3.5.2–3.5.4) reveal that JWST has the capability of probing water-rich atmospheres around TRAPPIST-1b, c, and d, either in transit spectra or in emission. Observations are particularly favorable for TRAPPIST-1b and c because the ISR received on the planets is high enough first to suppress clouds at the terminator, and thus to strongly increase the amplitude of water bands in transit spectra, and also to increase the temperatures, and thus the thermal emission of the planets. Measurements at 15 μm (MIRI F1500W filter) show, however, that a thick H₂O-dominated atmosphere is unlikely for TRAPPIST-1b (rejected at >3 σ), which leaves TRAPPIST-1c as the best candidate so far (compatible at ~1 σ). For TRAPPIST-1d the observations are more difficult not only because the planet is smaller and cooler, but also because water clouds blur the transit spectra. However, by accumulating enough transits (see Table 3), transit spectra could reveal the presence of water vapor. Together, comparative observations of the TRAPPIST-1 inner planets (planet c vs planet d) could be used to reveal the presence of nightside clouds in the case their composition turn out to be dominated by water vapor.

Although we used the TRAPPIST-1 planets as targets to compute observables, more easily accessible targets such as water-rich mini-Neptunes (Luque & Pallé 2022) could be used to detect this cloud effect, in particular by comparing observations of planets at different values of ISR. We identified a threshold ISR at around two times the insolation on Earth (the exact value depends on many factors that deserve to be explored in more detail) above which water clouds become significantly depleted at the terminator, significantly increasing the amplitude of absorption features in transit spectra. Based on the simulations we performed, we also identified another threshold ISR, around five times the insolation on Earth, above which water clouds entirely disappear. Although these ISR thresholds may vary from one planet to another, depending on its amount of water vapor, its size, and its type of star, for example, we looked for planets in this ISR range with the best observability metrics, using the Emission Spectroscopy Metrics (ESM) and Transmission Spectroscopy Metrics (TSM) of Kempton et al. (2018). In addition to the three inner planets of the TRAPPIST-1 (which are the planets with the best TSM and among the best ESM), we identified LP 791-18c (Crossfield et al. 2019; Peterson et al. 2023), TOI-270d (Günther et al. 2019), and LTT 3780c (Cloutier et al. 2020). These three planets have radii between 2.1 and 2.4 times that of the Earth, and their density is compatible with the presence of both water and hydrogen envelopes. If these planets happen to be enriched in water, then they could serve as a natural laboratory for testing the behavior of water clouds, using transit spectra, secondary eclipses, and phase curves. Previous simulations of K2-18b with the Generic PCM (Charnay et al. 2021) indeed show qualitatively similar cloud behavior for temperate mini-Neptunes with atmospheric composition of very high metallicity (> 300). At lower metallicity, clouds can form on the dayside essentially due to two factors (Charnay et al. 2015a,b); first, less water vapor induces less radiative heating on the dayside and less radiative cooling on the nightside, and second, the specific heat capacity decreases due to higher mixing ratios of hydrogen and helium, which further reduces the day-to-night temperature contrast.

Detecting the water cloud feedback using TRAPPIST-1 planets or any population of planets endowed with water-rich atmospheres would have profound consequences for our understanding of the habitability of exoplanets, but also of the Solar System planets (Earth vs. Venus).

4.3 Caveats

In this study we explored the ability of exoplanets to form their first oceans by simulating the condensation of the primordial water reservoir at the end of the magma ocean phase (Hamano et al. 2013; Lichtenberg et al. 2023).

In Turbet et al. (2021) and this new study, we performed several sensitivity studies, including varying the amount of CCN (Turbet et al. 2021), and varying the amount of water (this work; Turbet et al. 2021) and carbon dioxide (Turbet et al. 2021). However, to build confidence in these simulations, more sensitivity studies are needed. First, one can question the choice of parameterizations used in this work. Previous works (Yang et al. 2019b; Sergeev et al. 2020, 2022) have shown, for example, that the choice of the moist convection scheme has strong consequences on the properties of the substellar cloud layer predicted in the cold start simulations. Cloud formation in our hot start Generic PCM simulations is largely dominated by large-scale condensation (that is, water vapor condensation driven by large-scale air movements) and not by subgrid-scale condensation (driven by small convective cells, and parameterized using our moist convection scheme). This indicates that our results should be weakly sensitive to the choice of the subgrid moist convection parameterization, which is a strong indication that preferential nightside cloud formation is likely to be a robust mechanism across GCMs (see Sect. 4.1). However, this deserves to be tested in detail with other GCMs and with cloud-resolving simulations (Tan et al. 2021). Second, one can question the choice of atmospheric compositions used in this work. It has been shown that various atmospheric compositions could be expected at the end of the magma ocean phase, depending on the redox state (Bower et al. 2022). Modeling a wider range of atmospheric compositions, more representative of all the possible post-magma scenarios, is therefore a necessary step to validate our results.

Moreover, even if the primordial water reservoir never condensed, it is possible that late volcanic degassing or volatile delivery by cometary impacts could bring water that contributes to forming surface oceans (Moore & Cowan 2020). This would probably require the planet to have lost its atmosphere at some point, temporarily transforming the planet into a “land planet” (Kodama et al. 2015), otherwise the delivered water would probably remain vaporized. This scenario also deserves to be explored in detail, at least in part with the help of dedicated 3D numerical climate simulations.

5 Conclusions

In this study we used the 3D Generic planetary climate model to perform a large grid of simulations for H₂O-dominated atmospheres planets orbiting around a wide variety of main sequence stars. The simulations were designed to reproduce the conditions of early ocean formation on rocky planets due to the condensation of the primordial water reservoir at the end of the magma ocean phase (Hamano et al. 2013; Lichtenberg et al. 2023).

So far, all global climate model (GCM) simulations, assuming the entire water reservoir is initially fully condensed on the surface, have shown that low-altitude dayside convective clouds stabilize oceans even at high incoming stellar radiation (ISR; Yang et al. 2013, 2014, 2019a,b; Kopparapu et al. 2016, 2017; Way et al. 2016; Sergeev et al. 2020, 2022). In contrast, our new GCM simulations, assuming the entire water reservoir is initially fully vaporized in the atmosphere, show that high-altitude nightside stratospheric clouds prevent surface water condensation from happening even at low ISR. This behavior was already observed by Turbet et al. (2021) in the case of early Earth and Venus, and we thus extend it here to synchronously rotating planets around all types of main sequence stars (M, K, and G types). More generally, this confirms and extends the results of Turbet et al. (2021) that the formation of nightside clouds seems to occur independently of the rotation rate assumed. Future work could test to what rotation rate this result holds, given that there should be a threshold rotation rate above which Coriolis forces inhibit day-night stratospheric transport. This is particularly relevant for young planets, which are likely to have a much faster rotation rate than that determined by synchronous rotation. We propose to name this family of planets “mochi planets”, with the planetary surface represented by red bean paste, and the cloud layer represented by sticky rice from which someone has taken a bite.

Based on these results, we introduced a water condensation zone (WCZ), whose inner edge is located at lower ISR than the inner edge of the habitable zone (HZ). In other words, we found that the ISR required to form oceans by condensation is always significantly lower than that required to vaporize oceans. The contrast is particularly striking with all the recent 3D GCM results that assume most water is initially on the planetary surface (Yang et al. 2013, 2014; Kopparapu et al. 2016, 2017; Way et al. 2016). Interestingly, the inner edge of the WCZ lies quite close to (but still at lower ISR, due to the greenhouse effect of nightside clouds) the widely used inner edge of the HZ of Kopparapu et al. (2013), based on 1D cloud-free calculations. This is mainly because water vapor profiles are similar between 1D and 3D climate models, and because the absence of dayside clouds makes the atmosphere behave as if it were cloud-free from the point of view of shortwave heating.

We confirm with 3D GCM simulations the results of Selsis et al. (2023) that in thick H₂O-dominated atmospheres, the vertical thermal profile diverges from that calculated using 1D inverse radiative-convective models (Kopparapu et al. 2013; Turbet et al. 2019; Boukrouche et al. 2021). This has important consequences for the lifetime of the magma ocean period, and for the calculations of mass–radius relationships for irradiated planets (Turbet et al. 2020b).

Last but not least, we show, taking the inner planets of the TRAPPIST-1 system as an example, that the JWST has the capability to probe water-rich atmospheres. We identify a threshold insolation, around two times the insolation on Earth (the exact value depends on many factors), above which water clouds disappear from the terminator, significantly increasing the amplitude of absorption features in transit spectra. For the TRAPPIST-1 system, we predict that the transition occurs between TRAPPIST-1c and d. JWST observations using transit spectroscopy, but also secondary eclipses and thermal phase curves, can be used first to probe the existence of water-rich atmospheres on TRAPPIST-1 inner planets, and then to test the existence of the nightside cloud feedback. Although we used the TRAPPIST-1 planets as a test bed, other accessible targets such as water-rich mini-Neptunes (Luque & Pallé 2022) could be used to detect this cloud effect. This would have profound consequences for our understanding of the habitability of exo-planets, but also of that of the Solar System planets (Earth vs. Venus).

Acknowledgements

This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Grant Agreement No. 832738/ESCAPE. This work has been carried out within the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation. M.T. acknowledges the financial support of the SNSF. M.T. thanks the Gruber Foundation for its generous support to this research. M.T. acknowledges support from the Tremplin 2022 program of the Faculty of Science and Engineering of Sorbonne University. M.G., M.T. and F.S. acknowledge support from the BELSPO program BRAIN-be 2.0 (Belgian Research Action through Interdisciplinary Networks) contract B2/212/B1/PORTAL. This work was performed using the High-Performance Computing (HPC) resources of Centre Informatique National de l’Enseignement Supérieur (CINES) under the allocations No. A0060110391, A0080110391, A0100110391, A0120110391 made by Grand Équipement National de Calcul Intensif (GENCI). A total of ~1.3 M CPU hours were used for this project on the OCCIGEN supercomputer, resulting in ~2.3 t eq. of CO₂ emissions. T.J.F. acknowledges the support of the NASA GSFC Sellers Exoplanet Environments Collaboration (SEEC), which is funded by the NASA Planetary Science Divisions Internal Scientist Funding Model. J.L. and F.S. acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 679030/WHIPLASH), and from the french state: CNES, Programme National de Planétologie (PNP), the ANR (ANR-20-CE49-0009: SOUND), and in the framework of the Investments for the Future programme IdEx, Université de Bordeaux/RRI ORIGINS. The authors thank the Generic PCM team for the teamwork development and improvement of the model. This research has made use of NASA’s Astrophysics Data System. M.T. thanks the LMD Planeto team as well as Nadège Lagarde (LAB) and Xavier Delfosse (IPAG) for useful discussions. We thank the anonymous reviewer for the very insightful comments and suggestions, which improved the quality of the manuscript. The Generic-PCM (and documentation on how to use the model) can be downloaded from our SVN repository https://svn.lmd.jussieu.fr/Planeto/trunk/LMDZ.GENERIC/). More information and documentation on the model are available on http://www-planets.lmd.jussieu.fr. The GCM outputs of the TRAPPIST-1 planet simulations are provided here: https://zenodo.org/record/5627945. Other data underlying this article can be shared upon reasonable request to the corresponding author.

Appendix A Additional tables

We provide here in Tables A.1 the main list of 3D global climate model (GCM) simulations performed with the Generic PCM in this study.

Appendix B Additional observables

B.1 Transit spectra of TRAPPIST-1b, c, and d

We provide here additional figures of transit spectra computed for TRAPPIST-1b, c, and d, for H₂O-rich atmospheres (see Fig. B.1) and CO₂-rich atmospheres (see Fig. B.2). Transit spectra include realistic error bars calculated based on JWST Cycle 1 planned observations with NIRSpec and NIRISS.

B.2 Emission spectra and thermal phase curves of TRAPPIST-1b, c, and d

We provide here additional figures of emission spectra and thermal phase curves, for each individual GCM simulation of TRAPPIST-1b, c, and d. Each plot shows (on the left) the emission spectra at all phase angles and (on the right) the associated thermal phase curve, integrated over all wavelengths.

Fig. B.1 Transit spectra computed for the three innermost planets of the TRAPPIST-1 system calculated from 3D GCM simulations, assuming a 10 bar H2O atmosphere (+1 bar N2). The red and blue spectra are calculated for the wavelength coverage of NIRISS SOSS and NIRSpec PRISM, respectively. The 1σ error bars are calculated using the JWST Cycle 1 observation program (for more details, see Section 3.5.2). The transit spectra are in transit depth units (the transit depth is equal to zero at the surface).

Fig. B.2 Same as Fig. B.1, but for 10 bar CO2 atmospheres.

Fig. B.3 Thermal emission spectra at various phase angles (left). Thermal phase curve (right) integrated over all wavelengths.

Fig. B.4 Thermal emission spectra at various phase angles (left). Thermal phase curve (right) integrated over all wavelengths.

Fig. B.5 Thermal emission spectra at various phase angles (left). Thermal phase curve (right) integrated over all wavelengths.

Fig. B.6 Thermal emission spectra at various phase angles (left). Thermal phase curve (right) integrated over all wavelengths.

Fig. B.7 Thermal emission spectra at various phase angles (left). Thermal phase curve (right) integrated over all wavelengths.

Fig. B.8 Thermal emission spectra at various phase angles (left). Thermal phase curve (right) integrated over all wavelengths.

Fig. B.9 Thermal emission spectra at various phase angles (left). Thermal phase curve (right) integrated over all wavelengths.

Appendix C The effect of atmospheric expansion

Post-runaway H₂O-dominated atmospheres have been shown to be significantly expanded (Turbet et al. 2019, 2020b). Goldblatt (2015) showed that in extended atmospheres such as H₂O-rich atmospheres, the shortwave fluxes (i.e., emitted by the star) and longwave fluxes (i.e., emitted by the planet) can be significantly modified compared to the usual approximation. The absorbed and emitted fluxes are usually calculated by taking the observed radius of the planet as the reference cross section of the absorbed and emitted fluxes. To evaluate the error made by this approximation, we used 1D inverse radiative-convective calculations of H₂O-rich atmospheres (Turbet et al. 2019, 2020b) and calculated, for each spectral band in the radiative transfer, the altitude, and therefore the radius at which the opacity is equal to 1. We calculated a shortwave radius R_SW (and a longwave radius R_LW) by weighting the spectral radius by the stellar spectrum (and the OTR spectrum, respectively). Figure C.1 shows the scaling factor ratio as a function of surface temperature of the H₂O-rich atmosphere. The scaling factor ratio is defined as ${\left(\frac{R_{\text{LW}}}{R_{\text{SW}}}\right)}^2$ and provides an estimate of the error made on the TOA radiative budget. There are three important points to note. First, the scaling factor ratio is always close to (and within a few percent of) 1. This indicates that this expansion affects the water condensation limit to a few percent as well. The scaling factor ratio being almost always ≥ 1, this indicates that the ISR of the water condensation limit was slightly underestimated in our calculations. Second, the evolution of the scaling factor as a function of surface temperature is more sophisticated here than was calculated in Goldblatt (2015). This is because we calculated the τ =1 radius for each spectral band of the model, while Goldblatt (2015) calculated the radius of a wavelength-integrated τ = 1. As a result, and unlike what Goldblatt (2015) calculated, the scaling factor ratio evolves non-monotonically as a function of surface temperature. At low surface temperatures, the scaling factor ratio increases because R_LW increases much faster (the emission comes from a given atmospheric pressure, as shown in Boukrouche et al. 2021, so increasing the total atmospheric pressure increases R_LW) than R_SW (the stellar flux always reaches the surface). At intermediate surface temperatures, the scaling factor ratio reaches a plateau, because most of the shortwave flux (and longwave flux) are directly absorbed (and emitted, respectively) by the atmosphere. Increasing the total surface pressure shifts equally R_LW) and R_SW). At high surface temperatures, the scaling factor ratio decreases because the surface is hot enough that it starts to emit in the visible atmospheric windows of H₂O. Third, the scaling factor varies significantly depending on the total water vapor pressure assumed, as well as the type of host star. This is well illustrated in Fig. C.1. We note however that the behavior of the scaling factor ratio as a function of surface temperature is qualitatively similar for all the cases we explored.

It is important to note however that Fig. C.1 likely overestimates the effect of the atmospheric expansion on the TOA radiative budget. In this calculation we used vertical thermal profiles from 1D inverse radiative-convective models, which overestimate temperatures in the low atmosphere, atmospheric scale height, and thus atmospheric expansion (Selsis et al. 2023). When high water content is reached, the temperature profile in the lower atmosphere reaches an isotherm (see, e.g., Fig. 13c), thus limiting the atmospheric expansion and the effect on the TOA radiative budget.

Overall, the effect of atmospheric expansion is expected to play a minor role on our calculations of the water condensation zone.

Fig. C.1 Scaling factor ratio ${\left(\frac{R_{\text{LW}}}{R_{\text{SW}}}\right)}^2$ as a function of surface temperature, with RLW the mean planetary radius at which the OTR is emitted and RSW the mean planetary radius at which the ISR is absorbed. Each panel corresponds to a given water content (from 1 to 330 bar of H2O). The five colors correspond to the five different host stars.

References

All Tables

All Figures

Fig. 1 Water cloud column (top), albedo (middle), and thermal emission to space (bottom) horizontal maps, for an initially hot and steamy 2 R⊕ planet (with a 10 bar H2O + 1 bar N2 atmosphere) synchronously rotating (fixed rotation period of 60 Earth days) around an M3 star. The simulations shown here were forced to ISR (defined here and throughout the manuscript as the global mean average value, i.e., a factor of 4 lower than the stellar flux at the substellar point) ranging from 500 W m−2 or 1.47 F⊕ (right) down to 300 W m−2 or 0.88 F⊕ (left), and the outputs were averaged over 250 Earth days. The associated simulations names are M3-YANG-1, 2, 3, 4, and 5. For reference, insolation is 340.5 W m−2 on present-day Earth.

In the text

Fig. 2 Water cloud horizontal and vertical distributions. Water cloud Col. (a) in latitude vs. longitude coordinates and (b) tidally locked (TL) coordinates (b). Water cloud content vertical distribution (c) as a function of TL latitude (averaged over all TL longitudes). The plots were computed for an initially hot and steamy 2 R⊕ planet (with a 10 bar H2O + 1 bar N2 atmosphere) synchronously rotating (fixed rotation period of 60 Earth days) around an M3 star, and forced to an ISR of 500 W m−2 (simulation M3-YANG-1). The outputs were averaged over 250 Earth days. For reference, insolation is 340.5 W m−2 on present-day Earth. Wind fields are superimposed on the plots. Winds are plotted at a pressure of ~50 Pa in panels a and b, and averaged over all TL longitudes in panel c.

In the text

Fig. 3 Hysteresis loops and conditions of ocean formation for an initially hot and steamy 2 R⊕ planet (with a 10 bar H2O + 1 bar N2 atmosphere) synchronously rotating (fixed rotation period of 60 Earth days) around an M3 star (simulation M3-YANG-1). The figure depicts the evolution of surface temperature (top left), stratospheric water mixing ratio (top right), bond albedo (bottom left), and greenhouse effect (bottom right) as a function of the incoming stellar flux. The blue branches correspond to simulation results (Yang et al. 2013) that assume water is initially condensed on the surface. The red branches (this study) correspond to initially hot and steamy simulation results that assume water is initially in vapor form in the atmosphere.

In the text

Fig. 4 Water cloud column (top), albedo (middle), and thermal emission to space (bottom) horizontal maps, for an initially hot and steamy 1 R⊕ planet (with a 10 bar H2O + 1 bar N2 atmosphere) orbiting around different types of stars: TRAPPIST-1 (T★ ~ 2600 K), Proxima Centauri (T★ ~ 3050 K), an M3 star (T★ ~ 3400 K), a K5 star (T★ ~ 4400 K), and the Sun (T★ ~ 5780 K). The associated simulations names are T1-1, Pcen-1, M3-1, K5-1, and SUN-1, respectively. The simulations shown here were forced to an ISR of 500 W m−2 (i.e., 1.47 F⊕), and the outputs were averaged over 100 Earth days. A synchronous rotation for all the simulations was assumed, except for those of the Sun. In this case the results were taken directly from Turbet et al. (2021): rotation period of 1 Earth day, orbital period of 365 Earth days, and obliquity of 23.5°. For this simulation only, the maps were calculated in the heliocentric frame (i.e., keeping the subsolar point at 0° longitude and 0° latitude), and using an average of five Earth years as in Turbet et al. (2021).

In the text

Fig. 5 Conditions of ocean formation for initially hot and steamy 1 R⊕ planets (with a 10 bar H2O + 1 bar N2 atmosphere) orbiting around different types of stars: TRAPPIST-1 (T★ ~ 2600 K), Proxima Centauri (T★ ~ 3050 K), an M3 star (T★ ~ 3400 K), a K5 star (T★ ~ 4400 K) and the Sun (T★ ~ 5780 K). The associated simulations names are T1-1 to 8, Pcen-1 to 7, M3-1 to 7, K5-1 to 6 and SUN-1 to 6, respectively. The figure depicts the evolution of surface temperature (top) and bond albedo (bottom) as a function of the incoming stellar flux. We made the assumption of a synchronous rotation for all the simulations, except for those of the Sun (where we took the results directly from Turbet et al. (2021) in which rotation period is equal to 1 Earth day, orbital period is equal to 365 Earth days, and obliquity is equal to 23.5°). Horizontal dashed lines indicate the bond albedos of planets in a runaway greenhouse effect, obtained from 1D cloud-free calculations (Kopparapu et al. 2013). The dashed arrows indicate the ISR threshold limit at which the radiative TOA budget becomes imbalanced, leading to a runaway water condensation.

In the text

Fig. 6 Vertical profiles (global mean, envelope defined by global maximum and minimum values) of atmospheric temperatures (left), water cloud mixing ratio (middle), and absorbed stellar fluxes (right). The profiles were computed for an initially hot and steamy 1 R⊕ planet (with a 10 bar H2O + 1 bar N2 atmosphere) orbiting around different types of stars: TRAPPIST-1 (T★ ~ 2600 K), Proxima Centauri (T★ ~ 3050 K), an M3 star (T★ ~ 3400 K), a K5 star (T★ ~ 4400 K), and the Sun (T★ ~ 5780 K). The curves were drawn from the same GCM simulations as in Fig. 4. The associated simulations names are T1-1, Pcen-1, M3-1, K5-1, and SUN-1, respectively.

In the text

Fig. 7 Conditions of ocean formation for initially hot and steamy 1 R⊕ planets synchronously rotating around an M3 star (T★ ~ 3400 K). The figure depicts the evolution of surface temperature (top) and bond albedo (bottom) as a function of the incoming stellar flux, for four different bulk atmospheric compositions: 10 bar H2O + 1 bar N2 (M3), 3 bar H2O + 1 bar N2 (M3-3bar), 1 bar H2O + 1 bar N2 (M3-1bar), and 0.3 bar H2O + 1 bar N2 (M3-0.3bar).

In the text

Fig. 8 Various inner edge boundaries of the habitable zone (HZ) for different types of stars, based on the results of 1D calculations (Kopparapu et al. 2013) and 3D calculations for slow (Kopparapu et al. 2017; Way & Del Genio 2020) and fast rotators (Leconte et al. 2013a; Wolf & Toon 2015). The black solid and dashed lines were calculated based on 1D climate calculations, for the runaway greenhouse effect and the moist greenhouse effect, respectively (Kopparapu et al. 2013). The gray solid and dashed lines were calculated based on 3D climate calculations of tidally locked planets, for the runaway greenhouse effect and the moist greenhouse effect, respectively (Kopparapu et al. 2017). The red bold line (and red brackets) indicates the newly calculated water condensation limit. All the currently known planets and exoplanets (with a mass ≤ 5 M⊕ and/or a radius ≤ 1.6 R⊕) in or near the inner boundary of the HZ are shown (black squares). Also shown (gray brackets) are the incoming stellar radiation (ISR) they received at the ZAMS based on their age estimates, using the grid of stellar models of Baraffe et al. (2015). The water condensation limit (red bold line) should be compared with the ISR that planets received at the ZAMS (gray brackets).

In the text

Fig. 9 Cloud forcings for initially hot and steamy 1 R⊕ planets (with a 10 bar H2O + 1 bar N2 atmosphere) orbiting around different types of stars: TRAPPIST-1 (T★ ~ 2600 K), Proxima Centauri (T★ ~ 3050 K), an M3 star (T★ ~ 3400 K), a K5 star (T★ ~ 4400 K), and the Sun (T★ ~ 5780 K). The associated simulations names are T1-1 to 8, Pcen-1 to 7, M3-1 to 7, K5-1 to 6 and SUN-1 to 6, respectively. The figure depicts the greenhouse effect (top) of clouds, the amount of incident radiation (middle) reflected back by the clouds (the more negative the value, the greater the reflected flux), and the net radiative effect (bottom) of clouds (positive values mean warming), as a function of the incoming stellar flux. In all these simulations, clouds produce a strong atmospheric warming.

In the text

Fig. 10 Various positions of the water condensation limit, depending on the age of the planetary system considered. The standard water condensation limit (black dashed line) is valid at the ZAMS. It can be used as is to evaluate the condensation of water on a planet, provided that the insolation of the planet is taken at the ZAMS (see, e.g., gray brackets in Fig. 8). The other water condensation limits (solid colored lines) indicate, for a planet of age X (here with X = 0.5, 1, 4.5, and 10 Gyr), that if that planet is to the left of the age X curve, then it has never been able, or never will be able, to condense its water at the surface. These curves were calculated from the standard water condensation curve (black dashed line), and shifted using the grid of stellar models of Baraffe et al. (2015). Numerical values of all the curves are provided in Table 2.

In the text

Fig. 11 Water cloud column (top), albedo (middle), and thermal emission to space (bottom) horizontal maps, for TRAPPIST-1b (left), TRAPPIST-1c (middle), and TRAPPIST-1d (right). The associated simulations names are T1b, T1c, and T1d, respectively. The simulations assume initially hot and steamy tidally locked planets (with a 10 bar H2O + 1 bar N2 atmosphere), and the outputs were averaged over 100 Earth days. The dashed white lines indicate the location of the terminator.

In the text

Fig. 12 Temperature vertical profiles for TRAPPIST-1b (left), TRAPPIST-1c (middle), and TRAPPIST-1d (right). The associated simulation names are T1b, T1c, and T1d, respectively. The solid lines are global mean, temporal average (over 100 Earth days) temperature profiles calculated from GCM simulations. The shaded lines are snapshots taken at different latitude and longitude points to show the range of profiles within a simulation. The dashed lines are calculated with a 1D inverse radiative-convective model (Turbet et al. 2019, 2020b) similar to those used in Kopparapu et al. (2013) to calculate the inner edge of the HZ.

In the text

Fig. 13 Temperature vertical profiles (a) comparing TRAPPIST-1b, c, and d (simulations T1b, T1c, and T1d), (b) several internal heat fluxes for TRAPPIST-1b (simulations T1b, T1b-Fgeo1, T1b-Fgeo5, and T1b-Fgeo25), and (c) several water partial pressures for TRAPPIST-1b (simulations T1b, T1b-3bar, T1b-1bar). The solid lines are global mean, temporal average (over 100 Earth days) temperature profiles calculated from GCM simulations. The dashed lines are calculated with a 1D inverse radiative-convective model (Turbet et al. 2019, 2020b) similar to those used in Kopparapu et al. (2013) to calculate the inner edge of the HZ.

In the text

Fig. 14 Various transit spectra computed with PSG for the three innermost planets of the TRAPPIST-1 system calculated from 3D GCM simulations (assuming H2O-dominated composition, and also CO2-dominated for comparison). In the upper panel a, transit spectra of TRAPPIST-1b for various compositions and pressures are compared; in the lower panel b, transit spectra of TRAPPIST-1b, c, and d for two compositions only (10 bar H2O + 1 bar N2; 10 bar CO2) are compared. The dashed lines indicate H2O cloud-free atmospheres, while solid lines assume H2O cloudy atmospheres. The dotted lines indicate CO2 atmospheres (including the effect of CO2 ice clouds, although extremely limited here; instead, CO2 condensation was cut off for TRAPPIST-1b and c simulations with 100 mbar of CO2 as this would cause the atmosphere to collapse on the nightside.). These are relative transit depths; the curves are shifted with arbitrary offsets to compare them more easily.

In the text

Fig. 15 Emission spectra of TRAPPIST-1b (upper panel) and c (lower panel), calculated at the secondary eclipse. The spectra are calculated using GCM simulations of H2O-rich (blue to gray) and CO2-rich atmospheres (red to orange). Emission spectra are also computed for airless planets (in bold black). The 1σ error bars are calculated using the JWST Cycle 1 observation program (for more details, see Sect. 3.5.3).

In the text

	Fig. 16 Same as Fig. 15, but with the emission converted into brightness temperature units.
In the text

	Fig. 17 Thermal phase curves of TRAPPIST-1b (upper panel) and c (lower panel), integrated in the F1500W filter (centered around 15 μm), for various atmospheric scenarios: H2O-rich (blue to gray), CO2-rich atmospheres (red to orange), airless planet (bold black). The H2O-rich atmospheres also include 1 bar of N2.
In the text

Fig. B.1 Transit spectra computed for the three innermost planets of the TRAPPIST-1 system calculated from 3D GCM simulations, assuming a 10 bar H2O atmosphere (+1 bar N2). The red and blue spectra are calculated for the wavelength coverage of NIRISS SOSS and NIRSpec PRISM, respectively. The 1σ error bars are calculated using the JWST Cycle 1 observation program (for more details, see Section 3.5.2). The transit spectra are in transit depth units (the transit depth is equal to zero at the surface).

In the text

	Fig. B.2 Same as Fig. B.1, but for 10 bar CO2 atmospheres.
In the text

	Fig. B.3 Thermal emission spectra at various phase angles (left). Thermal phase curve (right) integrated over all wavelengths.
In the text

	Fig. B.4 Thermal emission spectra at various phase angles (left). Thermal phase curve (right) integrated over all wavelengths.
In the text

	Fig. B.5 Thermal emission spectra at various phase angles (left). Thermal phase curve (right) integrated over all wavelengths.
In the text

	Fig. B.6 Thermal emission spectra at various phase angles (left). Thermal phase curve (right) integrated over all wavelengths.
In the text

	Fig. B.7 Thermal emission spectra at various phase angles (left). Thermal phase curve (right) integrated over all wavelengths.
In the text

	Fig. B.8 Thermal emission spectra at various phase angles (left). Thermal phase curve (right) integrated over all wavelengths.
In the text

	Fig. B.9 Thermal emission spectra at various phase angles (left). Thermal phase curve (right) integrated over all wavelengths.
In the text

	Fig. C.1 Scaling factor ratio ${\left(\frac{R_{\text{LW}}}{R_{\text{SW}}}\right)}^2$ as a function of surface temperature, with RLW the mean planetary radius at which the OTR is emitted and RSW the mean planetary radius at which the ISR is absorbed. Each panel corresponds to a given water content (from 1 to 330 bar of H2O). The five colors correspond to the five different host stars.
In the text