GLIESE 581D IS THE FIRST DISCOVERED TERRESTRIAL-MASS EXOPLANET IN THE HABITABLE ZONE

The local red dwarf Gliese 581 (20.3 lt-yr from the Sun, M = 0.31 M_Sun, L = 0.0135 L_Sun, spectral type M3V) (Hawley et al. 1997) has received intense interest over the last decade due to the low mass exoplanets discovered around it. As of early 2011, it has been reported to host up to six planets (Udry et al. 2007; Mayor et al. 2009; Vogt et al. 2010). One of these, GJ581g, was announced in 2010 September and estimated to be in the habitable zone (the orbital range in which a planet's atmosphere can warm the surface sufficiently to allow surface liquid water; Kasting et al. 1993; Pierrehumbert 2011). However, its discovery has been strongly disputed by other researchers, including the team responsible for finding the other four planets in the system (Kerr 2010; Tuomi 2011). For the moment, therefore, GJ581g remains unconfirmed.

GJ581d, in contrast, which was first discovered in 2007 and has a minimum mass between 5.6 and 7.1 M_Earth, has now been robustly confirmed by radial velocity observations (Udry et al. 2007; Mayor et al. 2009; Vogt et al. 2010). Due to its greater distance from the host star, GJ581d was initially regarded as unlikely to have surface liquid water unless strong warming mechanisms due to, e.g., CO₂ clouds (Forget & Pierrehumbert 1997; Selsis et al. 2007) present in its atmosphere. Recently, simple one-dimensional radiative-convective studies (Wordsworth et al. 2010b; von Paris et al. 2010; Kaltenegger et al. 2011) have suggested that a dense atmosphere could provide a significant greenhouse effect on GJ581d. However, the planet's tidal evolution poses a key problem for its habitability.

As it is most likely either in a pseudo-synchronous state with a rotation period that is a function of the eccentricity, or in spin-orbit resonance like Mercury in our solar system (Leconte et al. 2010; Heller et al. 2011), it should have extremely low insolation at its poles and possibly a permanent night side. Regions of low or zero insolation on a planet can act as cold traps where volatiles such as H₂O and CO₂ freeze out on the surface. A few previous studies (Joshi et al. 1997; Joshi 2003) have examined atmospheric collapse in three dimensions with simplified radiative transfer, but only for Earth-like atmospheric pressures or lower (0.1 to 1.5 bar). For low values of stellar insolation and large planetary radii, even dense CO₂ atmospheres will be prone to collapse, which could rule out a stable water cycle altogether for a super-Earth like GJ581d. To conclusively evaluate whether GJ581d is in the habitable zone, therefore, three-dimensional simulations using accurate radiative transfer are necessary.

Here, we present global climate model (GCM) simulations we performed to assess this issue. In Section 2, we describe the model used. In Section 3 we describe our results, and in Section 4 we discuss implications and propose future observational tests for the simulated habitable scenarios.

In our simulations, we made the initial hypothesis that GJ581d has a climate dominated by the greenhouse effects of CO₂ and/or H₂O, as is the case for all rocky planets with atmospheres in the solar system (Venus, Earth, and Mars). To assess the influence of water on the climate independently, we considered two classes of initial conditions: a rocky planet with no water and an ocean planet, where the surface is treated as an infinite water source. CO₂ was taken as the primary constituent of the atmosphere and H₂O was allowed to vary freely, with surface ice/liquid and cloud formation (including radiative effects) taken into account for either gas when necessary. Restricting the composition of the atmosphere to two species in this way allows us to determine conservative conditions for habitability, as it neglects the warming due to other greenhouse gases like CH₄ or buffer gases like N₂ or Ar (von Paris et al. 2010; Goldblatt et al. 2009; Li et al. 2009).

The simulations were performed using a new type of GCM that we developed specifically for exoplanet and paleoclimate studies. It uses radiative transfer data generated directly from high resolution spectra, which allows the accurate simulation of climates for essentially any atmospheric cocktail of gases, aerosols and clouds for which optical data exists. The dynamical core of the model was adapted from the LMD Mars GCM, which uses an enstrophy and total angular momentum conserving finite difference scheme (Sadourny 1975; Forget et al. 1999). Scale-selective hyperdiffusion was used in the horizontal plane for stability. The planetary boundary layer was parameterized using the method of Mellor & Yamada (1982) to calculate turbulent mixing, with the latent heat of H₂O also taken into account in the surface temperature calculations when necessary. A standard roughness coefficient of z₀ = 1 × 10⁻² m was used for both rocky and ocean surfaces for simplicity, although we verified that our results were insensitive to variations in this parameter. Spatial resolution of 32 × 32 × 20 in longitude, latitude, and altitude was used for all simulations; we performed one comparison test at the highest rotation rate with 64 × 64 × 20 resolution and found that the differences were small.

Our radiative transfer scheme was similar to that which we developed previously for one-dimensional simulations (Wordsworth et al. 2010a, 2010b). For a given mixture of atmospheric gases, we computed high-resolution spectra over a range of temperatures, pressures, and gas mixing ratios. For this study we used a 6×9×7 temperature, pressure, and H₂O volume mixing ratio grid with values T = {100, 150, ..., 350} K, p = {10⁻³, 10⁻², ..., 10⁵} mbar, and $q_{{\rm H_2O}}=\lbrace 10^{-7}, 10^{-6}, \ldots, 10^{-1} \rbrace$, respectively. The correlated-k method was used to produce a smaller database of coefficients suitable for fast calculation in a GCM. The model used 38 spectral bands in the longwave and 36 in the shortwave, and 16 points for the g-space integration, where g is the cumulated distribution function of the absorption data for each band. In most simulations CO₂ was assumed to be the main constituent of the atmosphere, except the locally habitable ice planet experiments (see Section 3), where N₂ was used. CO₂ collision-induced absorption was included using a parameterization based on the most recent theoretical and experimental studies (Wordsworth et al. 2010a; Gruszka & Borysow 1998; Baranov et al. 2004). For the stellar spectrum, we used the Virtual Planet Laboratory AD Leo data (Segura et al. 2005). AD Leo is an M-class star with effective temperature T_eff = 3400 K, which is acceptably close to the most recent estimates of T_eff = 3498 ± 56 K (von Braun et al. 2011) for the purposes of climate modeling. A two-stream scheme (Toon et al. 1989) was used to account for the radiative effects of both clouds and Rayleigh scattering, as in Wordsworth et al. (2010b).

In the water cycle and cloud modeling, care was taken to ensure that the parameterizations used were based on physical principles and not tuned to Earth-specific conditions. When this was not possible (as for e.g., the density of condensable cloud nuclei in the atmosphere N_c), we tested the sensitivity of our results to variations in those parameters. Three tracer species were used in our simulations: CO₂ ice, H₂O ice, and H₂O vapor. Tracers were freely advected in the atmosphere, subject to changes due to sublimation/evaporation and condensation and interaction with the surface. For both gases, condensation was assumed to occur when the atmospheric temperature dropped below the saturation temperature. Local mean CO₂ and H₂O cloud particle sizes were determined from the amount of condensed material and the density of condensable nuclei N_c. This parameter was set to 10⁵ kg⁻¹ in most of our simulations; we tested the effect of varying it over the range 10⁴–10⁶ kg⁻¹ and found that the maximum difference in mean surface temperature after 60 orbits was less than 5 K. As a further test of the robustness of our results, we also performed some tests with cloud radiative effects removed altogether (see Section 3).

Ice particles of both species were sedimented according to Stokes law (Forget et al. 1999). Below the stratosphere, adjustment was made to relax temperatures due to convection and/or condensation of CO₂ and H₂O. For H₂O, moist and large-scale convection were taken into account following Manabe & Wetherald (1967). Precipitation of H₂O due to coagulation was also included using a simple threshold parameterization (Emanuel & Ivkovi-Rothman 1999).

On the surface, the local albedo varied according to the composition (rocky, ocean, CO₂ or H₂O ice; see Table 1). In the wet simulations, ice formation (melting) was assumed to occur when the surface temperature was lower (higher) than 273 K, and temperature changes due to the latent heat of fusion were taken into account. In all cases, the simulations were run until collapse/glaciation occurred or steady states of thermal equilibrium were reached. The time taken to reach thermal equilibrium can be estimated from the atmospheric radiative relaxation timescale (Goody & Yung 1989)

$$\begin{equation} \tau _r = \frac{c_p p_{s}}{\sigma g T_{e}^3}, \end{equation} \tag{ 1 }$$

where c_p, p_s, g, T_e and σ are the specific heat capacity of the atmosphere, the mean surface pressure, the surface gravity, the atmospheric emission temperature, and the Stefan–Boltzmann constant, respectively. Taking T_e = 200 K, c_p = 850 J K⁻¹ kg⁻¹, p_s = 30 bar, and g = 16.6 m s⁻² yields τ_r ∼ 4000 Earth days (60 GJ581d orbits). This was close to the timescales we observed in the model by plotting time series of mean surface temperature.

We used the minimum mass for GJ581d given by Mayor et al. (2009) instead of the smaller value (M_min = 5.6 M_Earth) proposed by Vogt et al. (2010). We took the actual mass of GJ581d to be $M=M_{{\rm min}} \slash \mbox{sin } 60^\circ = 8.2 \;M_{{\rm Earth}}$, given that the statistically most probable value for the inclination angle is 60°. The radius and gravity for rocky and ocean/ice cases (Table 1) were then determined from theoretical models (Sotin et al. 2007). In the latter case, the assumed bulk composition of the planet was 50% H₂O. Model tests using M = M_min and M = 1.6M_min (the latter value comes from dynamical stability considerations (Mayor et al. 2009)) did not reveal significant differences from our main results.

To produce emission spectra, we recorded top-of-atmosphere longitude–latitude maps of outgoing fluxes computed by the GCM over one orbit. The inclination angle of the orbit relative to the observer was assumed to be 60°. An isotropic (Lambertian) distribution of specific intensities at the top of the atmosphere was assumed; comparison with a line-by-line radiative transfer code at wavelengths where the limb-darkening was most pronounced revealed that the disk-integrated flux error due to this effect was below 5%.

We performed simulations with 5, 10, 20, and 30 bar atmospheric pressure and 1:1, 1:2, and 1:10 orbit-rotation resonances for both rocky and ocean planets (see Table 1). Eccentricity was set to zero in the simulations we present here, with the stellar flux set to the true value at a = 0.22 and not increased to account for orbital averaging as in our one-dimensional study (Wordsworth et al. 2010b). However, we also performed tests with e = 0.38 and found that eccentricity was not critical to the results. Similarly, we assumed zero obliquity to assess the likelihood of atmospheric collapse in the most severe cases, but sensitivity tests showed that its influence on the climate was second order at the high pressures where the atmosphere remained stable.

CO₂ gas had a powerful greenhouse warming effect in our simulations because GJ581d orbits a red dwarf star. In the solar system, Rayleigh scattering reflects the bluer incoming starlight much more effectively, and hence a planet receiving the same flux as GJ581d would be uninhabitable for any CO₂ pressure. In our rocky simulations, for pressures below ∼10 bar, the atmosphere was indeed unstable and began to condense on the dark side and/or poles of the planet (Figure 2(a)). However, for denser atmospheres, we found that horizontal heat transport and greenhouse warming became effective enough to remove the threat of collapse and allow surface temperatures above the melting point of water (Figure 1). Rotation rate affected the atmospheric stability through both the insolation (the planet has permanent day and night sides in the most extreme 1:1 resonance case) and the horizontal heat transport (faster rotating atmospheres were less efficient at transporting heat poleward).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Even in the stable simulations, some CO₂ usually condensed in the middle atmosphere (Figure 3), leading to CO₂ ice cloud formation as occurs on present-day Mars (Montmessin et al. 2007). CO₂ clouds typically had a net warming effect (up to 12 K), depending on the cloud microphysical assumptions, due to the scattering of infrared radiation from the ground (Forget & Pierrehumbert 1997). To check the robustness of our results, we also performed tests with the (unrealistic) assumption of no CO₂ cloud radiative effects, and found that the atmospheres were colder, but still stable above around 10 bar.

Zoom In Zoom Out Reset image size

In the ocean planet simulations, we neglected all oceanic horizontal heat transport. This assumption led us to overestimate global temperature differences and hence the probability of atmospheric collapse, in keeping with the aim of a conservative habitability estimate. We found that atmospheric H₂O vapor greatly increased warming, while H₂O cloud formation low in the atmosphere tended to cool the planet by increasing the planetary albedo (Figure 3). This led to a transition in the climate as total pressure increased. At 10 bar and below, cooling effects dominated and runaway glaciation occurred, followed by atmospheric collapse. From 20 bar, the positive feedback of water vapor on greenhouse warming raised mean temperatures significantly compared with the rocky case (Figure 2(b)). In the intermediate region, climate stability was difficult to assess as surface temperature trends were extremely small (as low as a few K per 100 orbits in some cases). There, our results were somewhat sensitive to our microphysical assumptions; for example, lowering the precipitation threshold l₀ (see Table 1) resulted in optically thinner H₂O clouds and stable, cooler climates at lower total pressures. Nonetheless, at higher pressures we found stable, hot climates with little global temperature variations even in the tidally locked cases, regardless of the choice of microphysical parameters.

We investigated the effect of starting with a surface covered by H₂O ice. For the dense atmospheres where the climate was stable, this only weakly changed the planetary albedo, and the ice melted for pressures of 20 bar or more. Hence the runaway (H₂O) glaciation that is possible on Earth (Budyko 1969) would be unlikely to occur on an ocean-planet GJ581d. Our results here contrast with those of Spiegel et al. (2009), who used a one-dimensional energy balance model to investigate the effects of obliquity on planetary habitability, because they neglected the dependence of planetary albedo on the total pressure in dense atmospheres.

We also investigated the possibility that an ice-covered, tidally locked GJ581d with permanent day and night sides could be locally habitable on the day side due to partial melting. However, in this scenario, the day side only stayed warm when the atmosphere was thin and hence inefficient at transporting heat. As a result, the planet's dark side became cold enough for the collapse of even an N₂ atmosphere. An ice planet with only a thin sublimation-driven H₂O atmosphere could have dayside temperatures above 273 K and continual transport of H₂O to the dark side, but the low atmospheric pressure would preclude liquid water except in extremely limited sub-surface regions. Clearly, dense, stable atmospheres offer better prospects for habitability.

Are CO₂ partial pressures of over 10 bar a realistic possibility for GJ581d? Given the planet's assumed gravity, this corresponds to a CO₂ column of only 4–6 bars on Earth. Even though Venus and Earth are smaller, their total CO₂ inventories scaled by planetary mass are around 10–100 times greater than this. On Venus, it is thought that a large fraction of this inventory is in the atmosphere (Bullock & Grinspoon 1996), while on Earth, the atmospheric partial pressure is regulated over geological timescales by the carbonate-silicate cycle (Walker et al. 1981) via plate tectonics. The geophysics of GJ581d is unknown, but if a similar mechanism were present there, its atmospheric CO₂ would stabilize above the level needed to maintain a liquid water cycle by negative feedback.

CO₂-H₂O-(buffer gas) atmospheres are clearly only a subset of all the possible scenarios for GJ581d. If small quantities of additional greenhouse gases such as CH₄ or SO₂ were present, they would increase warming further. However, if these gases were chemically unstable, they would need to be continually emitted by, e.g., volcanic or biological sources to have a long-term effect on the climate. More drastically, the planet could have a thick H₂–He envelope like Uranus or Neptune, or its atmosphere could have been removed entirely by increased stellar activity early in the life of the host star. Both scenarios would exclude surface liquid water (the former through an excess of surface pressure). GJ581d's large mass means that it can retain an atmosphere more easily than Earth or Venus, but the increased extreme ultraviolet (XUV) and ion fluxes from M-class stars early in their history means that removal processes can be much more efficient than in the solar system (Tanigawa & Ikoma 2007; Tian 2009). The escape of atmospheric hydrogen depends sensitively on both XUV flux levels and radiative cooling by H⁺₃ ions, both of which are poorly constrained for GJ581d (Selsis et al. 2007; Koskinen et al. 2007). Moreover, stripping by the stellar wind, which does not discriminate between light and heavy atoms, could have removed large amounts of any gas (including CO₂; Lammer et al. 2007). Detailed modeling of atmospheric escape may be able to constrain the effects of these processes on the atmosphere further.

Unlike the majority of the Kepler planetary candidates, Gliese 581d is relatively close to Earth, so in the future it will be possible to establish which scenario applies to it through direct spectroscopic observations. To determine the sensitivity required for this, we used the outgoing longwave radiation from our simulations to produce synthetic emission spectra for the habitable climates just discussed, along with those for a planet with no atmosphere (Figure 4). We did not calculate emission spectra for the H₂–He dominated case, but equilibrium chemistry calculations⁴ allowed us to establish that for an atmosphere where H₂ is the dominant species by mass, the abundance of CO₂ should be below the limit of detectability. For a rocky/icy planet where volatiles like CO₂ have collapsed on the surface and the atmosphere is thin or non-existent, the flux variations over one orbit will be large (gray region in Figure 4). Observations of CO₂ and H₂O absorption bands with low phase variations (red/blue regions in Figure 4) would hence be a strong indicator of the kind of stable, dense habitable atmospheres discussed in this Letter. From Figure 4, the flux sensitivity necessary for this in the infrared will be of order 10⁻²¹ W m⁻² μm⁻¹, corresponding to a planet/star ratio lower than 10⁻⁶ for wavelengths <20 μm. This is beyond the capabilities of current space and ground-based observatories, but it will become possible with future improvements in instrumental technology.

Zoom In Zoom Out Reset image size

F.S. acknowledges support from the European Research Council (Starting Grant 209622: E₃ARTHs).