© ESO, 2010

1. Introduction

In 2007, radial velocity measurements were used to discover two new planets in the Gl 581 system (Udry et al. 2007). These planets have captured much attention both in the community and among the general public, as their minimum masses were measured to be below 10 M_Earth, and they are close to the edges of their system’s nominal “habitable zone”, i.e., the loosely defined orbital region in which planets can sustain liquid water on their surfaces. The first planet, Gl 581c, which is closer to its star and was the first discovered, was initially estimated to be potentially habitable based on its equilibrium temperature T_eq = 320 K, using an Earth-like planetary albedo 0.29. In contrast, the second planet Gl 581d has an equilibrium temperature T_eq = 195 K for an albedo of 0.2, which suggests it may be too cold to sustain surface liquid water. However, these analyses neglect any possible warming of the surface due to the planet’s atmosphere.

In the first detailed assessment of the potential habitability of these planets, Selsis et al. (2007) reviewed a variety of factors that could influence their climates. They concluded that based on standard assumptions of atmospheric warming by a mixture of CO₂ and H₂O (with possible regulation of CO₂ via the carbonate-silicate cycle), Gl 581c was unlikely to be habitable, while for Gl 581d, the situation was much less clear. According to Kasting et al. (1993), the outer edge of the habitability zone is most likely the distance at which CO₂ condensation begins to occur on the surface of the planet. However, CO₂ condensation in the atmosphere leads to the formation of CO₂ clouds, which can cause a strong warming effect due to the scattering of infrared radiation (Forget & Pierrehumbert 1997). Hence it was concluded that further climate simulations were required.

To investigate the possible climate of Gl 581d under a range of conditions, we have performed one-dimensional radiative-convective calculations. In Sect. 2 we discuss the model we used, while in Sect. 3 we present our results for varying atmospheric compositions, surface albedos, gravity and cloud coverage. We also present some simple three-dimensional simulations that highlight the limitations of the one-dimensional globally averaged approach. In Sect. 4, we discuss the implications of our results and suggest directions for future research.

2. Method

Our radiative scheme is based on the correlated-k model, with the absorption data calculated directly from high resolution spectra. We used the program kspectrum¹ to produce line-by-line spectra from the HITRAN 2008 database (Rothman et al. 2009). These were then converted to correlated-k coefficients for use in the radiative transfer calculations. The CO₂ collision-induced absorption (CIA) was parameterized using the results of Baranov et al. (2004) and Gruszka & Borysow (1998), while the sub-lorentzian profiles of Perrin & Hartmann (1989) were used for the CO₂ far line absorption. For further details and justification of the method used, refer to Wordsworth et al. (2010). Note that at high CO₂ partial pressures (tens of bars), additional sources of CIA may appear that we have not been able to take into account (Eymet et al. 2009). However, these would tend to increase the surface temperature in our calculations. As we are interested in making a conservative estimate of habitability, we can safely neglect them here.

A matrix of coefficients was produced on a 6 × 9 temperature and log-pressure grid T =  { 100,150,...,350 }  K, p =  { 10^-3,10^-2,...,10⁵ }  mbar. For mixed CO₂-H₂O atmospheres, we produced correlated-k coefficients for H₂O volume mixing ratios q_H2O =  { 10^-7,10^-6,...,10^-1 } . We used 32 spectral bands in the longwave and 36 in the shortwave. Sixteen points were used for the g-space integration, where g is the cumulated distribution function of the absorption data for each band. For the stellar spectra, we used the Virtual Planet Laboratory AD Leo and solar data (Segura et al. 2005) for M-class and G-class stars, respectively. The solar zenith angle θ_z was set to 60°, and the stellar flux was further divided by two to give a spatially and temporally averaged value . As planetary albedo is a nonlinear function of zenith angle, the values obtained in 1D can differ slightly from those produced by a three-dimensional model. However, a 60° zenith angle generally overestimates the globally averaged value, and hence predicts lower surface temperatures. We therefore neglected the error due to this difference here. Here S₀ = 381.4 W m^-2 is the unaveraged stellar flux at the Gl 581d semi-major axis a, and e is the eccentricity. F_m is given along with a, e and the stellar luminosity L in Table 1. The factor involving e in the definition of F_m comes from the fact that the flux averaged over one orbit is higher in an eccentric orbit than in an equivalent circular one.

To account for the radiative effects of both clouds and Rayleigh scattering in our simulations, we used the Toon et al. (1989) scheme. Rayleigh scattering was included by the method described in Hansen & Travis (1974), while the properties of the clouds were computed from Mie theory. Refractive indices for the Mie calculation were taken from Hansen (2005) and Warren (1984) for the CO₂ and H₂O clouds, respectively. In the simulations with CO₂ clouds, a simple microphysical scheme was used to calculate the cloud opacity in each layer. Following Forget et al. (2004), the number of condensable nuclei per kilogram of air N_c was treated as a free parameter. At each timestep, the mean cloud particle radius at each level was calculated from N_c and the amount of condensed CO₂. This was then used to calculate the mean particle sedimentation rates and cloud radiative properties.

For water clouds, we used a simple approach. Fixed layers of particles of radius 10 μm and varying optical depth (defined at a reference wavelength 1.5 μm) were placed in the low troposphere, which we defined arbitrarily as the level at which p = 0.5p_s. In simulations with vertically varying amounts of water vapour, we used profiles of the form

(1)

where q_H2O is the water vapour mixing ratio, q_sat the water vapour saturation ratio, and RH the relative humidity. The latter was defined as in Manabe & Wetherald (1967)

(2)

until p/p_surf < 0.02, after which we set q_H2O = 0.

To model convection in the lower atmosphere we used adjustment to keep the lapse rate Γ stable. In dry ideal atmospheres, Γ =  − g/c_p. As dense CO₂ is a non-ideal gas, Γ can deviate from this value in the lower troposphere when the pressure is sufficiently high. We assess the importance of this deviation in the Appendix. The effects of moist convection on the lapse rate were not included. While moist lapse rates are shallower than dry ones, which results in lower surface temperatures, the quantity of unsaturated water vapour in cold atmospheres is low. The error due to this approximation was small in our calculations for surface temperatures below 273 K, and hence unimportant to our assessment of the minimum CO₂ partial pressure necessary for habitability. In contrast, the inclusion of CO₂ condensation above the troposphere was extremely important. We took this effect into account using the algorithm described in Forget et al. (1998) and vapour pressure curves derived from Lide (2000).

Thirty vertical levels in standard sigma coordinates were used, and all simulations were iterated in time until a steady state had been reached (this took up to 50 Earth years of simulation time for the highest pressure runs). In all cases the energy balance (incoming vs. outgoing radiation) of the model was checked at the end of each run. A list of all parameters used in the simulations is given in Table 1.

3. Results

We have studied the effects of a range of climate parameters on the mean surface temperature of Gl 581d. For clarity, this section is subdivided according to the effect studied. When not explicitly stated otherwise in the text, the parameters used for each simulation are those in bold in Table 1.

3.1. Stellar spectrum

To better understand the differences in climate caused by the fact that Gliese 581 is an M-class star, we first performed simulations comparing G-class (Sol) and M-class (AD Leo) stellar spectra. In both cases we normalized the total fluxes to the same value F_m (see Table 1). Figure 1 shows the equilibrium mean surface temperatures obtained as a function of surface pressure for these two cases.

Fig. 1 Surface temperature vs. surface pressure for clear pure CO2 atmospheres under G-class (Sol) and M-class (AD Leo) stellar insolation. Dashed and dotted lines show CO2 and H2O phase transitions, respectively. In the G-class case, the line is discontinued because the atmosphere condenses on the surface for pressures greater than a few bar.

As can be seen, the clear pure CO₂ atmospheres under a G-class star collapse on the surface for pressures greater than about 3 bar, but when the star is M-class, temperatures continue to increase, reaching the water melting point at just over 10 bar and a maximum value at around 30 bar. The essential reason for this difference is that Rayleigh scattering, which has an optical depth τ_R ∝ λ^-4, has a much weaker effect on the red-shifted M-class stellar spectrum.

In addition to reducing the effect of Rayleigh scattering, the red shift in the stellar spectrum leads to increased warming in the upper stratosphere. This effect is clear from the vertical temperature profiles, which are plotted in Fig. 2. As can be seen, near-IR CO₂ absorption leads to a temperature inversion near 0.05 bar, which becomes even more pronounced if CO₂ condensation is neglected in the calculation (Fig. 2b).

Fig. 2 Steady-state vertical temperature profiles for clear pure CO2 atmospheres of different surface pressures under M-class stellar insolation a) with and b) without CO2 condensation included. The dashed line shows the CO2 saturation pressure.

To show the importance of these separate effects more clearly, in Fig. 3 we have plotted the incident stellar spectra and the radiation reaching the ground for cases with and without Rayleigh scattering included. As can be seen, Rayleigh scattering dominates for wavelengths  < 1 μm, while gaseous absorption dominates above 1 μm.

Fig. 3 Incident flux in the visible at the top of the atmosphere (black) and at the ground, with (blue) and without (red) Rayleigh scattering included. In this example the atmosphere consisted of 40 bars pure CO2 (no clouds), while the stellar spectrum was G-class (Sol) for a) and M-class (AD Leo) for b).

As AD Leo has a mean temperature of around 3400 K, the planetary albedo values we obtain are lower than those for the blackbody 3700 K case that was used to define the habitable boundary in Kasting et al. (1993). This leads to a significant increase in the calculated surface temperature. Note that the temperature of Gl 581 is approximately 3200 K, so the planetary albedos we calculate using the AD Leo spectrum should lead to conservative habitability estimates.

3.2. CO₂ clouds

If CO₂ condensation occurs in the upper atmosphere of a planet, it should cause CO₂ cloud formation. This effect is observed for example in the present-day Martian mesosphere (Montmessin et al. 2007). While there are still many unknowns associated with the microphysics of CO₂ clouds, a basic estimate of their effect can be gained through the scheme described in Sect. 2. In Fig. 4a, the surface temperature vs. surface pressure is plotted for simulations with radiatively active CO₂ clouds and varying values of the condensation parameter N_c. The infrared scattering effect described by Forget & Pierrehumbert (1997) is expected to be slightly less efficient when the star is M-class, because a greater portion of the incident stellar radiation is also scattered back to space. Nonetheless, the clouds still increase the surface warming by an amount that increases with N_c, up to a theoretical maximum of around 30 K at 40 bar for N_c = 10⁵ kg^-1. Figure 4b shows that at high pressures, the planetary albedo begins to increase significantly due to CO₂ cloud coverage. This ultimately limits the maximum warming effect that the clouds can provide, as the downward scattering of thermal radiation from the planet’s surface becomes less important than the direct reflection of stellar radiation to space.

3.3. Surface gravity

Sotin et al. (2007) proposed a relationship between planetary mass and radius (r/r_E) = α(M/M_E)^β, with α and β equal to 1.0 and 0.274 for rocky planets and 1.262 and 0.275 for ocean planets, respectively. Given this relation and the current uncertainties in Gl 581d’s mass, we can expect the planet’s surface gravity to be in the range 10 − 30 m s^-2. In Fig. 5 we have plotted surface temperature vs. pressure for simulations with g = 10, 20 and 30 m s^-2. For a given atmospheric pressure, the CO₂ column amount p_s/g (and therefore the total mass of the atmosphere) decreases with g. Hence the primary effect of increasing g, predictably, is to cool the surface.

Fig. 4 a) Surface temperature vs. surface pressure for pure CO2 atmospheres under M-class insolation with CO2 cloud radiative effects included. Dashed and dotted lines show CO2 and H2O phase transitions, respectively. b) Planetary albedo vs. surface pressure for the same simulations as in a).

Fig. 5 Surface temperature vs. surface pressure for clear pure CO2 atmospheres under M-class insolation with varying surface gravity. Surface temperatures are plotted vs. a) pressure and b) CO2 column density. Dashed and dotted lines in a) show CO2 and H2O phase transitions, respectively.

Variations in g also change the adiabatic lapse rate, while leaving the CO₂ saturation pressure unaffected. Hence for a given CO₂ column amount the temperature profile varies with the gravity, which also influences the climate. These changes are important in determining the point at which maximum greenhouse warming occurs (see e.g. the g = 30 ms^-2 case in Fig. 5a), but they do not strongly affect the surface temperature for a given CO₂ column amount at lower pressures. This can be seen in Fig. 5b, where the same three curves are plotted vs. p_s/g. The atmosphere begins to collapse before the surface temperature reaches the water melting point only in the most conservative g = 30 m s^-2 case.

3.4. Surface albedo

Surface albedo has a dramatic effect on the surface temperatures of planets with thin atmospheres, but its importance might be expected to be less in dense or cloudy atmospheres. We tested the effect of variations in surface albedo in our calculations for A_s = 0.1, 0.2 and 0.3. In Figs. 6a and b, we have plotted the surface temperature vs. pressure for clear-sky and N_c = 10⁵ kg^-1 CO₂ cloud cases, respectively.

Fig. 6 Surface temperature vs. surface pressure for pure CO2 atmospheres a) without clouds and b) with CO2 clouds, Nc = 105 kg-1, under M-class insolation with varying surface albedo. Dashed and dotted lines show CO2 and H2O phase transitions, respectively.

As can be seen, the difference in climate between A_s = 0.1 and 0.3 is large in the clear-sky case. While the lower surface albedo leads to T_s > 273 K at 10 bar with a maximum T_s = 315 K at 40 bar, the higher surface albedo causes the atmosphere to collapse on the surface for p_s > 10 bar. However, in the cases with CO₂ cloud cover (Fig. 6b), the differences in surface temperature are generally smaller, and for all albedo values, surface temperatures above the water melting point are achieved before the surface pressure reaches 10 bar. Hence if the atmosphere contains enough CO₂ to increase surface temperatures above 273 K and clouds of some kind are present, the effect of surface albedo on the climate should not be critical.

3.5. N₂ partial pressure

N₂ is an major constituent of the atmospheres of several rocky bodies in the Solar System. We tested the effects it would have on the climate of Gl 581d by repeating the basic calculation of Sect. 3.1 with added partial pressures of 1, 5 and 20 bars of N₂. As there are no important N₂ absorption bands in the visible or infrared, its primary effect on the climate is through a) additional Rayleigh scattering and b) changes in the specific heat capacity and mean molar mass of the atmosphere. We included both these effects in our calculations. N₂ also has a warming effect in the Earth’s atmosphere due to the broadening of CO₂ and H₂O spectral lines (Goldblatt et al. 2009). However, this effect is relatively small (~5 K for Earth conditions) and we are interested in conservative estimates of warming. Hence we neglected it here.

Fig. 7 Surface temperature vs. surface pressure for clear mixed CO2-N2 atmospheres under M-class insolation. Dashed and dotted lines show CO2 and H2O phase transitions, respectively. b) Planetary albedo vs. surface pressure for the same simulations as in a). In these figures, the x-axis shows the CO2 partial pressure.

As can be seen from Fig. 7, the effect of N₂ on the surface temperature is relatively small. Only for p_N2 = 20 bar is the planetary albedo significantly increased above the value for pure CO₂ due to increased Rayleigh scattering. In this extreme case, the temperature fails to reach the water melting point, with the atmosphere collapsing after p_s = 20 bar. However, if N₂ warming were included, this difference would presumably be less significant.

3.6. H₂O vapour and clouds

Water plays an important and complex role in the radiative budget of Earth, through both direct vapour absorption and cloud absorption/scattering. Here we have studied its effects via the simple approach described in Sect. 2. Figure 8 shows runs performed with the relative humidity RH = 0 at all levels, RH = 1 at all levels, and RH defined as in (2). As can be seen, increasing water vapour has the predictable effect of increasing the total warming, by an amount that increases with the surface temperature. This shifts the surface pressure at which we reach the water melting point down to 5 bar. Note that we neglect the effects of moist convection on the tropospheric lapse rate, so the CO₂+H₂O results in Fig. 8 at the highest pressures likely overestimate the temperature somewhat.

Fig. 8 Surface temperature vs. surface pressure for clear mixed CO2-H2O atmospheres under M-class insolation. Dashed and dotted lines show CO2 and H2O phase transitions, respectively.

Figure 9 shows simulations performed with fixed H₂O cloud layers in the troposphere (see Sect. 2) and varying vapor RH humidity as defined in (2). The planetary albedo increases with the cloud optical depth, but the amount by which it does so is lower at high pressures, because then atmospheric absorption/scattering above the clouds is more important. Counterintuitively, however, the clouds increase the surface temperature at low pressures and decrease it at high pressures. This effect can be explained by the fact that the atmosphere is relatively clear in some infrared spectral regions at low pressures. While an H₂O cloud layer reflects starlight, increasing the planetary albedo, it will also absorb outgoing infrared radiation from the surface, effectively acting as a “thermal blanket” and increasing the surface temperature. At high pressures, however, the lower atmosphere absorbs strongly in almost all infrared bands due to increased water vapour and the effects of CO₂ CIA. Hence the clouds have a minimal effect on the outgoing infrared radiation, but still increase the planetary albedo, which leads to a net cooling.

Fig. 9 a) Surface temperature vs. surface pressure for mixed CO2-H2O atmospheres under M-class insolation with no CO2 clouds and fixed H2O cloud coverage. Dashed and dotted lines show CO2 and H2O phase transitions, respectively. b) Planetary albedo vs. surface pressure for the same simulations as in a). In all cases optical depth is defined at the reference wavelength (1.5 μm).

For higher values of H₂O cloud optical depth, we found that at intermediate pressures, the atmosphere collapsed on the surface (see also Fig. 10). This suggests that sufficient H₂O cloud cover at low altitudes may limit the habitability of Gl 581d. However, we used 100% fractional cloud coverage in these simulations, while on Earth, global low-level H₂O cloud coverage is generally less than 50%, with reported optical depths of 4.7 at 0.6 μm (Kitzmann et al. 2010, derived from results of Rossow & Schiffer 1999; and Warren et al. 2007). Hence our results for large optical depths may overestimate the cooling effect of the clouds at high pressures. For a more comprehensive study of the effects of H₂O clouds than we have given here, three-dimensional climate simulations including cloud microphysics are necessary.

3.7. Most conservative scenario

What happens if we combine the most conservative estimates for surface gravity, surface albedo, N₂ abundance and CO₂ cloud cover? Fig. 10 shows a simulation that was performed with g = 30 m s^-2, A_s = 0.3, p_N2 = 20 bar, minimal CO₂ cloud coverage (N_c = 10³ kg^-1) and water vapour present according to (2). The H₂O cloud opacity was varied between τ = 1 and 3. As can be seen, in all three cases the atmosphere collapses at CO₂ pressures between 3 and 20 bar. Nonetheless, it can be seen that at higher pressures, the atmosphere is once again stable, with surface temperatures well above the water melting point at p_CO2 =  40 bar for these values of the H₂O cloud opacity. Given the initial assumption of CO₂ partial pressure of tens of bars, therefore, the outlook for the habitability of Gl 581d would appear rather good.

Fig. 10 Surface temperature vs. surface pressure for the “most conservative case” described in the text. Dashed and dotted lines show CO2 and H2O phase transitions, respectively, while crosses, circles and stars indicate surface temperatures for atmospheres with H2O clouds of opacity τ = 1, 2 and 3, respectively. In this figure, the x-axis shows the CO2 partial pressure.

3.8. Three-dimensional equilibrium temperatures: airless planet case

While one-dimensional simulations of the type described here are useful because of their simplicity, they make the inevitable assumption of globally averaged stellar insolation. As discussed in Selsis et al. (2007), Gl 581d is likely to be in a tidally resonant orbit, which means that there will be dramatic differences in insolation across the surface of the planet. To get an idea of the effect this could have on the climate, we performed some simple equilibrium temperature calculations in three dimensions.

In Fig. 11, we have plotted the temporal evolution of surface temperature on the 0° and 180° longitude sides of an airless planet with the same orbit as Gl 581d, zero obliquity, tidal resonances of 1:1 (left) and 2:1 (right), surface albedo A_s = 0.2 and initial uniform surface temperature 250 K. Note that the “dark side” temperature in the plot on the left (solid line) is still decreasing toward 0 K in the simulation due to the finite thermal inertia that was given to the surface.

Fig. 11 Snapshot surface temperature on the points of a hypothetical tidally resonant Gliese 581d with minimal atmosphere and surface albedo A = 0.2 that are sub- (dashed line) and anti-stellar (solid line) at periastron. Upper and lower dotted lines indicate the equilibrium sub-stellar temperature at periastron and apastron, respectively. The initial surface temperature in the simulations was 250 K at time = 0 days.

The dotted lines in the figure correspond to the equilibrium sub-stellar temperature at periastron

(3)

and apastron

(4)

As can be seen, at one location the local surface temperature can vary by as much as 300 K over the course of one orbit. This variability reduces for planets with atmospheres as a function of the surface pressure – on Venus, for example, where the diurnal period is 171 Earth days and the pressure is 93 bars, the only significant surface temperature variation occurs due to topographic effects (Mueller et al. 2008). However, the decrease in variability of T_s with pressure will be dependent in general on the specific atmospheric and orbital parameters of each planet. Hence further modeling is required before global temperature variations can be estimated for the hypothetical Gl 581d atmospheres we have studied here.

The second point of interest raised by this simple 3D model is that locally, the temperature on the planet rises well above the water melting point during the “warm season” (substellar point at periastron), despite the fact that there is no atmospheric warming. This raises an interesting question: could a mainly frozen planet in this orbit with a thin (p_s > 6.1 mbar) atmosphere contain a zone where liquid water forms seasonally? Two extreme cases can be envisaged, depending on the total amount of water present. If the fraction of water on the planet is small, the regions of low solar insolation are likely to act as cold traps in which water ice forms, leaving the zone in which liquid water is possible dry. If the planet’s surface is covered by H₂O, however, there will be no limit to the amount of local evaporation that can occur. In this case, the local climate could potentially be quite interesting, although the high albedo of water ice would reduce the surface temperatures seen in Fig. 11.

4. Discussion

Through one-dimensional radiative-convective climate modeling, we have found that as little as 5 bars of CO₂ (~2.5 bar equivalent column amount under Earth gravity) may be sufficient to maintain a global mean temperature above the melting point of water. For CO₂ partial pressures greater than 40 bar (~13 bar equivalent column amount under Earth gravity), we found that the mean surface temperature of Gliese 581d is high enough to allow surface liquid water even when the most conservative values of gravity, surface albedo and N₂ partial pressure are chosen, although strong dependence on the total H₂O cloud opacity was observed.

We tested the dependence of surface temperature on a wide range of atmospheric and geophysical parameters. Because the M-class spectrum of Gliese 581 is red-shifted with respect to that of the Sun, heating due to near-IR absorption in the upper atmosphere increases in our simulations, while Rayleigh scattering is much less important. This is in agreement with the general habitability study of Kasting et al. (1993), but we find the difference to be even stronger than was reported there, due to the low effective temperature of Gl 581. We also estimate an modest warming effect due to CO₂ clouds (up to 30 K at p_CO2 = 50 bar), dependent on details of the cloud microphysics, although the efficiency of this process is reduced compared to the case for G-class stars.

We find that increasing gravity primarily decreases the CO₂ column amount at which the maximum surface temperature occurs. Unsurprisingly, increased surface albedo decreases the mean surface temperature, although if there is cloud coverage, this effect is relatively small for the high CO₂ pressures at which the surface starts to become habitable. Increasing the partial pressure of N₂ has little effect below partial pressures p_N2 of tens of bars. Finally, water vapour increases the surface temperature by a moderate amount below T_s = 273 K, after which its role becomes increasingly important. We briefly assessed the impact of H₂O clouds, and found that they warmed (cooled) the surface at low (high) pressures, due to the interplay between their effects on incoming stellar and outgoing thermal radiation.

As discussed by Selsis et al. (2007), N₂-CO₂-H₂O atmospheres likely only represent a fraction of the possible range for Earth-like/super-Earth extrasolar planets. While the presence of radiatively inactive gases such as argon would be of little importance to the calculations discussed here, small amounts of other greenhouse gases such as methane would result in increased absorption, if such gases could remain chemically stable in the atmosphere. Unless the near-IR absorption of these gases were large enough to cause a strong temperature inversion in the upper atmosphere, this would further increase the surface temperature of the planet. The presence of other clouds or haze layers could alter the radiative balance in less predictable ways, depending on their scattering/absorption properties and height in the atmosphere. Unfortunately, effects of this kind are difficult to constrain further until direct atmospheric observations become available.

The results discussed here clearly have interesting implications, as they show that Gl 581d could be the first discovered habitable exoplanet. They are also testable by future observations – dense CO₂ atmospheres have recognizable spectral signatures that could be detected, for example, using the proposed Darwin or TPF missions. However, there is an important limitation to this study: all simulations performed were one-dimensional. As we noted in Sect. 3.8, calculations of the globally averaged surface temperature neglect variations due to local changes in stellar insolation (as well as topography and other effects). This has serious implications for dense, relatively cold CO₂ atmospheres, as if the dark side temperature of the planet is too low, CO₂ will condense on the surface. To test whether a dense atmosphere on Gl 518d would be stable against the collapse of CO₂ in this way, we plan to repeat the calculations reported here in the future using a three-dimensional climate model.

Appendix A

To assess the importance of the non-ideal behaviour of CO₂ gas at high pressures in our calculations, we used the modified Redlich-Kwong equation of state proposed by Li et al. (2006). The Redlich-Kwong equation takes the form

(5)

where p, T, v and R are pressure, temperature, specific volume and the ideal gas constant, respectively. Various definitions of a and b are possible; in Li et al. (2006)a = a₁(T) + a₂(T)p, with a₁ and a₂ polynomial functions of temperature, while b = 0.08664RT_c/p_c, with T_c and p_c the critical point temperature and pressure for CO₂. Equation (5) can be written as a cubic equation in terms of v and hence is solvable analytically. The inverse compressibility factor β = Z^-1 = RT/pv calculated from (5) agrees approximately with the values stated in Kasting (1991) for the triple point and saturated CO₂ at 0 °C. This is unsurprising, as the equation used there was of the same basic form.

To get a first-order estimate of the effect of non-ideality on the tropospheric lapse rate (and hence on the climate), we approximated (5) to a virial form as

(6)

with

(7)

Given the definition of specific heat capacity in Kasting (1991)

(8)

it is clear that in the virial approximation, to leading order

(9)

Various formulations for c_p(0,T) were considered, including the quadratic function given in Kasting (1991) and the power law function of Lebonnois et al. (2010). We found that the latter expression, which is defined as , had the greatest accuracy vs. reference data (Lide 2000). Using this and calculating B″(T) numerically by finite differences, we found for T = 250 K and p =  0, 50 bar that c_p = 795.9 and 825.0 J K^-1 kg^-1, respectively. Climate calculations performed with fixed specific heat capacities of 750 and 850 J K^-1 kg^-1 showed maximum temperature differences of only 2.2 K over

the entire range 0.2–50 bar, with the latter value of c_p giving colder results. Hence we decided to ignore variations in c_p with altitude in our calculations, and used the upper limit c_p = 850 J K^-1 kg^-1 everywhere in the pure CO₂ calculations. For mixed CO₂-N₂ atmospheres, a weighted mean value for c_p was used, with the heat capacity of N₂ taken to be a constant 1003 J K^-1 kg^-1.

Acknowledgments

R. Wordsworth acknowledges funding from l’Institut de France, Fondation Cino et Simone Del Duca. F. Selsis acknowledges support from the European Research Council (Starting Grant 209622: E₃ARTHs).

References

All Tables

All Figures

	Fig. 1 Surface temperature vs. surface pressure for clear pure CO2 atmospheres under G-class (Sol) and M-class (AD Leo) stellar insolation. Dashed and dotted lines show CO2 and H2O phase transitions, respectively. In the G-class case, the line is discontinued because the atmosphere condenses on the surface for pressures greater than a few bar.
In the text

	Fig. 2 Steady-state vertical temperature profiles for clear pure CO2 atmospheres of different surface pressures under M-class stellar insolation a) with and b) without CO2 condensation included. The dashed line shows the CO2 saturation pressure.
In the text

	Fig. 3 Incident flux in the visible at the top of the atmosphere (black) and at the ground, with (blue) and without (red) Rayleigh scattering included. In this example the atmosphere consisted of 40 bars pure CO2 (no clouds), while the stellar spectrum was G-class (Sol) for a) and M-class (AD Leo) for b).
In the text

	Fig. 4 a) Surface temperature vs. surface pressure for pure CO2 atmospheres under M-class insolation with CO2 cloud radiative effects included. Dashed and dotted lines show CO2 and H2O phase transitions, respectively. b) Planetary albedo vs. surface pressure for the same simulations as in a).
In the text

	Fig. 5 Surface temperature vs. surface pressure for clear pure CO2 atmospheres under M-class insolation with varying surface gravity. Surface temperatures are plotted vs. a) pressure and b) CO2 column density. Dashed and dotted lines in a) show CO2 and H2O phase transitions, respectively.
In the text

	Fig. 6 Surface temperature vs. surface pressure for pure CO2 atmospheres a) without clouds and b) with CO2 clouds, Nc = 105 kg-1, under M-class insolation with varying surface albedo. Dashed and dotted lines show CO2 and H2O phase transitions, respectively.
In the text

	Fig. 7 Surface temperature vs. surface pressure for clear mixed CO2-N2 atmospheres under M-class insolation. Dashed and dotted lines show CO2 and H2O phase transitions, respectively. b) Planetary albedo vs. surface pressure for the same simulations as in a). In these figures, the x-axis shows the CO2 partial pressure.
In the text

	Fig. 8 Surface temperature vs. surface pressure for clear mixed CO2-H2O atmospheres under M-class insolation. Dashed and dotted lines show CO2 and H2O phase transitions, respectively.
In the text

	Fig. 9 a) Surface temperature vs. surface pressure for mixed CO2-H2O atmospheres under M-class insolation with no CO2 clouds and fixed H2O cloud coverage. Dashed and dotted lines show CO2 and H2O phase transitions, respectively. b) Planetary albedo vs. surface pressure for the same simulations as in a). In all cases optical depth is defined at the reference wavelength (1.5 μm).
In the text

	Fig. 10 Surface temperature vs. surface pressure for the “most conservative case” described in the text. Dashed and dotted lines show CO2 and H2O phase transitions, respectively, while crosses, circles and stars indicate surface temperatures for atmospheres with H2O clouds of opacity τ = 1, 2 and 3, respectively. In this figure, the x-axis shows the CO2 partial pressure.
In the text

Fig. 11 Snapshot surface temperature on the points of a hypothetical tidally resonant Gliese 581d with minimal atmosphere and surface albedo A = 0.2 that are sub- (dashed line) and anti-stellar (solid line) at periastron. Upper and lower dotted lines indicate the equilibrium sub-stellar temperature at periastron and apastron, respectively. The initial surface temperature in the simulations was 250 K at time = 0 days.

In the text