© ESO, 2013

1. Beyond the runaway greenhouse limit

Because of the bias of planet characterization observatories toward strongly irradiated objects on short orbits, planets identified as near the inner edge of the habitable zone represent valuable targets. This does, however, raise numerous questions concerning our understanding of the processes delimiting this inner edge.

For a planet like the Earth, which from the point of view of climate can be considered as an aqua planet¹, a key constraint is provided by the positive feedback of water vapor on climate. When the insolation of the planet is increased, more water vapor is released in the atmosphere by the hotter surface, increasing the greenhouse effect. If on Earth today this effect is balanced by other processes, in particular the increase in infrared emission by the surface with temperature, it has been shown that above a critical absorbed stellar flux, the planet cannot reach radiative equilibrium balance and is heated until the surface can radiate at optical wavelength around a temperature of 1400 K (Kasting et al. 1984; Kasting 1988). All the water at the surface is vaporized. The planet is in a runaway greenhouse state. Roughly speaking, this radiation limit stems from the fact that the surface must increase its temperature to radiate more, but the amount of water vapor in the troposphere, hence the opacity of the latter, increases exponentially with this surface temperature (because of the Clausius-Clapeyron relation). Then the thermal emission of the opaque atmosphere is emitted at the level at which the optical depth is close to unity, a level whose temperature reaches an asymptotic value, hence the limited flux² (Nakajima et al. 1992; Pierrehumbert 2010).

While large uncertainties remain regarding the radiative effect of water vapor continuum and the nonlinear, 3D effect of dynamics and clouds, the critical value of the average absorbed stellar radiation, (where F_⋆ is the stellar flux at the substellar point and is the planetary bond albedo), should conservatively range between 240 W/m² (the mean absorbed stellar radiation on Earth) and 350 W/m² (Abe et al. 2011). Terrestrial planets emitting such (low) IR fluxes will be difficult to characterize in the near future.

Fortunately, the runaway greenhouse limit discussed above is not as fundamental as it may appear. Indeed, the theory outlined above assumes that a sufficient³ reservoir of liquid water is available everywhere on the surface (Selsis et al. 2007). While it would seem from 1D modeling that liquid water would disappear at a lower surface temperature for a planet with a limited water inventory, it goes the opposite way in 3D. Global climate models show that such “land planets” can have a wider habitable zone (Abe et al. 2011). Near the inner edge, this property stems from the fact that atmospheric transport of water from warm to cold regions of the planets is not counteracted by large-scale surface runoff (Abe et al. 2005). Heavily irradiated regions are then drier and can emit more IR flux than the limit found for a saturated atmosphere.

In principle, it is thus possible to find planets that absorb and re-emit more flux than the runaway greenhouse threshold and for which liquid or solid water can be thermodynamically stable at the surface, provided that they have accreted limited water supplies (or lost most of them). As shown hereafter, this could be the case for Gl 581 c (Udry et al. 2007; Mayor et al. 2009; Forveille et al. 2011) and HD 85512 b (Pepe et al. 2011), the only two confirmed planets having masses that are compatible with a terrestrial planet (M_p ≲ 7−8  M_⊕; although defining such a mass limit is questionable) and that lie beyond the inner edge of the classical habitable zone but receive a moderate stellar flux, F_⋆/4 < 1000 W/m². These planets lie even closer to the star that the inner edge of the limits set by Abe et al. (2011) for land planets.

However, for those close-in planets whose rotation states have been strongly affected by tidal dissipation and which are either in a state of synchronous, pseudo-synchronous⁴, or resonant rotation with near zero obliquity (Goldreich & Peale 1966; Heller et al. 2011; Makarov et al. 2012), efficient and permanent cold traps present on the dark hemisphere or at the poles could irreversibly capture all the available water in a permanent ice cap. To tackle this issue, we performed several sets of global climate model simulations for prototype land planets with the characteristic masses and (strong) incoming stellar fluxes of Gl 581 c and HD 85512 b. Three-dimensional climate models are needed to assess the habitability of land planets in such scenarios because horizontal inhomogeneities in impinging flux and water distribution are the key features governing the climate. This may explain why this scenario has received little attention so far.

One of the difficulties lies in the lack, as will be extensively discussed here, of a precise flux limit above which a runaway greenhouse state is inevitable. Instead, for a wide range of parameters, two stable equilibrium climate regimes exist because of the very strong positive feedback of water vapor. The actual regime is determined not only by the incoming stellar flux and spectrum, or by the planet’s atmosphere, but also by the amount of water available and its initial distribution. Consequently, instead of trying to cover the whole parameter space, here we have performed simulations for two already detected planets. We show that qualitatively robust features appear in our simulations and argue that the processes governing the climate discussed here are quite general.

In Sect. 2 we introduce the various components of our numerical climate model and present the basic climatology of our model in a completely dry case (Sect. 3). We show that the atmospheric circulation can settle into two different regimes – namely, super-rotation or stellar/anti-stellar circulation – and that the transition of one regime to the next occurs when the equatorial Rossby deformation radius significantly exceeds the planetary radius. We also show that strong discrepancies can occur between 1D and globally averaged 3D results. Then we show that when water is present, the climate regime reached depends on the initial conditions and can exhibit a bistability (Sect. 4). We argue that this bistability is a dynamical manifestation of the interplay between the runaway greenhouse instability on the hot hemisphere and the cold trapping arising either on the dark side or near the poles. In Sect. 5, we discuss the long-term stability of the different regimes found here and show that, if atmospheric escape is efficient, it could provide a stabilizing feedback for the climate. In Sect. 6, we suggest that if the planet has been able to accumulate a large amount of ice in its cold traps such that ice draining by glaciers is significant, liquid water could be present at the edge of the ice cap where glaciers melt. We also discuss the possible presence of subsurface water under a thick ice cap. Finally, in Sect. 7 we present synthetic spectra and discuss the observable signatures of the two climate regimes discussed above.

2. Method

Our simulations were performed with an upgraded version of the LMD generic global climate model (GCM) specifically developed for the study of extrasolar planets (Wordsworth et al. 2010, 2011; Selsis et al. 2011) and paleoclimates (Wordsworth et al. 2013; Forget et al. 2013). The model uses the 3D dynamical core of the LMDZ 3 GCM (Hourdin et al. 2006), based on a finite-difference formulation of the primitive equations of geophysical fluid dynamics. A spatial resolution of 64  ×  48 × 20 in longitude, latitude, and altitude is used for the simulations.

2.1. Radiative transfer

The method used to produce our radiative transfer model is similar to Wordsworth et al. (2011). For a given mixture of atmospheric gases, here N₂ with 376 ppmv of CO₂ and a variable amount of water vapor, we computed high-resolution spectra over a range of temperatures and pressures using the HITRAN 2008 database (Rothman et al. 2009). Because the CO₂ concentration depends on the efficiency of complex weathering processes, it is mostly unconstrained for hot synchronous land planets. The choice of an Earth-like concentration is thus arbitrary, but our prototype atmosphere should be representative of a background atmosphere with a low greenhouse effect. For this study we used temperature and pressure grids with values T =  {110,170,...,710}  K,p =  {10^-3,10^-2,...,10⁵}  mbar. The H₂O volume mixing ratio could vary in the range {10^-7,10^-6,...,1}. The H₂O lines were truncated at 25 cm^-1, while the water vapor continuum was included using the CKD model (Clough et al. 1989). We also account for opacity due to N₂−N₂ collision-induced absorption (Borysow & Frommhold 1986; Richard et al. 2011).

The correlated-k method was then used to produce a smaller database of coefficients suitable for fast calculation in a GCM. Thanks to the linearity of the Schwarzschild equation of radiative transfer, the contribution of the thermal emission and downwelling stellar radiation can be treated separately, even in the same spectral channel. We therefore do not assume any spectral separation between the stellar and planetary emission wavelengths. For thermal emission, the model uses 19 spectral bands, and the two-stream equations are solved using the hemispheric mean approximation (Toon et al. 1989). Infrared scattering by clouds is also taken into account. Absorption and scattering of the downwelling stellar radiation is treated with the δ-Eddington approximation within 18 bands. Sixteen points are used for the -space integration, where is the cumulated distribution function of the absorption data for each band.

As in Wordsworth et al. (2010), the emission spectrum of our M dwarf, Gl 581, was modeled using the Virtual Planet Laboratory AD Leo data (Segura et al. 2005). For the K star HD 85512, we used the synthetic spectrum of a K dwarf with an effective temperature of 4700 K, g = 10^2.5 m s^-2 and [M/H] =  −0.5 dex from the database of Allard et al. (2000).

2.2. Water cycle

In the atmosphere, we follow the evolution of water in its vapor and condensed phases. These tracers are advected by the dynamical core, mixed by turbulence and dry and moist convection. Much care has been devoted to develop a robust and numerically stable water cycle scheme that is accurate both in the trace gas (water vapor mass mixing ratio q_v ≪ 1) and dominant gas (q_v ~ 1) limit. In particular, the atmospheric mass and surface pressure variation (and the related vertical transport of tracers through pressure levels) due to any evaporation/sublimation/condensation of water vapor is taken into account with a scheme similar to the one developed by Forget et al. (2006) for CO₂ atmospheres.

Cloud formation is treated using the prognostic equations of Le Treut & Li (1991). For each column and level, this scheme provides the cloud fraction f_c and the mass mixing ratio of condensed water q_l, which are both functions of q_v and the saturation vapor mass mixing ratio q_s. In addition, when part of a column reaches both 100% saturation and a super-moist-adiabatic lapse rate, moist convective adjustment is performed following Manabe & Wetherald (1967), and the cloud fraction is set to unity. In the cloudy part of the sky, the condensed water is assumed to form a number N_c of cloud particles per unit mass of moist air. This quantity, which represents the number density per unit mass of activated cloud condensation nuclei, is assumed to be uniform throughout the atmosphere (Forget et al. 2013). Then, the effective radius of the cloud particles is given by r_eff = (^3 q_l/_{4 π ρc Nc})^1/3, where ρ_c is the density of condensed water (10³ kg/m³ for liquid and 920 kg/m³ for ice). Precipitation is computed with the scheme of Boucher et al. (1995). Since this scheme explicitly considers the dependence on gravity, cloud particle radii, and background air density, it should remain valid for a wide range of situations. Finally, the total cloud fraction of the column is assumed to be the cloud fraction of the level with the thickest cloud, and radiative transfer is computed both in the cloudy and clear sky regions. Fluxes are then linearly weighted between the two regions.

On the surface, ice can also have a radiative effect by linearly increasing the albedo of the ground to A_ice (see Table 1; in the present version of the code, both surface and ice albedos are independent of the wavelength) until the ice surface density exceeds a certain threshold (here 30 kg m^-2). In the general case, the values taken for both of these parameters can have a significant impact on the mean radiative balance of a planet. In addition, providing an accurate ice mean albedo can prove difficult, because it depends on the ice cover, light incidence, state of the ice/snow, but also on the stellar spectrum (ice absorbing more efficiently infrared light; Joshi & Haberle 2012). However, in our case of a strongly irradiated, spin synchronized land planet where water is in limited supply, most of the day side is above 0°C, and ice can form only in regions receiving little or no visible flux. As a result, and as confirmed by our tests, the value chosen for A_ice has a negligible impact on the results presented hereafter.

3. Climate on synchronous planets: what 1D cannot tell you

So far, most studies concerned with the determination of the limits of the habitable zone have used 1D calculations. While these single-column models can provide reasonable answers for planets with (very) dense atmospheres and/or a rapid rotation that limits large-scale temperature contrasts, we show here that only a 3D model can treat the case of synchronized or slowly rotating exoplanets properly. To that purpose we quantify the discrepancies that can arise between 1D and 3D models both analytically and numerically.

3.1. Equilibrium temperatures vs physical temperatures

When considering the energetic balance of an object, one of the first quantities coming to mind is the equilibrium temperature, i.e. the temperature that a blackbody would need to have to radiate the absorbed flux. If one assumes local thermal equilibrium for each atmospheric column, one can define the local equilibrium temperature $\begin{array}{ccc}{\mathit{T}}_{\mathrm{eq}\mathit{,}{\mathit{\mu }}_{\mathit{\star }}}\mathrm{\equiv }{\left(\frac{\mathrm{(}\mathrm{1}\mathrm{-}\mathit{A}\mathrm{)}\hspace{0.17em}{\mathit{F}}_{\mathit{\star }}\hspace{0.17em}{\mathit{\mu }}_{\mathit{\star }}}{\mathit{\sigma }}\right)}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathit{,}& & \end{array}$(1)where F_⋆ is the stellar flux at the substellar point at the time considered, A the local bond albedo, μ_⋆ the cosine of the stellar zenith angle (assumed equal to zero when the star is below the horizon), and σ the Stefan-Boltzmann constant.

Relaxing the local equilibrium condition but still assuming global thermal equilibrium of the planet, one can define the effective equilibrium temperature $\begin{array}{ccc}{\mathit{T}}_{\mathrm{eq}}\mathrm{\equiv }{\left(\frac{\mathrm{1}}{\mathit{S}\hspace{0.17em}\mathit{\sigma }}\mathrm{\int }\mathrm{(}\mathrm{1}\mathrm{-}\mathit{A}\mathrm{)}\hspace{0.17em}{\mathit{F}}_{\mathit{\star }}\hspace{0.17em}{\mathit{\mu }}_{\mathit{\star }}\mathrm{d}\mathit{S}\right)}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathrm{\equiv }{\left(\frac{\mathrm{(}\mathrm{1}\mathrm{-}\mathit{A̅}\mathrm{)}\hspace{0.17em}{\mathit{F}}_{\mathit{\star }}}{\mathrm{4}\hspace{0.17em}\mathit{\sigma }}\right)}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathit{,}& & \end{array}$(2)where S is the planet’s surface area and the effective, planetary bond albedo.

It is important to note that, because T_eq is derived from global energetic considerations, it is not necessarily close to the mean surface temperature, , and is often not. This is for two main reasons:

• The most widely recognized of these reasons is that, when an absorbing atmosphere is present, the emitted flux does not directly come from the ground. Because of this greenhouse effect, the surface temperature of an isolated atmospheric column can be higher than its equilibrium temperature. This explains why the mean surface temperatures of Venus and Earth are 737 K and 288 K, while their effective equilibrium temperatures are around 231 K and 255 K (assuming an albedo of 0.75 and 0.29).

• The second reason is the nonlinearity of the blackbody bolometric emission with temperature. Because this emissions scales as T⁴, it can be demonstrated that the spatial average of the local equilibrium temperature (Eq. (1)) is necessarily lower than the effective equilibrium temperature defined by Eq. (2)⁵. Rigorously, this stems from the fact that T_eq,μ⋆ is a concave function of the absorbed flux that, with the help of the Rogers-Hölder inequality (Rogers 1888; Hölder 1889), yields $\begin{array}{ccc}\frac{\mathrm{1}}{\mathit{S}}\mathrm{\int }{\mathit{T}}_{\mathrm{eq}\mathit{,}{\mathit{\mu }}_{\mathit{\star }}}\mathrm{d}\mathit{S}⩽{\left(\frac{\mathrm{1}}{\mathit{S}}\mathrm{\int }{\mathit{T}}_{\mathrm{eq}\mathit{,}{\mathit{\mu }}_{\mathit{\star }}}^{\mathrm{4}}\mathrm{d}\mathit{S}\right)}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathrm{\equiv }{\mathit{T}}_{\mathrm{eq}}\mathit{.}& & \end{array}$(3)In other words, regions receiving twice as much flux do not need to be twice as hot to reach equilibrium, and cold regions have more weight in the mean temperature.

To give an idea of how strong this second effect can be, we consider the limiting case of an airless, synchronous planet in thermal equilibrium with a constant surface albedo, . In this case, the local surface temperature is equal to the local equilibrium temperature, and the mean surface temperature is given by $\begin{array}{ccc}\mathit{T̅}\mathrm{s}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathrm{d}\mathit{\phi }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}\mathrm{d}{\mathit{\mu }}_{\mathit{\star }}{\left(\frac{\mathrm{(}\mathrm{1}\mathrm{-}\mathit{A̅}\mathrm{)}\hspace{0.17em}{\mathit{F}}_{\mathit{\star }}\hspace{0.17em}{\mathit{\mu }}_{\mathit{\star }}}{\mathit{\sigma }}\right)}^{\mathrm{1}\mathit{/}\mathrm{4}}\\ & \mathrm{=}& \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}{\left(\frac{\mathrm{(}\mathrm{1}\mathrm{-}\mathit{A̅}\mathrm{)}\hspace{0.17em}{\mathit{F}}_{\mathit{\star }}}{\mathit{\sigma }}\right)}^{\mathrm{1}\mathit{/}\mathrm{4}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathrm{d}\mathit{\phi }\hspace{0.17em}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}{\mathit{\mu }}_{\mathit{\star }}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathrm{d}{\mathit{\mu }}_{\mathit{\star }}\mathit{.}\end{array}$(4)Performing the integrals, and using Eq. (2), we find $\begin{array}{ccc}\mathit{T̅}\mathrm{s}\mathrm{=}\frac{\mathrm{2}\sqrt{\mathrm{2}}}{\mathrm{5}}{\mathit{T}}_{\mathrm{eq}}\mathrm{\approx }\mathrm{0.57}\hspace{0.17em}{\mathit{T}}_{\mathrm{eq}}\mathit{.}& & \end{array}$(5)While any energy redistribution by the atmosphere will tend to lower this difference⁶, one needs to remember that even the mean surface temperature can be much lower than the effective equilibrium temperature.

These simple considerations show that equilibrium temperature can be misleading and should be used with great care. Indeed, even if 1D models can take the greenhouse effect into account, they cannot predict energy redistribution on the night side. Consequently, not only can they not predict temperatures on the night side, but they also predict fairly inaccurate mean temperatures. Thus, while equilibrium temperature remains a handy tool for comparing the temperature regime from one planet to another, it is not accurate enough to draw conclusions on the climate of any given object.

3.2. Climate on a dry synchronous planet

We now present a first set of 3D simulations for Gl 581 c and HD 85512 b that are modeled as completely dry planets. The goal of these idealized simulations is to gain insight into the various circulation patterns arising on synchronously rotating terrestrial planets and into how the circulation regime depends on the atmospheric mass (mean surface pressure) and planetary rotation rate.

3.2.1. Surface temperatures

These dry runs were performed using the parameters from Table 1, without a water cycle. For various surface pressures, the model was initialized from a rest state with uniform surface temperature and integrated until statistical equilibrium was reached.

The temperature maps for three sample surface pressures for the two planets are shown in Fig. 1. The most apparent feature is the expected dichotomy between the sunlit (between −90° and 90° longitude) and dark hemispheres. A surprising outcome of these simulations is the large temperature difference between the substellar point and the night side compared to the one in synchronous Earth simulations of Joshi (2003) and Merlis & Schneider (2010) or to the equator-pole temperature difference on Earth. To understand these differences, we have performed an additional suite of simulations of fast rotating or tidally synchronized Earth with our full radiative transfer or with a gray gas approximation (as in Merlis & Schneider 2010) and with or without a water cycle. The model recovers both Earth present climatology and results similar to Merlis & Schneider (2010).

Thanks to this suite of simulations, the causes of this strong temperature contrast are found to be (more or less by order of importance):

• 1.

  Latent cooling/heating. When water is abundantly present, the temperature at the substellar point is controlled mostly thermodynamically and is much lower than in the dry case. Latent energy transported by the atmosphere to the night side also contributes significantly to the warming of the latter.

• 2.

  Non-gray opacities. In our case dry case, greenhouse is solely due to the 15 μm band of CO₂ (especially on the night side). Thus, using a gray opacity of the background atmosphere with an optical depth near unity, as done in Merlis & Schneider (2010), considerably underestimates the cooling of the surface and thus overestimates the night side temperature.

• 3.

  A different planet. When its bulk density is constant, both the size and the surface gravity increase with the planet’s size. Apart from the dynamical effect due to the variation in Rossby number (extensively discussed below), for a given surface pressure (p_s), both the mass of the atmosphere and its radiative timescale (τ_rad) decrease with gravity (and with incoming flux for the latter) and that the advective timescale (τ_adv) increases with the radius. For a given surface pressure, the energy transport, thought to scale with the atmospheric mass and τ_rad/τ_adv, thus tends to be less efficient on a more massive and more irradiated object.

Indeed, as expected, the temperature difference decreases with increasing total pressure, because the strength of the advection of sensible heat to the night side increases with atmospheric mass. Both because of this more efficient redistribution and of the greenhouse effect of the atmosphere, the globally averaged mean surface temperature increases with surface pressure, as shown in Fig. 2.

For comparison, we also computed the surface temperature yielded by the 1D version of our model. In this case, the atmospheric column is illuminated with a constant flux equal to the globally averaged impinging flux (1/4 of the substellar value) and with a stellar zenith angle of 60°. Apart from this prescription for the flux and the fact that large-scale dynamics is turned off, our 1D code is exactly the same as our 3D one. As can be seen in Fig. 2, even with a 10 bar atmosphere, the 1D model is still 30 K hotter than the globally averaged 3D temperature. This difference increases up to 70 K for the 100 mb case because redistribution is less efficient (see Sect. 3.1). A surface albedo of 0.3 is assumed here, similar to that in desert regions on Earth. The sensitivity of the surface temperature when albedo is varied between 0.1 and 0.5 is shown in Fig. 2.

A troubling and quite counterintuitive outcome of these simulations is that, as shown in Fig. 2, although Gl 581 c receives more flux and has a higher mean surface temperature than HD 85512 b, its minimum temperatures are always lower. As can be seen in Fig. 1, the models for Gl 581 c present very cold regions located at mid latitudes and west of the western terminator (also called morning limb), whereas the HD 85512 b models have a rather uniform (and higher in average) temperature over the whole night side.

3.2.2. Circulation regime

Further investigations have shown that this difference is a signature of the distinct circulation regimes arising on the two planets. Indeed, as shown in Fig. 3, our simulations of the atmosphere of Gl 581 c develop eastward jets (super-rotation) as found in many GCM of tidally locked exoplanets (Showman et al. 2008, 2011, 2013; Thrastarson & Cho 2010; Heng & Vogt 2011; Heng et al. 2011; Showman & Polvani 2011). Conversely, the atmosphere of HD 85512 b seems to settle predominantly into a regime of stellar/anti stellar circulation where high-altitude winds blow from the day side to night side in a symmetric pattern with respect to the substellar point (upper right panel of Fig. 3; similar to the slowly rotating case of Merlis & Schneider 2010).

To understand these distinct circulation regimes, one needs to understand the mechanism responsible for the onset of super-rotation on tidally locked exoplanets. Showman & Polvani (2011) show that the strong longitudinal variations in radiative heating causes the formation of standing, planetary-scale equatorial Rossby and Kelvin waves. By their interaction with the mean flow, these waves transport eastward momentum from the mid latitudes toward the equator, creating and maintaining an equatorial jet. As a result, they predict that the strength of the jet increases with the day/night temperature difference and that the latitudinal width of this jet scales as the equatorial Rossby deformation radius $\begin{array}{ccc}{\mathit{L}}_{\mathrm{Ro}}\mathrm{\equiv }\sqrt{\frac{\mathit{N}\hspace{0.17em}\mathit{H}}{\mathit{\beta }}}\mathit{,}& & \end{array}$(6)where N is the Brunt-Väisälä frequency, H the pressure scale height of the atmosphere (N H being the typical speed of gravity waves), and β ≡ ∂f/∂y the latitudinal derivative of the Coriolis parameter, f ≡ 2 Ω sinθ, with Ω the planet rotation rate and θ the latitude (β = 2 Ω/R_p at the equator). Roughly, the equatorial Rossby deformation radius is the typical length scale over which an internal (gravity) wave created at the equator is significantly affected by the Coriolis force. This allows us to define a dimensionless equatorial Rossby deformation length $\begin{array}{ccc}\mathrm{ℒ}\mathrm{\equiv }\frac{{\mathit{L}}_{\mathrm{Ro}}}{{\mathit{R}}_{\hspace{0.17em}\mathrm{p}}}\mathrm{=}\sqrt{\frac{\mathit{N}\hspace{0.17em}\mathit{H}}{\mathrm{2}\hspace{0.17em}\mathrm{\Omega }\hspace{0.17em}{\mathit{R}}_{\hspace{0.17em}\mathrm{p}}}}\mathrm{·}& & \end{array}$(7)This latter prediction is directly relevant to our case. Indeed, although they have rather similar sizes and gravity, Gl 581 c rotates around its star in ~13 Earth days (d) whereas HD 85512 b, which orbits a brighter K star, has a ~58 d orbital period. If spin orbit synchronization has occurred, the Coriolis force is thus expected to have a much weaker effect on the circulation of HD 85512 b. To quantify this, we evaluate ℒ in our case. In a dry stably stratified atmosphere with a lapse rate dT/dz, the Brunt-Väisälä frequency is given by $\begin{array}{ccc}{\mathit{N}}^{\mathrm{2}}\mathrm{=}\frac{\mathit{g}}{\mathit{T}}\left(\frac{\mathit{g}}{{\mathit{c}}_{\mathrm{p}}}\mathrm{+}\frac{\mathrm{d}\mathit{T}}{\mathrm{d}\hspace{0.17em}\mathit{z}}\right)\mathit{,}& & \end{array}$(8)c_p being the heat capacity of the gas at constant pressure. The pressure scale height is given by $\begin{array}{ccc}\mathit{H}\mathrm{\equiv }\frac{{\mathit{k}}_{\hspace{0.17em}\mathrm{B}}\mathit{T}}{{\mathit{m}}_{\mathrm{a}}\hspace{0.17em}\mathit{g}}\mathit{,}& & \end{array}$(9)where k_B is the Boltzmann constant and m_a the mean molecular weight (in kg) of the air. Assuming an isothermal atmosphere, this yields $\begin{array}{ccc}\mathrm{ℒ}\mathrm{=}\sqrt{\frac{{\mathit{k}}_{\hspace{0.17em}\mathrm{B}}}{{\mathit{m}}_{\mathrm{a}}\hspace{0.17em}{\mathit{c}}_{\mathrm{p}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\frac{{\mathit{T}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{\mathrm{2}\hspace{0.17em}\mathrm{\Omega }\hspace{0.17em}{\mathit{R}}_{\hspace{0.17em}\mathrm{p}}}}\mathrm{·}& & \end{array}$(10)Assuming an N₂-dominated atmosphere, using the equilibrium temperature as a proxy for T and disregarding the latitude factor, one finds ℒ = 1.1 for Gl 581 c and 2.5 for HD 85512 b. In our slowly rotating case, the planet is thus too small to contain a planetary Rossby wave, and the mechanism described above is too weak to maintain super-rotation. This also explains why the jet is planet-wide in the more rapidly rotating case of Gl 581 c. To confirm this analysis, we ran a simulation of Gl 581 c where we, everything else being equal, reduced the rotation rate by a factor of five so that the ℒ is the same as for HD 85512 b (Ω = Ω_orb/5 case; middle column of Fig. 3). As expected from the dimensional analysis, in this slowly rotating case, no jets appear, and we recover a stellar/antistellar circulation. Another argument in favor of this analysis is that the location and shape of the cold regions and wind vortices present in the Gl 581 c case (see Figs. 1 and 3) are strongly reminiscent of Rossby-wave gyres seen in the Matsuno-Gill standing wave pattern (Matsuno 1966; Gill 1980; Showman & Polvani 2011).

Another way to look at this is to consider a multiway force balance (Showman et al. 2013). One way to tilt eddy velocities northwest-southeast in the northern hemisphere and southwest-northeast in the southern hemisphere, driving super-rotation, is to have a force balance between pressure-gradient (∇_pΦ ≡ k_BΔTΔlnp/(m_aR_p), where ΔT is the horizontal day-night temperature difference and Δlnp is logarithmic pressure range over which this temperature difference is maintained), Coriolis (≈f u where u is the wind speed), and advection (≈u²/R_p) acceleration (see Showman et al. (2013) for a discussion on the effect of drag). In this case one should have $\begin{array}{ccc}{\mathrm{\nabla }}_{\mathrm{p}}\mathrm{\Phi }\mathrm{~}\mathit{f}\hspace{0.17em}\mathit{u}\mathrm{~}\frac{{\mathit{u}}^{\mathrm{2}}}{{\mathit{R}}_{\hspace{0.17em}\mathrm{p}}}\mathit{,}& & \end{array}$(11)which leads to u ~ f R_p and eventually $\mathit{f}\mathrm{~}\sqrt{{\mathrm{\nabla }}_{\mathrm{p}}\mathrm{\Phi }\mathit{/}{\mathit{R}}_{\hspace{0.17em}\mathrm{p}}}$.

If $\mathit{f}\mathrm{\ll }\sqrt{{\mathrm{\nabla }}_{\mathrm{p}}\mathrm{\Phi }\mathit{/}{\mathit{R}}_{\hspace{0.17em}\mathrm{p}}}$ this three-way force balance should reduce to a two-way force balance between advection and pressure-gradient accelerations, which in turn leads to symmetric day-night flow. A numerical estimate with the same parameters as above with ΔT = 100  K and Δlnp = 1 yields a critical rotation period of about 10 d, close to Gl 581 c rotation period, but much smaller than the one of HD 85512 b. This confirms that HD 85512 b should have a much more pronounced day-night circulation pattern than Gl 581 c, but it also tells us that the Coriolis force is already quite weak in the force balance of Gl 581 c, explaining the significant meridional wind speeds seen in Fig. 3.

There are a few general differences with previous close-in giant planet 3D simulations, however. First, in all our simulations, low-altitude winds converge toward the substellar point. This results from the strong depression caused there by a net mass flux horizontal divergence high in the atmosphere driving a large-scale upward motion. Second, even at levels where eastward winds are the strongest (in the simulations with super-rotation), the wind speeds exhibit a minimum (possibly negative) along the equator, westward from the substellar point. At a fixed altitude, this results from a pressure maximum occurring there because of the higher temperature, and consequently higher pressure scale height.

Another feature visible in Fig. 4 is the dependence of the strength of the zonal jet on the atmospheric mean surface pressure (lower pressure causing faster winds). A detailed analysis is beyond the scope of this study. Nevertheless, these simulations seem to confirm that zonal wind speeds increase with forcing amplitude (Showman & Polvani 2011). Indeed, an increasing pressure entails a decreasing thermal contrast between the two hemispheres and a smaller forcing, hence the lower wind speeds.

Surprisingly, the presence of a strong equatorial jet lowers the heat redistribution efficiency by the atmosphere. By focusing the flow around the equator, the Coriolis force prevents an efficient redistribution to the midlatitudes on the night side. To quantify this redistribution in a simple way, we computed the ratio η of the mean outgoing thermal flux that is emitted by the night side over the mean outgoing thermal flux from the day side (see Fig. 5). With this definition, no redistribution should yield η = 0, whereas a very homogeneous emission would produce η = 1. As expected, the redistribution increases with total pressure. At constant pressure, redistribution seems to be controlled by the equatorial Rossby number and is more important in the slowly rotating cases (HD 85512 b and Ω = Ω_orb/5 Gl 581 c case).

Further idealized simulations are needed to fully understand the mechanisms controlling both the circulation pattern and the energy redistribution on slowly rotating, synchronous terrestrial planets.

3.2.3. Thermal profile

Finally, in Fig. 6, we present typical temperature profiles obtained from our simulations (here in the 1 bar case, but profiles obtained for different pressures are qualitatively similar). For comparison, Fig. 6 also shows the profile obtained from the corresponding global 1D simulation, and the spatially averaged profile from our 3D run.

One can directly see that, in the lower troposphere, temperature decreases with altitude on the day side following an adiabatic profile and increases on the night side. This temperature inversion on the night side is caused by the redistribution of heat, which is most efficient near 2−4 km above the surface where winds are not damped by the turbulent boundary layer. The resulting mean temperature profile cannot be reproduced by the 1D model, which is too hot near the surface and too cold by more than 50 K in the upper troposphere. The discrepancy tends to decrease in the upper stratosphere.

4. Bistable climate: runaway vs cold trapping

As inferred from 1D models (e.g. Kasting et al. 1993; Selsis et al. 2007) and confirmed by our 3D simulations (not shown), if water were present in sufficient amount everywhere on the surface (aquaplanet case), both our prototype planets (Gl 581 c and HD 85512 b) would enter a runaway greenhouse state (Kasting et al. 1984). The large flux that they receive thus prevents the existence of a deep global ocean on their surfaces.

However, as is visible in Fig. 1, because of the rather inefficient redistribution, a large fraction of the surface present moderate temperatures for which either liquid water or ice would be thermodynamically stable at the surface. In the 10b case, the whole night side exhibits mean temperatures that are similar to the annual mean in Alaska. It is thus very tempting to declare these planets habitable because there is always a fraction of the surface where liquid water could flow. While tempting, one should not jump to this conclusion because of the two following points.

• Even if water is thermodynamically stable, it is always evaporating. Consequently, the water vapor is transported by the atmosphere until it condenses and precipitates (rains or snows) onto the ground in cold regions. This mechanism tends to dry warmer regions and to accumulate condensed water in cold traps. If these cold traps are permanent, or if no mechanism is physically removing the water from them – e.g. ice flow due to the gravitational flattening of an ice cap, equilibration of sea level (see Sect. 6) – this transport is irreversible.

• Because it is a strong greenhouse gas, water vapor has an important positive feedback. The presence of a small amount of water can dramatically increase surface temperatures and even trigger the runaway greenhouse instability.

Therefore, even if useful, “dry” simulations cannot be used to assess the climate of a planet where water is potentially present, and the distribution of the latter⁷. In this section we present a set of simulations that take the full water cycle into account but with a limited water inventory.

We show that, as already suggested by Abe et al. (2011), climate regimes on heavily irradiated land planets result from the competition of the two aforementioned processes that enables the existence of two stable equilibrium states: a “collapsed” state where almost all the water is trapped in permanent ice caps on the night side or near the poles, and a runaway greenhouse state where all the water is in the atmosphere. Furthermore, we show that the collapsed state can exist at much higher fluxes than inferred by Abe et al. (2011) and that there is no precise flux limit, because it can depend on both the atmosphere mass and water amount.

4.1. Synchronous case

Our numerical experiment is as follows. For our two planets and for various pressures, we start from the equilibrium state of the corresponding “dry” run described in Sect. 3.2 to which a given amount of water vapor is added. To do so, a constant mass mixing ratio of water vapor, q_v, is prescribed everywhere in the atmosphere, and surface pressure is properly rescaled to take the additional mass into account. Two typical examples of the temporal evolution of the surface temperature and water vapor content in these runs is presented in Fig. 7.

While idealized, this numerical experiment has the advantage of not introducing additional parameters to characterize the initial distribution of water. It also provides us with an upper limit on the impact of a sudden water vapor release as vapor is present in the upper atmosphere where its greenhouse effect is maximum. From another point of view, this experiment roughly mimics a planet-wide release of water by accretion of a large asteroid or planetesimal and can give us insight into the evolution of the water inventory after a sudden impact.

Because of the sudden addition of water vapor, the atmosphere is disturbed out of radiative balance. Water vapor significantly reduces the planetary albedo by absorbing a significant fraction of the stellar incoming radiation, and it also reduces the outgoing thermal flux. The temperature thus rapidly increases until a quasi-static thermal equilibrium that is consistent with the amount of water vapor is reached in a few tens of days (see Fig. 7). At the same time, water vapor condenses, forms clouds, and precipitates onto the ground near the coldest regions of the surface. As temperature increases, the amount of precipitation can decrease because the re-evaporation of the falling water droplets or ice grains becomes more efficient. The evaporation of surface water also becomes more efficient.

Then, the system can follow two different paths. If the initial amount of water vapor and the associated greenhouse effect are large enough, evaporation will overcome precipitation, and the atmosphere will stay in a stable hot equilibrium state where all the water is vaporized (with the exception of a few clouds; Fig. 7 and top panel in Fig. 8). Because this state is maintained by the greenhouse effect of water vapor and any amount of water added at the surface is eventually vaporized, this state is the usual runaway greenhouse state (Nakajima et al. 1992). In contrast, for lower water-vapor concentrations, the temperature increase is not large enough, and the water vapor collapses in the cold traps (Fig. 7 and bottom panel in Fig. 8). The amount of water vapor decreases and the surface slowly cools down again until most of the water is present in ice caps in the colder regions of the surface. An example of the ice surface density of a collapsed state at the end of the simulation is presented in Fig. 9. Because sublimation/condensation processes are constantly occurring, this ice surface density may continue to evolve, but on a much longer time scale that is difficult to model (see Wordsworth et al. 2013 for a discussion). Then, water vapor is fairly well mixed within all the atmosphere, and its mixing ratio is mainly controlled by the saturation mass mixing ratio just above the warmest region of the ice caps. Depending on the simulations, this mass mixing ratio varies between 10^-4 and 10^-5.

To investigate the factor determining the selection of one of these two regimes, we ran a large grid of simulations, varying both the total surface pressure and the initial amount of water vapor. Since the main effect of water vapor is its impact on the radiative budget, the quantity of interest is the total water vapor column mass (also called water vapor path) $\begin{array}{ccc}{\mathit{m}}_{\mathrm{v}}^{\mathit{\star }}\mathrm{\equiv }{\mathrm{\int }}_{\mathit{z}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}{\mathit{q}}_{\mathrm{v}}\hspace{0.17em}\mathit{\rho }\hspace{0.17em}\mathrm{d}\mathit{z}\mathrm{=}{\mathrm{\int }}_{\mathit{p}\mathrm{=}\mathrm{0}}^{{\mathit{p}}_{\mathrm{s}}}{\mathit{q}}_{\mathrm{v}}\hspace{0.17em}\frac{\mathrm{d}\mathit{p}}{\mathit{g}}\mathit{,}& & \end{array}$(12)which is linked to the spectral optical depth of water vapor (${\mathit{\tau }}_{\mathrm{v}}^{\mathit{\nu }}$) by the relation $\mathrm{d}{\mathit{\tau }}_{\mathrm{v}}^{\mathit{\nu }}\mathrm{=}{\mathit{\kappa }}_{\mathrm{v}}^{\mathit{\nu }}{\mathit{q}}_{\mathrm{v}}\mathit{\rho }\mathrm{d}\mathit{z}\mathrm{=}{\mathit{\kappa }}_{\mathrm{v}}^{\mathit{\nu }}\mathrm{d}{\mathit{m}}_{\mathrm{v}}^{\mathit{\star }}$.

Results for the case of Gl 581 c are summarized in Fig. 10. At a given pressure, a small initial amount of water vapor results in a collapsed case. However, there is always a critical, initial water vapor path, which is a function of the background pressure, above which the runaway greenhouse instability is triggered and develops. Because redistribution efficiency increases with atmospheric mass, cold trapping becomes weaker and the critical initial water vapor path decreases with background atmosphere surface pressure.

We also carried out simulations with a large amount of water (more than the critical water vapor path discussed above), but distributed only at the surface in the cold traps. For a background atmosphere less massive than ~1−5 bar, these simulations are found to be stable against runaway. This demonstrates that, for a given total amount of water larger than the critical water vapor path, a “moist” bistability exists in the system. Owing to the inhomogeneous insolation, the 1D notion of critical flux above which runaway greenhouse becomes inevitable is no longer valid, highlighting the dynamical nature of the runaway greenhouse that needs a large enough initial amount of water vapor to be triggered.

When the insolation is decreased, a higher amount of water vapor is needed to trigger the runaway greenhouse instability, and this at every background pressure. As would be expected, planets receiving less flux are more prone to fall into a collapse state. However, unexpectedly, results for HD 85512 b do not follow this trend. At every background pressure, a lower amount of water vapor is needed to trigger the runaway greenhouse instability, as shown in Fig. 10. This is due to the peculiar fact that, because of the higher redistribution efficiency caused by the lower rotation rate of the planet, cold traps are warmer (see Sect. 3.2 and Fig. 2 for details). Since evaporation is mainly controlled by the saturation pressure of water above the cold traps, it is more efficient, destabilizing the climate.

In addition to the work of Abe et al. (2011), which shows that runaway could be deferred to higher incoming fluxes on land planets, our simulations demonstrate that, even for a given flux, such heavily irradiated land planets can exhibit two stable climates in equilibrium. Furthermore, we have shown that the state reached by the atmosphere depends not only on the impinging flux, but also on the mass of the atmosphere, on the initial water distribution, and on the circulation regime (controlled by the forcing, the size of the planet, and its rotation rate). This bistability is due to the dynamical competition between cold trapping and greenhouse effect. An important consequence of this bistability is that the climate present on planets beyond the inner edge of the habitable zone will not only depend on their background atmosphere’s mass and composition, but also on the history of the water delivery and escape, as discussed below.

4.2. The 3:2 resonance

Because of the strong tidal interaction with their host star and of the small angular momentum contained in the planetary spin, rocky planets within or closer than the habitable zone of M and K stars are thought to be rapidly synchronized. Indeed, the standard theory of equilibrium tides (Darwin 1880; Hut 1981; Leconte et al. 2010) predicts that a planet should synchronize on a timescale equal to $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{syn}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}{\mathit{r}}_{\mathrm{g}}^{\mathrm{2}}\frac{{\mathit{M}}_{\hspace{0.17em}\mathrm{p}}\hspace{0.17em}{\mathit{a}}^{\mathrm{6}}}{\mathit{G}{\mathit{M}}_{\mathit{\star }}^{\mathrm{2}}{\mathit{R}}_{\hspace{0.17em}\mathrm{p}}^{\mathrm{3}}{\mathit{k}}_{\hspace{0.17em}\mathrm{2}\mathit{,}\mathrm{p}}\hspace{0.17em}\mathrm{\Delta }{\mathit{t}}_{\hspace{0.17em}\mathrm{p}}}\mathit{,}& & \end{array}$(13)where r_g is the dimensionless gyration radius (${\mathit{r}}_{\mathrm{g}}^{\mathrm{2}}\mathrm{=}\mathrm{2}\mathit{/}\mathrm{5}$ for a homogeneous interior; Leconte et al. 2011), G is the universal gravitational constant, k_2,p the tidal Love number of degree 2, and Δt_p a time lag that characterizes the efficiency of the tidal dissipation into the planet’s interior (the higher Δt_p, the higher the dissipation). Other variables are defined in Table 1. While the magnitude of the tidal dissipation in massive terrestrial planets remains largely unconstrained (Hansen 2010; Bolmont et al. 2011), one can have a rough idea of the orders of magnitude involved by using the time lag derived for the Earth, k_2,p Δt_p = 0.305 × 629 s (Neron de Surgy & Laskar 1997). This yields a synchronization timescale of 25 000 yr for Gl 581 c and 10 Myr for HD 85512 b, preventing these planets from maintaining an initial obliquity and fast rotation rate (Heller et al. 2011).

As demonstrated by Makarov et al. (2012) in the case Gl 581 d, which is more eccentric and farther away from its parent star than Gl 581 c, tidal synchronization is not the only possible spin state attainable by the planet. As for Mercury, if the planet started from an initially rapidly rotating state, it could have been trapped in multiple spin orbit resonances during its quick tidal spin down because of its eccentric orbit. Because Gl 581 c also has a non-negligible eccentricity of 0.07 ± 0.06 (Forveille et al. 2011), it has almost a 50% probability of being captured in 3:2 resonance (depending on parameter used) before reaching the synchronization (Makarov & Efroimsky, priv. comm., 2012). Higher order resonances are too weak.

To assess the impact of an asynchronous rotation, we performed the analysis described above in the 3:2 resonance case. The rotation rate of the planet is thus three halfs of the orbital mean motion and the solar day lasts two orbital periods. As expected, depending on the thermal inertia of the ground and on the length of the day, the night side is hotter than in the synchronized case. More importantly, there is no permanent cold trap at low latitudes. However, as for Mercury (Lawrence et al. 2013), because of the low obliquity, a small area near the poles remains cold enough to play this role, and ice can be deposited there. As the magnitude of the large-scale transport of water vapor is lower towards the pole compared to the night side (zonal winds are stronger than meridional ones), the trapping is less efficient. As a result, the transition towards the runaway greenhouse state discussed in the previous section is triggered by a lower initial water vapor path. The timescale needed to condense all the water vapor near the poles is also much longer.

5. Evolution of the water inventory

How stable are the two equilibrium states described above in the long term? It is possible that the climate regime exhibited by an irradiated land planet changes during its life depending on the rate of water delivery and escape of the atmosphere and/or water: triggering runaway when water delivery (or water ice sublimation by impacts) is too large and collapsing water when either the background atmosphere or the water vapor mass decreases below a certain threshold. We now briefly put some limits on the fluxes involved.

5.1. Atmospheric escape

When water vapor is present in the upper atmosphere, XUV radiations can photo dissociate it and the light hydrogen atoms can escape⁸. An upper limit on this escape flux is provided by the energy limited flux (Watson et al. 1981) approximated by $\begin{array}{ccc}{\mathit{F}}_{\mathrm{el}}\mathrm{=}\mathit{ϵ}\frac{{\mathit{R}}_{\hspace{0.17em}\mathrm{p}}{\mathit{F}}_{\mathrm{XUV}}}{\mathit{G}{\mathit{M}}_{\hspace{0.17em}\mathrm{p}}}\mathrm{\left(}\mathrm{kg}\hspace{0.17em}{\mathrm{m}}^{-2}\hspace{0.17em}{\mathrm{s}}^{-1}\mathrm{\right)}\mathit{,}& & \end{array}$(14)where F_XUV is the energetic flux received by the planet and ϵ the heating efficiency. The problem then lies in estimating the flux of energetic photons, which varies with time (see Selsis et al. 2007 for a discussion). Estimates of F_el are shown in Fig. 11 where the XUV flux parametrization of Sanz-Forcada et al. (2010) and a heating efficiency of ϵ = 0.15 have been used⁹. In some cases, this formula can also be used to estimate the escape of the background atmosphere itself.

In atmospheres where the amount of water vapor in the upper atmosphere is limited, the escape flux can be limited by the diffusion of water vapor in air. Following Abe et al. (2011), the diffusion-limited escape flux is written as $\begin{array}{ccc}{\mathit{F}}_{\mathrm{dl}}\mathrm{=}{\mathit{f}}_{\mathrm{str}}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\right)}\hspace{0.17em}{\mathit{b}}_{\mathrm{ia}}\frac{\mathrm{(}{\mathit{m}}_{\mathrm{a}}\mathrm{-}{\mathit{m}}_{\mathrm{i}}\mathrm{)}\hspace{0.17em}\mathit{g}}{{\mathit{k}}_{\hspace{0.17em}\mathrm{B}}{\mathit{T}}_{\mathrm{str}}}\mathit{,}& & \end{array}$(15)where f_str(H₂) is the mixing ratio of hydrogen in all forms in the stratosphere, T_str is the temperature there, and m_a and m_i are the particle mass of air and the species considered, and b_ia is the binary diffusion coefficient of i in air (for H₂ in air b_ia = 1.9 × 10²¹(T/300 K)^0.75 m^-1s^-1).

In the collapsed state, the longitudinal and vertical cold traps are efficient. Temperature and water-vapor mixing ratios are thus low in the stratosphere (T_str ≈ 220 K and f_str(H₂) < 10^-4; see Fig. 11). In the runaway state, the stratosphere temperature is much higher (≈300  K), and water is rather well mixed in the atmosphere. The mixing ratio can be directly calculated from the total water amount. To show an example, we plotted the case of f_str(H₂) ≈ 10^-1 in Fig. 11.

Of course, to have a comprehensive view of the problem, hydrodynamic and nonthermal (due to stellar winds) escape should also be accounted for. However, some conclusions can be made. First, it seems that escape of the background atmosphere could be efficient during the first few hundred million years to one billion year. This supports the possibility of a thin atmosphere. Second, as expected, escape of water vapor is always greater in the runaway state because of the higher temperature and humidity in the higher atmosphere (Kasting et al. 1984). Third, in the collapse state, cold traps limit the humidity of the upper atmosphere, so escape is most likely diffusion-limited to a rather low value.

Interestingly enough, the sudden increase of water loss when runaway is triggered can provide a strong negative feedback on the water content of land planets beyond the classical inner edge of the habitable zone. Indeed, at least at early times where the XUV flux is strong enough, if a significant fraction of water is accreted suddenly and triggers the runaway, water escape will act to lower the water vapor amount until collapse ensues. This should happen when the atmospheric water vapor content is within order of magnitude of, although lower than, the critical water column amount shown in Fig. 10¹⁰. Then, cold-trapped water has a much longer lifetime. This could particularly play a role during the early phases of accretion, where a runaway state is probable considering the strong energy and water vapor flux due to the accretion of large planetesimals (Raymond et al. 2007), a possibly denser atmosphere, and a faster rotation. If the planet receives a moderate water amount, the water could be lost in a few tens to a hundred million years, leaving time to the planet to cool down and accrete a significant amount of ice through accretion of comets and asteroids later on, either progressively or during a “late heavy bombardment” type of event for example (Gomes et al. 2005). The possibility a such a scenario is supported by the debris disk has been recently imaged in the Gl 581 system (Lestrade et al. 2012).

5.2. Cold trapping rate

A constraint on the accumulation of ice from such a late water delivery is that the delivery rate must be low enough to allow trapping to take place without triggering the runaway. To first order, our simulations presented in Sect. 4 can be used to quantify this cold trapping rate.

We consider a case where the water vapor collapses. As seen in Fig. 7, except for the very early times of the simulations, the water vapor content of the atmosphere exponentially tends towards its (low) equilibrium value. We can thus define an effective trapping rate $\mathit{ṁ}\begin{array}{c}\mathit{\star }\\ \mathrm{v}\end{array}\mathrm{=}{\mathit{m}}_{\mathrm{v}}^{\mathit{\star }}\mathit{/}\mathrm{(}\mathrm{2}\hspace{0.17em}{\mathit{\tau }}_{\mathrm{c}}\mathrm{)}$ (in kg s^-1m^-2), where ${\mathit{m}}_{\mathrm{v}}^{\mathit{\star }}$ is the initial water-vapor column mass and τ_c the time taken by the atmosphere to decrease the water vapor amount by a factor of two. This is basically a rough estimate of the critical rate of water delivery above which runaway greenhouse would ensue.

The maximum trapping rate obtained in our simulations as a function of the background atmosphere surface pressure is given in Fig. 12. As expected, trapping is less efficient when the atmospheric mass is higher because cold traps are warmer. The exponential decrease can be understood considering that re-evaporation efficiency linearly decreases with saturation pressure that exponentially decreases with temperature following Clausius-Clapeyron law. Then, the maximum amount of ice that can be accumulated really depends on the duration of the phase during which water is delivered to the planet, but the quite high trapping rate found in our simulations suggest that if water accretion occurs steadily, the final water inventory may not be limited by the trapping efficiency but by the rate of water delivery after the planet formation¹¹.

6. On the possibility of liquid water

Our simulations show that, if the atmosphere is not too thick and the rate of water delivery high enough (while still avoiding runaway), a large amount of water can accumulate in the cold traps. However, the low temperatures present there prevent liquid water stability. Are these land planets possibly habitable? Can liquid water flow on their surface for an extended period of time?

A possibility is that both volcanism and meteoritic impacts could be present on the night side, creating liquid water, as has been recently proposed to explain the evidence of flowing liquid water during the early Martian climate (Segura et al. 2002; Halevy et al. 2007; Wordsworth et al. 2013). Liquid water produced that way is, however, only episodic.

In this section, we briefly explain various mechanisms that could produce long-lived liquid water (as summarized in Fig. 13). In particular, we argue that, considering the thickness of the ice cap that could be present, gravity-driven ice flows similar to those taking place on Earth could transport ice toward warm regions where it would melt. Liquid water could thus be constantly present near the edge of the ice cap. We also discuss the possibility of having subsurface liquid water produced by the high pressures reached at the bottom of an ice cap combined with a significant geothermal flux.

We also explored the quasi-snowball scenario proposed by Selsis et al. (2007) for Gl 581 c where a low equilibrium temperature is provided by a very high ice albedo (>0.9) and temperature slightly exceeds 0°C at the substellar point. Notwithstanding the difficulty of having such high albedo ices around such a red star (Joshi & Haberle 2012), our simulations tend to show that this scenario can be discarded. Indeed, the very strong positive feedback of both water vapor and ice albedo render this solution very unstable and no equilibrium state was found.

6.1. Ice flows

As described in Sect. 4, if water accumulates on the surface, the atmosphere continuously acts to transport it towards the coldest regions where it is the most stable (more precipitations and less evaporation). As this accumulation seems possible mostly for low background atmospheric pressures, locations always exist on the surface where temperature is below freezing (see Figs. 1 and 2) and only ice is stable.

If atmospheric transport were the only transport mechanism of water on Earth, the water would also accumulate onto the poles, and liquid water would only be present when the Sun melted the top of the ice cap during summer. However, other processes come into play to limit the thickness of the polar ice caps. Most importantly, because ice is not completely rigid, gravity tends to drain ice from the upper regions (or regions where the ice cap is thicker) towards lower regions creating ice flows, as seen in Greenland and Antarctica, among other places. These ice flows physically remove ice from the cold traps to transport it towards regions where it can potentially melt, at the edge of glaciers, for example.

Assuming that enough water has been able to accumulate on the surface without triggering the runaway greenhouse instability, such ice flows should eventually occur. Thus, liquid water produced by melting should be present when the ice is drained to the day side where temperatures are above 0°C. Quantifying the amount of ice necessary and the average amount of liquid water is, however, much more uncertain, and maybe not relevant considering our present (lack of) knowledge concerning extrasolar planets surface properties and topography. It should be said that in Greenland, an ice sheet of around 2 km is sufficient to transport ice over distances of about 500−1000 km and create ice displacements up to 2−3 km per year near glaciers (Rignot et al. 2011; Rignot & Mouginot 2012). Since gravity (which is the main driver) and typical horizontal length scales both scale linearly with radius if density is kept constant, this order of magnitude should remain valid for more massive land planets.

As discussed in Sect. 5, ice thicknesses of a few tens of kilometers could possibly be reached. It is thus possible to imagine a stationary climate state where i) a thick enough ice cap produces ice flows and liquid water at the edge of warmer regions; ii) this liquid water evaporates; and iii) the resulting water vapor is transported back near the coldest regions where it precipitates. A question that remains is whether the additional amount of water vapor that will be released in the atmosphere can be large enough to trigger the runaway greenhouse instability that would cause the melting and evaporation of the whole ice cap.

To assess the possibility of this scenario, we ran another set of simulations with a thick ice cap. Because we did not want to add any free parameter (such as the total amount of ice) and we want a worst-case scenario (the most favorable to triggering the runaway), we modified our base model so that the surface could act as an infinite reservoir of water whose property (liquid or solid) depends solely on the ground temperature (above or below 0 °C). In particular, we did not consider the latent heat needed to melt water that would have a stabilizing effect on the climate by lowering the water-vapor evaporation rate.

Another problem comes from the fact that freezing regions change as the amount of water vapor increases in the atmosphere. In our numerical experiment, we thus updated the surface properties at each time step as follows. A dry surface grid cell becomes wet (a reservoir of water) if it has a neighbor with ice on its surface: the ice cap extends. Since we want to model a case where only the edge of the ice cap melts, if a wet surface grid cell has no neighbor where water ice is present, it is dried. Then the model is initialized from the dry state from Sect. 3.2 where a seed of ice has been implanted in the cold trap, and ran until either a stationary equilibrium state is reached or the surface is completely dried by the runaway greenhouse effect. Test simulations show that the final state reached does not depend on the initial distribution of ice.

For Gl 581 c, our simulations show that above five bars, runaway is inevitable. Below one bar, a stationary state is always found where ice and liquid water are both present on the surface at the same time (see Fig. 14). Depending on the topography near the ice cap edge, this could correspond to a situation where lakes are constantly replenished by the melting water or where rivulets moisten a more extended region. Interestingly, as shown in the middle and left hand panels of Fig. 14, ice does not necessarily cover all the night side because of the radiative effect of water vapor which tends to warm the surface (see Fig. 15), and the inhomogeneous redistribution of energy discussed in Sect. 3.2. In particular, melting first occurs eastward from the day side near the equator where the surface is heated by the hot air coming from the day side through the jet, as seen in the 400 mbar case of Fig. 14.

For HD 85512 b, again, we find that runaway occurs at lower background surface pressures. The transition occurs between 400 and 500 mb. Also, because of the more homogeneous temperature distribution of the night side, the transition is not as progressive as shown in Fig. 14 for Gl 581 c. If the 400 mb case reaches a steady state with a full ice coverage of the night side, the 500mb case is almost dry.

These simulations demonstrate that, even in our most stringent case, the scenario of a thick ice cap, vigorous ice flows and a significant amount of liquid water could be present without triggering the runaway greenhouse effect. Indeed, in these simulations, most of the grid cells next to the ice cap are covered by liquid water (see Fig. 14), and the climate remains stable. Given our resolution, this corresponds to a wet region with a width of about 1000 km, which is probably an overestimation of the horizontal extension of the region where melting will occur. Interestingly, when there is liquid water on the surface in our model, part of the wet surface is always receiving some sunlight, because a state with liquid water only on the night side is not stable. In this case, water vapor condenses and the water freezes . One must bear in mind that, without further information, we did not include any runoff of liquid water at the surface, assuming that the small amount of flowing water would be evaporated before crossing these 1000 km. A more complete model of the ice sheet, ice flows and runoff would be necessary to completely settle this matter.

6.2. Subsurface liquid water?

As for icy moons such as Europa and Enceladus (Sotin et al. 1998; Pappalardo et al. 1999) and subsurface lakes such as the Vostok lake on Earth (Oswald & Robin 1973), if a thick ice cap exists, the possible presence of a significant amount of subsurface liquid water at the bottom of the ice layer is also worth considering. Indeed, because of the decrease in ice melting temperature with pressure and the increase in temperature with depth caused by any geothermal flux, it is easier, in a way, to melt water below a thick ice sheet than on the surface.

To assess this possibility, we estimated the ice thickness needed to melt water beneath it as function of the geothermal flux. Considering that vertical energy transport in the subsurface is mainly performed by thermal conduction, and that quasi-static equilibrium is reached, the temperature gradient is given by $\begin{array}{ccc}\frac{\mathrm{d}\mathit{T}}{\mathrm{d}\hspace{0.17em}\mathit{z}}\mathrm{=}\mathrm{-}\frac{{\mathit{\phi }}_{\mathrm{g}}}{{\mathit{K}}_{\mathrm{th}}}\mathit{,}& & \end{array}$(16)where φ_g is the (positive) geothermal flux and K_th is the thermal conductivity of water ice that ranges between 2.22 W/m/K at 0°C and 3.48 W/m/K at − 100 [entity ! #x20 ! ] °C. Thanks to hydrostatic equilibrium, and assuming a constant density and conductivity for ice (here 920 kg/m³ and 0.3 W/m/K respectively), the temperature pressure profile is given by $\begin{array}{ccc}\mathit{T}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}{\mathit{T}}_{\mathrm{s}}\mathrm{+}\frac{{\mathit{\phi }}_{\mathrm{g}}}{{\mathit{\rho }}_{\mathrm{ice}}\mathit{g}{\mathit{K}}_{\mathrm{th}}}\left(\mathit{p}\mathrm{-}{\mathit{p}}_{\mathrm{s}}\right)\mathit{,}& & \end{array}$(17)where T_s and p_s are the local surface temperature and pressure and ρ_ice the mass density of ice. Following Wagner et al. (1994), we parametrize the melting curve of water ice by $\begin{array}{ccc}{\mathit{p}}_{\mathrm{melt}}\mathrm{\left(}\mathit{T}\mathrm{\right)}& \mathrm{=}& {\mathit{p}}_{\mathrm{0}}[\mathrm{1}\mathrm{-}\mathrm{0.626}\mathrm{\times }{\mathrm{10}}^{\mathrm{6}}\left(\mathrm{1}\mathrm{-}{\left(\frac{\mathit{T}}{{\mathit{T}}_{\mathrm{0}}}\right)}^{-3}\right)\\ & & \mathrm{+}\hspace{0.17em}\mathrm{0.197135}\mathrm{\times }{\mathrm{10}}^{\mathrm{6}}\left(\mathrm{1}\mathrm{-}{\left(\frac{\mathit{T}}{{\mathit{T}}_{\mathrm{0}}}\right)}^{\mathrm{21.2}}\right)]\mathit{,}\end{array}$(18)where T₀ = 273.16 K and p₀ = 611.657 Pa are the temperature and pressure of the triple point of water. This parametrization is valid only between 251.165 and 273.16 K because liquid water is not stable below this temperature at any pressure.

Then, for a given geothermal flux and surface temperature, subsurface liquid water can exist if the curves parametrized by Eqs. (17) and (18) intersect each other. By solving this implicit equation, we can find the pressure at which melting occurs and define the minimal ice cap thickness needed to reach this pressure, l_ice = (p − p_s)/(ρ_iceg) (for the atmosphere considered, p_s can safely be neglected in the above equations). Numerical estimates of this thickness as a function of surface temperature and geothermal flux are shown in Fig. 16 for a typical gravity of 18 m/s².

One can see that even if no geothermal flux is present, subsurface liquid water could exist for warm enough regions. The thickness needed would be rather large (~3−10 km), which seems unlikely. Depending on the amount of internal flux, however, the minimal thickness can be much lower. On Earth, radioactive decay and release of primordial heat generate an overall averaged geothermal flux of 92 mW/m² (Davies & Davies 2010). If we assume that the rate of volumetric heating is constant, this mean flux should scale linearly with the planet’s radius, and our more massive planet should produce between 1.6 and 1.8 times this flux. With such a flux, the thickness needed decreases below 1.5 km in the coldest regions but could be as low as 300 m.

Finally, another important source of internal flux for close-in planets is the friction created by the bodily tides raised by the host star. An estimate of this tidal dissipation in the synchronous case is given by Eq. (14) of Leconte et al. (2010) (see also Hut 1981). This heating is a function of the orbital eccentricity of the planet and of the efficiency of the tidal dissipation into the planet’s interior (the imaginary part of the Love number, k_2,p Δt_p; see Sect. 4.2). For an Earth-like dissipation and eccentricities shown in Table 1, the flux can be estimated to be ~600 mW/m² for Gl 581 c and would contribute in large part to the interior energy budget¹². As shown in Fig. 16, the thickness needed to produce liquid water would be much lower (<300  m). Because the measured eccentricity is highly uncertain, we computed the eccentricity needed to produce the internal fluxes shown in Fig. 16 (top scale of the figure). However, one must bear in mind that the dissipation constant is also highly uncertain, and the tidal flux computed here could easily change by an order of magnitude. For HD 85512 b the tidal flux is much lower (0.1 mW/m²) because the planet is much farther away from its parent K star.

Extended reservoirs of subsurface water can thus not be ruled out. The ice thickness needed, if realistic internal fluxes are present, is much lower than the one that could potentially be accumulated as discussed in Sect. 5.

7. Observable signature

Can observation allow us to distinguish between the different climate regimes presented here? To start answering this question, we consider various possible observables. If a given target’s flux is too small to use spectroscopy, broad band photometry (as done for Kepler-10 b (Batalha et al. 2011) and CoRoT-7 b (Léger et al. 2009)) can already provide constraints. Indeed, as shown in Fig. 17, the presence of water vapor in the atmosphere of a planet orbiting an M dwarf strongly decreases its Bond albedo because many water spectral bands peak at the same wavelengths as the stellar spectrum. However, while this method can hint at the presence of a strong absorber in the atmosphere of the planet observed – hinting at the presence of an atmosphere itself – a strong degeneracy exists between the various possible gases, so spectroscopy is needed. This signature will, of course, decrease for more massive, bluer stars.

To go further, we computed synthetic spectra of various models depicted above. We used the tool developed by Selsis et al. (2011) that can model the spectral phase curve of a planet in an arbitrary geometrical configuration directly from the output of the GCM code. Here we focus only on the spectrum that could be observed during the secondary transit of a transiting analog of Gl 581 c at 10 parsecs. An example of the spectra computed for the 200 mb runaway and collapse cases presented in Figs. 7 and 8 are shown in Fig. 18.

Although the intensity of the bolometric signal is weak, one can see that spectral differences between a runaway and a collapse case are striking. First, one can easily spot the absorption bands of water vapor in the reflection spectrum in the runaway case, which are not present in the dry collapse case. Indeed, when water vapor becomes a significant constituent of the atmosphere, the latter becomes thick in the H₂O molecular bands even in the near infrared. A large fraction of the incoming stellar radiation is thus absorbed in those bands before being reflected by the surface or Rayleigh scattering. Second, in the emission spectrum, one can see that the planet in runaway can only radiate trough the H₂O window regions around 4 and 10 μm. The surface that is heated by the strong resulting greenhouse effect thus radiates at shorter wavelengths.

If sufficient accuracy is met (see Selsis et al. 2011 for a more detailed discussion), further observation may thus provide strong constraints on the climate regime reached on such strongly irradiated land planets, since the presence of strong H₂O spectral features should be correlated with the shape of the thermal emission.

8. What’s next?

As for Gl 581 c or HD 85512 b, the habitability of planets receiving more flux than the classical runaway greenhouse has often been discarded on the basis of 1D models that predict that an ocean cannot be stable under these conditions.

Here, we have demonstrated that because of the inhomogeneity of the impinging flux and of the natural ability of the atmosphere to transport condensable species toward colder regions where it can be captured, the presence of surface water cannot be ruled out on the sole basis of the mean absorbed stellar flux. In fact, because of the strong positive feedbacks of water vapor, two stable equilibrium climate regimes exist for the atmosphere: a runaway greenhouse state and a state where all the water is trapped on the surface in cold regions. A key finding of this study is that the state in which the atmosphere settles depends not only on the incoming flux, but also on many other factors: atmospheric mass, water content and initial distribution, spin state, and circulation regime.

If it is not easy to set clear limits on the domain of existence of this bistable regime, however a few trends do emerge. First, the mechanism described above strongly relies on the existence of an efficient permanent (or slowly evolving) cold trap. The trapping in a spin orbit resonance or a very low obliquity thus seem to strongly favor the onset of the bistability. Second, except for changes in the circulation pattern, the domain of existence of the collapse state narrows down in the pressure/initial water vapor amount parameter space when the incoming flux increases. These two arguments suggest that the domain of the bistability should become smaller when the star becomes more massive because regions where tidal interactions causing resonances to occur and damping obliquity move at higher fluxes. Again, it is difficult to set a clear limit in terms of a spectral type above which this bistability cannot occur. Let us just say that our simulations support the existence of a bistability around M and K dwarfs, and that, if not for its thick atmosphere, Venus’ climate would be very close to the one of HD 85512 b, extending this statement to G dwarfs. The presence of ice on the poles of Mercury (Lawrence et al. 2013) also gives observational support to this scenario, at least within the limit of a vanishing atmosphere. The validity of our results thus seems to extend to a broad range of situations.

Interestingly, because all the parameters affecting the occurrence of a bistable state change during the life of the planet, the present climate of a given object may depend on the evolutionary path that it followed. In particular, the history of the water delivery and of the atmospheric escape may play a crucial role in determining whether these objects are dry rocks surrounded by a hot steam atmosphere or else harbor vast ice deposits on the night side and near the poles.

Another finding is that the atmospheric circulation present on synchronous extrasolar planets does not necessary settle into a state of super-rotation. Since the super-rotation seems to be powered by planetary-scale Rossby and Kelvin waves, we find that it weakens when the wavelength of these waves, i.e. the equatorial Rossby deformation radius, significantly exceeds the planetary radius. For small and/or very slowly rotating planets (synchronously rotating on a wide orbit, for example) the circulation seems to be dominated by a symmetrical stellar to anti stellar circulation at high altitudes.

To move further ahead, observational data are needed. Fortunately, exoplanet-characterizing observatories, such as EChO (Tinetti et al. 2012) and JWST, should soon be available. And even if the characterization of Earth analogs around M dwarfs proves to be very challenging, we will learn a lot from the observation of planets beyond the runaway greenhouse limit, such as transiting versions of the ones described here and even hotter. Combining spectral observations of the primary and secondary transits (and possibly phase curves) of a number of strongly irradiated terrestrial exoplanets should indeed help us unravel the variety of atmospheric masses, compositions, and histories, which are, at the present time, the most poorly constrained parameters used in exoplanetary atmosphere modeling.

---

Acknowledgments

J.L. thanks V. Makarov and M. Efroimsky for providing calculations on the spin state of Gl 581 c and O. Grasset for an inspiring discussion concerning icy moons and subsurface oceans. We thank our referee, A. Showman, for his insightful comments. J.L. received funding from the DIM ACAV. F.S. acknowledges support from the European Research Council (ERC Grant 209622: E₃ARTHs).

References

All Tables

All Figures