THE HABITABLE ZONE OF EARTH-LIKE PLANETS WITH DIFFERENT LEVELS OF ATMOSPHERIC PRESSURE

The degree of realism of the EBM is largely determined by the formalism adopted to describe the physical quantities C, D, I, S, and A. The recipes that we adopt are inspired by previous EBM work published by SMS08, WK97, and North et al. (1983). In our EBM, we have introduced several new features aimed at addressing some limitations of the method and broadening the range of application of the model.

An intrinsic limitation of diffusive EBMs is that the physics that drives the horizontal transport is more complex than implied by the diffusion term that appears in Equation (1). We have introduced a dependence of the diffusion coefficient D on the stellar zenith distance. In this way, the efficiency of horizontal heat transport becomes a function of latitude and orbital phase, as in the case of real planets. In a real planetary atmosphere, heat transport across latitudes is governed by large-scale fluid motions and their instabilities (such as baroclinic instability). Here, such effects are parameterized in an overly simplified way by assuming diffusive transport, albeit with a variable diffusion coefficient.

Previous EBM work has been applied to planets with Earth atmospheric pressure. Given the impact of atmospheric pressure on the climate and habitability, we have generalized the model for applications to terrestrial planets with different levels of surface pressure. With this aim, in addition to the pressure dependence of the diffusion coefficient, D, already considered by WK97, we have also taken into account the pressure dependence of the OLR, I, and of the effective thermal capacity, C. Our model allows permanent ice covers to be formed and calculates the cloud coverage according to the type of underlying surface. All the above features are presented in the Appendix, where we provide a full description of the method.

The diffusion equation (1) is a partial differential equation in two variables, x and t, that can be solved by separating the terms containing the temporal and spatial derivatives. The temporal part can be written in the form ∂T/∂t = f(t, T[x]), which can be solved with the Runge–Kutta method. In practice, we use the routine odeint (Press et al. 1992) characterized by an adaptive step size control. To solve the spatial part, f(t, T[x]), we use the Euler method. To this end, we discretize x in N equispaced cells, delimited by (N + 1) cell borders, and use a staggered grid, i.e., the quantities of interest are calculated at the center of the cells, whereas their derivatives are calculated at the borders. Boundary values are obtained assuming that the diffusion term and the partial derivative of T are null at the poles.

To validate the code, we used the same set of prescriptions of the physical quantities adopted by SMS08 and we verified the capability of the code to reproduce the results published in that paper. As an example, we show in Figure 1 the mean annual temperature–latitude profile of the Earth obtained in this way (dotted line), which is the same as shown in Figure 2 of SMS08 (solid line). We also verified that our code recovers the exactly same seasonal evolution of the zonal temperatures (Figure 3 in SMS08), as well as the "snowball" transition predicted to occur when the latitudinal diffusion is decreased by a factor of nine (Figure 4 in SMS08).

Zoom In Zoom Out Reset image size

After the adoption of our own set of physical recipes, described in the Appendix, the model was calibrated to reproduce the Earth climate data. The calibration was performed by adopting astronomical and planetary parameters appropriate for Earth and tuning the remaining parameters in such a way so as to reproduce the surface temperature and albedo of Earth (Figures 1 and 2).

Zoom In Zoom Out Reset image size

The set of Earth data adopted for the calibration is listed in Table 1. The incoming stellar radiation, S = S(t, x), is fully specified by the solar constant and the Earth orbital parameters and obliquity. The zonal values of Earth ocean fraction f_o, not shown in the table, were taken from Table III of WK97.

Table 1. Earth Data

Parameter	Value	Comment
q0	1361.6 W m−2	Solar constanta
a	1.000 AU	Semi-major axis
e	0.01671022	Orbital eccentricity
	23.43929	Obliquity
Am, ○	0.328	Mean annual global albedob
Tm, ○	287.44 K	Mean annual global surface temperaturec
pt, ○	1.0132 × 105 Pa	Total Earth surface pressure
CO2	380 ppmV	Volumetric mixing ratio of CO2
CH4	1.7 ppmV	Volumetric mixing ratio of CH4
cp, ○	1.005 × 103 J kg−1 K−1	Specific heat capacity of the atmosphere
m○	28.97	Mean molecular weight of the atmosphere

Notes. ^aFrom Kopp & Lean (2011). The measurement given in that paper, 1360.8 ± 0.5 W m⁻², was obtained during the 2008 solar minimum. The excursion between solar minimum and maximum quoted in the same paper amounts to 1.6 W m⁻². We have added half this excursion to the value measured at the minimum. ^bArea-weighted mean annual albedo of Earth measured from the average ERBE data for the period 1985–1989 (data taken from courseware of Pierrehumbert 2010). ^cArea-weighted mean annual surface temperature of Earth measured from ERA Interim data for the years 1979–2010 (see Dee et al. 2011).

Download table as:  ASCIITypeset image

The set of fiducial parameters that we adopt for our EBM is shown in Table 2. Some of these parameters are taken from previous work, cited in the last column of the table. Further information on the choice of the parameters is given in the Appendix.

Table 2. Fiducial Model Parameters

Parameter	Fiducial Value	Comment	Equation	Source
Catm, ○	10.1 × 106 J m−2 K−1	Effective thermal capacity of the Earth atmosphere	(A1)	Pierrehumbert (2010)
Cml50	210 × 106 J m−2 K−1	Effective thermal capacity of the oceans	(A2)	WK97
Csolid	1 × 106 J m−2 K−1	Effective thermal capacity of the solid surface	(A2)	This work
D○	0.600 W m−2 K−1	Diffusion coefficient	(A5)	Pierrehumbert (2010)
$\mathcal {R}$	6	Maximum excursion of diffusion efficiency	(A9)	This work
al	0.20	Surface albedo of lands	(A13)	WK97
ail	0.85	Surface albedo of ice on lands	(A13)	Pierrehumbert (2010)
aio	0.62	Surface albedo of ice on ocean	(A13)	This work
fcw	0.67	Cloud coverage on water	(A13)	This work
fcl	0.50	Cloud coverage on land	(A13)	This work
fci	0.50	Cloud coverage on ice	(A13)	This work

Download table as:  ASCIITypeset image

After the calibration procedure, our Earth model yields a mean annual global temperature T_{m, ○} = 287.8 K and a mean annual global albedo A_{m, ○} = 0.320, in fine agreement with the Earth values (Table 1). In the left panel of Figure 1, we compare the mean annual temperature–latitude profile predicted by our model (solid line) with experimental data (crosses). The agreement with the observations is excellent in most latitudes, with an area-weighted rms deviation of 2.2 K. However, the model fails to reproduce adequately the temperature profile at Antarctic latitudes. This is likely due to the peculiarities of Antarctica, such as its high elevation with respect to the sea level, that are not accounted for in EBMs. In any case, compared with previous EBM work, our model yields a better agreement at Antarctic latitudes. The predictions of the model of SMS08 are shown for comparison in Figure 1 (dashed line); the model by WK97 predicts a mean annual temperature of 260 K at the southernmost latitudes (see their Figure 4), well above the experimental data.

In the right panel of Figure 1, we compare the mean annual albedo–latitude profile of the Earth predicted by our EBM (solid line) with the experimental data obtained by ERBE for the period 1985–1989 (crosses). The agreement is reasonable but not perfect. The main reasons for disagreement are (1) the over-simplified physical prescriptions of the EBM, which do not consider the atmospheric scattering of the short-wavelength radiation, and (2) the limited amount of fine-tuning of the albedo parameters (e.g., a constant value of albedo for all land). For the sake of comparison, in the same figure we show the mean annual albedo profile obtained with the Earth model of SMS08 (dashed line), where the albedo is a simple function of temperature. At variance with this previous formulation, our parameters can be varied to model the albedo of planets with different types of surfaces. It is worth mentioning that the top-of-atmosphere recipes for the albedo given by WK97, which include a treatment of atmospheric scattering, cannot be applied to the present work. This is because the radiative–convective calculations of WK97 were performed by varying the concentration of CO₂, while we are interested in varying the total pressure keeping the composition fixed.

Careful inspection of the seasonal dynamics shows that in the northern hemisphere the temporal evolution of the model slightly lags that of the data (see Figure 2). This is probably due to an intrinsic limitation of one-dimensional EBMs, related to the calculation of the zonal thermal capacity: the zonal average of C tends to be overestimated in zones with ocean fractions above ≈0.1 because the effective thermal capacity of the oceans is much higher than that of the lands (WK97). As a consequence, our EBM overestimates the thermal inertia at the mid-northern latitudes of Earth (with ocean fractions f_o = 0.3–0.5) and predicts that the seasonal temperature should vary more slowly than is observed.

For the purpose of this investigation, an important feature of the temperature profile is the location of the ice line, which affects the calculation of planetary habitability. To test this feature, we compared the latitude φ_○ where T = 273 K in our EBM and in the Earth data. We find a mean annual difference 〈φ_{○, model} − φ_{○, data}〉 = +06 ± 32 in the southern hemisphere and −24 ± 104 in the northern hemisphere (average over 12 months). The small offset 〈φ_{○, model} − φ_{○, data}〉 in both hemispheres testifies in favor of the quality of the ice line calibration. The large scatter in the northern hemisphere is due to the above-discussed seasonal time lag resulting from the overestimate of the thermal capacity.

More details on the model calibration are given in the Appendix, where we show that, in addition to the temperature and albedo data, also the OLR data of Earth have been used to constrain our EBM (see Appendix A.3).

To investigate the effects of surface pressure on the location of the HZ, we generate a large number of models of an Earth-like planet by varying a and p while keeping constant the other orbital and planetary parameters. To focus our problem, we consider an atmosphere with an Earth-like composition but variable level of total pressure. In practice, we scale the partial pressures of non-condensable gases by a constant factor. In this way, the mixing ratios of non-condensable greenhouses, such as CO₂ or CH₄, are always equal to that of the Earth atmosphere. The mixing ratio of water vapor is instead set by the temperature and relative humidity. This scaling of the pressure is in line with the assumptions adopted in the calibration of the OLR described in Appendix A.3.

Each simulation starts with an assigned value of temperature in each zone, T_start, and is stopped when the mean global annual temperature, T_m, converges within a prefixed accuracy. In practice, we calculate the increment δT_m every 10 orbits and stop the simulation when |δT_m| < 0.01 K. The convergence is generally achieved in fewer than 100 orbits. Tests performed with widely different values of T_start > 273 K indicate that the simulations converge to the same solution. We adopted a value of T_start = 275 K, just above the threshold for ice formation, so that the simulation starts without ice coverage and without artificial triggering of snowball episodes. The choice of a low value of T_start is dictated by our interest in exploring a wide range of pressures: T_start = 275 K is below the water boiling point even at very low values of pressure, when the boiling point is just a few kelvin above the freezing point.

The starting value of total pressure p in the simulations is the total pressure of non-condensable gases, i.e., the "dry air" pressure. In the course of the simulation, p is updated by adding the current value of the partial pressure of water vapor. This is calculated as

$$\begin{equation} p{\rm H}_2{\rm O} = {\it rh} \times p_{{\rm sat},w}(T), \end{equation} \tag{ 2 }$$

where rh is the relative humidity and p_{sat, w}(T) is the saturated water vapor pressure.⁴ Also the specific heat capacity, c_p, and the mean molecular weight, m, of the atmosphere are recalculated at each orbit taking into account the contribution of water vapor. The thermal capacity and the diffusion coefficient, which depend on c_p and m, are also updated.

In addition to the regular exit based on the convergence criterion, the simulations are stopped when the mean planet temperature exceeds some prefixed limits (T_min, T_max). This forced exit without convergence is introduced to minimize the computing time when we run a large number of simulations, as in our case. We adopt T_min = 220 K, well below the limit of liquid water habitability. The value of T_max is based on the water-loss limit criterion explained below. When T_m is outside the interval (T_min, T_max), the simulation is stopped and the indices of habitability are set to zero.

The existence of a planetary water-loss limit is based on the following arguments. At a given value of relative humidity, pH₂O increases with temperature according to Equation (2) because the saturated pressure of water vapor increases with T. The dominant feedback mechanism of water vapor is the enhancement of the IR opacity, which tends to raise the temperature. In extreme cases, this positive feedback may lead to a complete water vaporization via a runaway greenhouse effect, followed by a loss of hydrogen in the upper atmosphere via EUV stellar radiation. All together, these effects indicate the existence of a temperature limit above which water is lost from the planet. In the case of Earth, the water-loss limit is predicted to occur at ≃ 340 K, while for a planet with p = 5 bar at ≃ 373 K (Selsis et al. 2007 and references therein). These values of temperature approximately lie at 90% of the liquid water temperature range calculated at the corresponding value of pressure. On the basis of these arguments, we adopt as a water-loss limit the value

$$\begin{equation} T_\mathrm{max}= T_{\rm ice}(p) + 0.9 \times [ T_{\rm vapor}(p) - T_{\rm ice}(p)], \end{equation} \tag{ 3 }$$

where T_ice(p) and T_vapor(p) are the melting and boiling points of water at pressure p (Lide 1997). The adoption of the water-loss limit (3) minimizes the difficulty of tracking the effects of water vapor in the atmosphere when the temperature is high (see Appendix A.3). The results on runaway water-loss events are conservative given our choice of adopting a "cold start" in the simulations.

For each set of input parameters, our simulations provide a matrix of surface temperatures calculated at discrete values of latitude and time, T(φ, t). We use different tools to analyze these results. To cast light on specific cases, we investigate the latitude–temperature profile and its seasonal evolution. To have a global view of a large set of results, we build up maps of mean annual global temperature and habitability.

We call mean planet temperature, T_m, the mean global annual temperature at the planet surface. To calculate this quantity, we average T(φ, t) on latitude and time. The average on latitude is weighted in area.

To estimate the surface habitability from the results of our simulations, we require that the temperature and the pressure lie within the range of liquid water in the phase diagram of H₂O. We consider values of total pressure above the triple point and below the critical pressure. Within this interval, we define a habitability function such that

$$\begin{equation} H(\varphi ,t) = \left\lbrace \begin{array}{@{} ll}1 & {\rm if}\; { T_\mathrm{ice}(p) \le T(\varphi ,t) \le T_\mathrm{vapor}(p)} \\ 0 & {\rm otherwise}. \end{array}\right. \end{equation} \tag{ 4 }$$

Several indices of surface habitability can be defined by integrating H(φ, t) in different ways. Following SMS08, we introduce the temporal habitability

$$\begin{equation} f_\mathrm{time}(\varphi) = {\int _{0}^P dt H(\varphi ,t) \over P}, \end{equation} \tag{ 5 }$$

i.e., the fraction of orbital period P in which the planet is habitable at a given latitude, and the regional habitability

$$\begin{equation} f_\mathrm{area}(t)= \frac{\int _{-{\pi /2}}^{+{\pi /2}} \, d\varphi \, [ H(\varphi ,t) \, \cos \varphi ]}{2}, \end{equation} \tag{ 6 }$$

i.e., the fraction of planet surface that is habitable at a given time. By integrating the habitability function both in latitude and time, we obtain the mean global annual habitability

$$\begin{equation} h = \frac{\int _{-{\pi /2}}^{+{\pi /2}} \, d\varphi \, \int _{0}^P \, dt \, [ H(\varphi ,t) \, \cos \varphi ] }{2P}. \end{equation} \tag{ 7 }$$

This represents the mean fraction of planet surface that is habitable during the orbital period.

In addition to these indices, already defined by SMS08, we introduce here the index of continuous habitability

$$\begin{equation} h_c = \frac{\int _{-{\pi /2}}^{+{\pi /2}} \, d\varphi \, f^{\prime }(\varphi) \, \cos \varphi }{2} , \end{equation} \tag{ 8 }$$

with

$$\begin{equation} f^{\prime }(\varphi) = \left\lbrace \begin{array}{@{} ll} 1 & {\rm if}\; {f_\mathrm{time}(\varphi) =1} \\ 0 & {\rm if}\; {f_\mathrm{time} (\varphi) <1}. \end{array}\right. \end{equation} \tag{ 9 }$$

The index h_c represents the fraction of planet surface that is continuously habitable during the orbital period. This index vanishes if all the latitude zones undergo a period of non-habitability in the course of the orbital period. By construction, it is always h_c ⩽ h.

In Figure 3, we show the map of habitability obtained from our simulations for an Earth-like planet in circular orbit around a Sun-like star. To obtain this figure, we run a total of 4032 simulations covering the interval of semi-major axis 0.65 AU ⩽ a ⩽ 1.35 AU, with a step δa = 0.01 AU, and the interval of pressure 10^−2.0 bar ⩽ p ⩽ 10^+0.8 bar, with a constant logarithmic step δlog p (bar) = 0.05. The values of the parameters that are kept constant in these simulations are listed in Table 3.

Zoom In Zoom Out Reset image size

Table 3. Model Parameters Fixed in the Simulations

Parameter	Adopted Value	Comment
N	37	Number of latitude zones
L⋆	L☉	Stellar luminositya
M⋆	M☉	Stellar massb
e	0.00	Orbital eccentricity
ω	0.00	Argument of the pericenter
	23.44	Axis obliquity
g	9.8 m s−2	Surface gravitational accelerationc
Prot	1 day	Rotation period
fo	0.70	Ocean fraction (constant in all latitude zones)

Notes. ^aThe solar luminosity L_☉ is calculated from the adopted value of solar constant, q₀ (Table 1). ^bThe stellar mass, in conjunction with the semi-major axis a, determines the orbital period adopted in each simulation. ^cThe surface gravitational acceleration is used in the radiative calculations of the OLR (Appendix A.3).

Download table as:  ASCIITypeset image

The map of Figure 3 shows the results obtained at each point of the plane (a, p) in terms of mean annual global temperature and habitability. The filled circles in the map indicate the positions in the plane (a, p) where h > 0; the value of total pressure associated to these symbols includes the partial pressure of water vapor updated in the course of the simulation. Crosses indicate positions on the plane where the simulations were forced to exit; the total pressure for these cases is the starting value of dry air pressure (see Section 3.1). Empty areas of the map, as well as crosses, indicate a location of non-habitability.

The simulations yield information not only on the degree of habitability, through the index h, but also on the "quality" of the habitability, through a detailed analysis of the seasonal and latitudinal variations of the temperature, as we shall see in the next section. However, some cautionary remarks must be made before interpreting these data.

The model has been calibrated using Earth climatological data. These data span a range of temperatures, roughly 220 K ≲ T ≲ 310 K, not sufficient to cover the broad diversity expected for exoplanets, even if we just consider those of terrestrial type. Given the fundamental role of temperature in the diffusion equation, one should be careful in using the physical quantities outside this range, where direct calibration is not possible. In this respect, the major reason of concern is the estimate of the OLR that has been done with radiative calculations (see Appendix A.3). The difficulty of calibrating the OLR outside the range of terrestrial temperatures makes uncertain the exact localization of the inner and outer edges of the HZ. In particular, the results with T_m ≳ 330 K should be treated with caution, given the strong effects of water vapor predicted to occur in this temperature range, which are not directly testable. These cases lie in the region of high pressure in Figure 3 (symbols color-coded in orange and red). The fact that in these cases also the pressure is quite different from the Earth value makes these results particularly uncertain. In the following discussion, we will consider these particular results to be purely tentative. We note that the difficulty of making climate predictions outside the parameter space sampled by the Earth is a common problem of any type of climate model, no matter how sophisticated. In this respect, simple models, like our EBM, help to obtain preliminary predictions to be tested by subsequent investigations.

The circumstellar HZ shown in Figure 3 shows several characteristics, in terms of mean planet temperature and habitability, that can be summarized as follows.

The radial extent of the HZ increases with p. The outer edge extends from 1.02 to 1.18 AU when the pressure rises from 0.1 to 3 bar. The inner edge approaches the star from 0.87 to 0.77 AU in the same pressure interval. No habitability is found below p ≃ 15 mbar.

The broadening of the HZ with increasing pressure is accompanied by an increase of the interval of mean planet temperatures spanned at constant p. At high pressures, most of the broadening of the HZ is contributed by the area of the plane where the solutions have mean temperatures T_m ≳ 60°C (i.e., T_m ≳ 333 K; orange and red symbols above the red line in the figure). If we focus on the interval of mean temperatures 0°C ≲ T_m ≲ 60°C (region with 273 K ≲ T_m ≲ 333 K between the magenta and red line), the broadening of the HZ is quite modest.

Remarkable differences exist between the low- and high-pressure regimes. At low pressures (p ≲ 0.3 bar) the habitability undergoes intense variations in the plane (a, p), with a general trend of increasing h with increasing p. At high pressure (p ≳ 1 bar) the habitability is approximately constant and high, with sudden transitions from h ≃ 1 inside the HZ, to h ≃ 0. The mean planet temperature also shows different characteristics between the low- and high-pressure regimes. Starting from pressure p ≳ 0.3 bar, the curves of equal temperature tend to move away from the star as the pressure increases. This behavior is not seen at lower pressures, where the HZ at a given temperature does not significantly change its distance from the star.

Another interesting feature of Figure 3 is the location of the line of constant mean planet temperature T_m = 273 K, indicated as a magenta solid line superimposed on the symbols of habitability. On the basis of the liquid water habitability criterion, one would expect a coincidence of this line with the outer edge of the HZ. This is true at p ≳ 2 bar, but not at lower values of pressure, the mismatch being quite large at the lowest values of pressure considered. The reason is that, using an EBM, we can determine whether some latitudinal zones have temperatures larger than zero, even when the mean planet temperature is lower. When this is the case, the planet is partly habitable.

Further insight on the differences between low and high pressures is offered by Figure 4, where we plot the surface habitability and temperature as a function of p. The three panels of the figure show the results obtained at three constant values of semi-major axis (i.e., insolation). In addition to the habitability h (lower solid curve) and the mean planet temperature T_m (upper solid curve), we show the minimum and maximum planet temperature dot-dashed and dashed curves, respectively). These temperatures are measured regardless of the time of year or of the latitude, i.e., the maximum temperature could be found at low latitudes and in summer, while the minimum, at high latitudes in winter. The general rise of mean temperature and habitability with increasing pressure is clear in each case. Moreover, one can see that the excursion between extreme temperatures is quite large at low pressure (ΔT ∼ 100–200 K), but becomes increasingly smaller as the pressure rises. Planets with high atmospheric pressure have a rather uniform surface temperature.

Zoom In Zoom Out Reset image size

The rich phenomenology described above results from the interaction of distinct physical factors that we now discuss.

The links between surface pressure and habitability are rather intricate. Pressure variations influence both the planet temperature and the extent of the liquid water temperature interval used to define the habitability. To disentangle these effects we first discuss the influence of pressure on the temperature and then on the habitability.

4.2.1. Surface Pressure and Planet Temperature

Variations of surface pressure affect the temperature in two ways. First, for a given atmospheric composition, the infrared optical depth of the atmosphere will increase with pressure. As a result, a rise of p will always lead to a rise of the greenhouse effect and temperature. Second, the horizontal heat transport increases with pressure. In our model, this is reflected by the linear increase with p of the diffusion coefficient D (Equation (A5), Appendix A.2). At variance with the first effect, it is not straightforward to predict how the temperature will react to a variation of the horizontal transport.

In the case of Earth, our EBM calculations predict a rise of the mean temperature with increasing D. This is due to the fact that the increased diffusion from the equator to the poles tends to reduce the polar ice covers and, as a consequence, to reduce the albedo and raise the temperature. However, our calculations predict the existence of particular types of climates in which a higher diffusion yields a lower mean temperature. An example is shown in the left panel of Figure 5, where we show the mean annual temperature–latitude profile for a planet at a = 1.0 AU with p = 10^−1.2 bar. Most of the planet is frozen, with the exception of an equatorial belt. In these conditions, an increase in the diffusion will decrease the equatorial temperature and extend the ice cover toward low latitudes; in turn, the increase of the ice cover will cool the planet via albedo feedback. Detailed EBM calculations of this type indicate that the initial ice coverage plays a key role in determining the response of T_m to a variation of the latitudinal transport. When the initial ice cover is modest, an increase of D heats the planet, as in the case of Earth. When the initial ice cover is somewhat above ≈50%, an increase of D might cool the planet, as in the example discussed above; in these cases the cooling may lead to a snowball transition driven by the ice-albedo feedback.

Zoom In Zoom Out Reset image size

The total effect of pressure variations on temperature will depend on the relative strengths of the greenhouse and latitudinal diffusion effects discussed above. If the greenhouse effect dominates, an increase of p will always raise the temperature at a given a, i.e., at a given insolation q = L_⋆/(4πa²). The results that we find in Figure 3 are in line with this expectation when the pressure is sufficiently high, specifically at p ≳ 10^−0.5 bar. However, when the atmospheric column density of greenhouse gases is sufficiently low, pressure variations will not significantly affect the optical depth and temperature. This probably explains the absence of a temperature rise with increasing p in the regime of very low pressures, visible in Figure 4. In this regime, an increase of pressure and diffusion may cool the planet.

As shown in Figure 4, pressure variations not only influence the mean temperature, but also the excursion between minimum and maximum planet temperature. The fact that this excursion becomes increasingly smaller with rising p is due to the increase of diffusion efficiency at high pressures. The high efficiency of the horizontal transport tends to cancel seasonal and latitudinal variations in the surface temperature of the planet. At the other extreme, planets with low atmospheric pressure do not benefit from this heat distribution and undergo large excursions in surface temperature at different latitudes and orbital phases.

For completeness, we mention that pressure variations also affect the thermal inertia of the atmosphere and, as a result, the timescale of adjustment of the planet temperature to seasonal and latitudinal variations of the radiative forcing. In our model, this effect is incorporated in Equation (A1) discussed in Appendix A.1. This effect becomes important only when the fraction of oceans is modest.

4.2.2. Surface Pressure and Planetary Habitability

In addition to the semi-major axis, the orbital eccentricity and the stellar properties can be measured in the framework of observational studies and are relevant for the climate and habitability of exoplanets. The results shown in Figure 3 have been derived for circular orbits and solar-type stars. Here we discuss the extent to which we can generalize these results for eccentric orbits and non-solar-type stars.

4.3.1. Eccentric Orbits

To investigate the dependence of the results on orbital eccentricity we produced a map of habitability in the plane (a, e) for an Earth-like planet with surface pressure p = 1 bar. For the ocean fraction we adopted f_o = 0.7 constant at all latitudes. All the remaining parameters were kept constant, with the values indicated in Table 3. A total number of 2544 simulations were run to cover the interval of semi-major axis 0.6 ⩽ a ⩽ 1.56 with a step δa = 0.02 and the interval of eccentricities 0 ⩽ e ⩽ 0.95 with a step δe = 0.02.

In Figure 6, we show the resulting map of habitability and mean planet temperature in the plane (a, e). A characteristic feature of this map is that the distance of the HZ increases with increasing eccentricity. To interpret this effect, we recall that the mean annual flux received by a planet in an elliptical orbit with semi-major axis a and eccentricity e varies according to the law (see, e.g., Williams & Pollard 2002)

$$\begin{equation} \langle q \rangle = { L_\star \over 4 \pi a^2 (1-e^2)^{1/2}}. \end{equation} \tag{ 10 }$$

Therefore, compared to a circular orbit of radius a and constant insolation q₀ = L_⋆/4πa², the mean annual flux in an eccentric orbit increases with e according to the relation 〈q〉 = q₀(1 − e²)^−1/2. In turn, this increase of mean flux is expected to raise the mean temperature T_m. This is indeed what we find in the results of the simulations. The rise of T_m at constant a can be appreciated in the figure, where the symbols are color-coded according to T_m. To test this effect in a quantitative way, we superimpose on the figure the curve of constant mean flux 〈q〉 = q₀, calculated from Equation (10) for L_⋆ = L_☉. One can see that the HZ follows the same type of functional dependence, e∝(1 − a⁻⁴)^1/2, of the curve calculated at constant flux. This result confirms that the increase of mean annual flux is the main effect that governs the shift of the HZ to larger distances from the star as the eccentricity increases.

Zoom In Zoom Out Reset image size

A second characteristic feature of Figure 6 is that the habitability tends to decrease with increasing eccentricity. This effect is more evident when we consider the map of continuous habitability, h_c, in the right panel of the figure. The effect is related to the large excursion of the instantaneous stellar flux along orbits that are very elongated. The maximum excursion of the flux grows as [(1 + e)/(1 − e)]², and therefore exceeds one order of magnitude when e > 0.5. As a consequence of this strong flux variation, the fraction of orbital period in which the planet is habitable at a given latitude must become increasingly smaller as the orbit becomes more elongated. This is equivalent to saying that f_time(φ) decreases and therefore also h and h_c decrease (Section 3.2) with increasing e. The effect on h_c must be stronger because this quantity depends on f_time(φ) via Equation (9). The comparison between the left and right panels of Figure 6 indicates the existence of an area of the plane (a, e), at high values of a and e, populated by planets that are habitable in small fractions of their orbit. This is demonstrated by the fact that such a population disappears when we consider the continuous habitability h_c. Apart from the existence of this area, it is clear from these figures that the radial extent of the HZ tends to decrease with increasing eccentricity.

On the basis of the above results, it is clear that if the eccentricity is relatively small, the radial extent of the HZ can be scaled according to the mean annual flux. To this end, we introduce the effective circular semi-major axis

$$\begin{equation} a_{{\rm eff}} = a (1-e^2)^{1/4}. \end{equation} \tag{ 11 }$$

Planets in eccentric orbits have a mean annual flux $\langle q \rangle = { L_\star / 4 \pi a_{{\rm eff}}^2}$. In Figure 7, we plot several curves of habitability calculated at different values of eccentricity using a_eff in the abscissa. One can see that, as long as the eccentricity is small (e ≲ 0.5) the curves of habitability versus a_eff calculated at different eccentricities show a good overlap. At high eccentricity the curves of habitability do not overlap. The strong effect of eccentricity on the continuous habitability h_c is evident in the right panel of Figure 7. However, even using the index h_c, the curves of habitability versus a_eff are almost independent of e when the orbits have low eccentricity.

Zoom In Zoom Out Reset image size

We conclude that, in studies of planetary habitability, it is possible to use a_eff as a proxy of a as long as the eccentricity is sufficiently small. In this way, we can explore the effects on habitability of physical quantities other than a and e while varying a single parameter that conveys the information of the orbital characteristics more relevant to the habitability. However, for highly elongated orbits with e > 0.5, one should perform a specific calculation taking into account both a and e.

The influence of orbital eccentricity on planet temperature and habitability has been previously investigated by Williams & Pollard (2002) and Dressing et al. (2010). For comparison with these studies, we calculated the mean global planet temperature as a function of orbital longitude for orbits with e = 0.1, 0.3, and 0.4. As a result, we find seasonal trends similar to those presented in Figure 2 of Williams & Pollard (2002) and Figure 3 of Dressing et al. (2010); temperature differences are present, but generally lower than 5 K.

For comparison with the habitability study of Dressing et al. (2010), we calculated the temporal habitability, f_time, versus latitude in orbits of various eccentricities and increasing semi-major axis. In this way, we used our model to reproduce the left panel of Figure 7 of Dressing et al. (2010). In common with these authors we find that (1) the range of habitable latitudes shrinks as a increases and (2) the temporal habitability undergoes a transition from high to low values above a critical value of a (the "outer edge"). At variance with Dressing et al., we find (1) a slightly larger value of the outer edge (e.g., 1.17 AU instead of 1.12 AU for the case e = 0.6) and (2) a tail of low temporal habitability (f_time ≲ 0.2) outside the "edge."

4.3.2. Non-solar-type Stars

At variance with the eccentricity and the stellar luminosity that can be derived from observational methods, many of the planetary parameters relevant for the climate and habitability are unconstrained by the observations of extrasolar planets. This is true, for instance, for rotation period, axis obliquity, planet geography, and surface albedo. Here we discuss how variations of such unconstrained quantities, in combination with variations of surface pressure, may affect the habitability of exoplanets.

4.4.1. Planet Rotation Period

In Figure 8, we plot the curves of habitability h versus a obtained for three rotation periods: P_rot = 1/3, 1, and 3 days (dashed, solid, and dotted curves, respectively). The four panels of the figure correspond to p = 0.1, 0.3, 1.0, and 3.0 bar, as indicated in the labels.

Zoom In Zoom Out Reset image size

Planet rotation affects the habitability curves. The area subtended by the curves tends to increase and the shape tends to become top-flatted as the rotation period increases. The interpretation of this effect is as follows. An increase of P_rot yields a quadratic increase of the diffusion via Equation (A5) because $D \propto \Omega ^{-2} \propto P_{\rm rot}^{2}$. This effect tends to homogenize the surface temperature, particularly in the high-pressure regime (bottom panels), since the diffusion coefficient increases with p as well. An homogeneous temperature will yield abrupt transitions between a fully habitable and a fully non-habitable situation, the extreme case considered here being the top-flatted, box-shaped habitability curve at P_rot = 3 days and p = 3.0 bar in the bottom right panel of the figure. At the other extreme of low rotation periods and low p, the habitability curves show instead a pronounced peak. Planets in these conditions experience significant variations in their habitability even with a modest change of a (i.e., of insolation). When the pressure is low (top-left panel) the habitability is relatively low. An increase of the rotation period helps to transport the horizontal heat and tends to raise the peak and increase the area subtended by the curve.

The comparison between the results shown in the four panels of Figure 8 indicates that the width and centroids of the curves of habitability tend to change together as pressure varies. All together, variations of the rotation period do not dramatically affect the general features of the pressure-dependent HZ.

Previous studies on the influence of the rotation period on planetary climate and habitability have been presented by SMS08. In common with these studies, we find a decrease in the planet temperature and habitability with increasing rotational velocity. However, our calculations do not yield the runaway transition to a snowball climate that was found by SMS08 by adopting P_rot = 1/3 day for the Earth. The different choice of the albedo prescriptions is a plausible reason for this different result. As one can see in the right panel of Figure 1, the albedo–latitude profile adopted by SMS08 is characterized by a sharp transition at high latitudes. Our model profile, which is in better agreement with the experimental data, does not show such a feature. The sharp transition of the albedo adopted by SMS08 is probably responsible for the very strong ice-albedo feedback found by these authors.

4.4.2. Axis Obliquity

In Figure 9, we plot the curves of habitability h versus a calculated at four values of obliquity: = 0°, 30°, 60°, and 90° (dashed, solid, dash-dotted, and dotted curves, respectively). As in the previous figure, the four panels correspond to p = 0.1, 0.3, 1.0, and 3.0 bar.

Zoom In Zoom Out Reset image size

The curves of habitability show strong variations with obliquity, in particular when the pressure is low. The general trend, at all pressures, is that the habitability increases with increasing when the obliquity is low ( ≲ 60°) and decreases with increasing when the obliquity is high ( ≳ 90°). This complex behavior can be explained in the following way. The configuration at = 0° favors the formation of permanent ice caps in the polar regions, where the star is always at large zenith distance. As the obliquity starts to increase from zero, a larger fraction of polar regions undergo a period of stellar irradiation at low zenith distance in some phase of the orbit. This tends to reduce the ice caps and therefore to increase the habitability. However, when the obliquity becomes quite large, a permanent ice belt starts to build up in the equatorial zones, leading to a decrease in the habitability. To understand why the equator is colder than the poles we consider the extreme case = 90°. In this case, the maximum insolation of a polar region occurs when the pole faces the star; the instantaneous insolation at the pole does not change during the planet rotation, so that the mean diurnal insolation is S = q. The maximum insolation of an equatorial region instead occurs when the rotation axis is perpendicular to the star–planet direction. In this case, the instantaneous insolation at the equator undergoes the night–day cycle and the mean diurnal flux is S = q/π (see Equation (A22) for this specific configuration, in which δ = 0, φ = 0, and H = π/2). Our calculations indicate that the equatorial zones start to build up permanent ice when the obliquity increases from = 60° to = 90°. The exact seasonal evolution is strongly dependent on the thermal inertia of the climate components.

The formation of an equatorial ice belt was first predicted by WK97 for the case = 90° and p = 1 bar. Our calculations indicate that the effect is stronger at lower pressures. The obliquity effects tend to disappear as the pressure increases. The bottom panels of Figure 9 show that at p = 1 bar, the variations are moderate and at p = 3 bar, the influence of obliquity becomes modest. This is due to the high efficiency of the horizontal transport at high p, which tends to cancel temperature gradients on the planet surface, preventing the formation of polar ice caps or an equatorial ice belt at extreme values of obliquity.

The width and centroids of the curves of habitability calculated at different obliquities tend to change together as pressure varies in the four panels of Figure 9. In this respect, we can conclude, as in the case of the rotation period, that variations of the obliquity do not affect the overall characteristics of the pressure-dependent HZ.

Previous studies on the effects of obliquity have been performed by WK97 and Spiegel et al. (2009, hereafter SMS09). For comparison with WK97, we calculated the zonal surface temperature, T, as a function of orbital longitude, L_S, for an Earth with obliquity = 90°. By considering the temperatures in five latitude zones centered on latitudes −85°, −45°, +5°, +45°, and +85°, we obtain trends of zonal T versus L_S very similar to those shown in Figure 1(B) of WK97, with temperature differences generally below 10 K.

For comparison with SMS09, we calculated the temperature–latitude profile of a planet with the same north polar continent considered by these authors, for three different obliquities ( = 235, 60°, and 90°). Our model predicts temperature–latitude profiles similar to those visible in Figure 8 of SMS09, with an important difference: in the case = 90° the temporal evolution that we find is similar to that found by SMS09 only in the first ≃ 10 orbits of the simulation. Afterward, a snowball transition starts to occur in the simulation of SMS09 (bottom-right panel of their Figure 8), but not in our simulation. This result confirms that our climate model is relatively stable against snowball transitions, as we have discussed at the end of Section 4.4.1. Apart from this fact, the influence of obliquity on climate predicted by our EBM is similar to that predicted by previous work.

4.4.3. Planet Geography

Climate EBMs incorporate planet geography in a schematic way, the main parameter being the zonal coverage of oceans, f_o, which at the same time determines the coverage of lands, f_l = 1 − f_o. To explore the effects of geography on habitability we performed two types of tests. First we compared planets with different global coverage of oceans keeping f_o constant in all latitude zones. Then we compared planets with a different location of continents (i.e., polar versus equatorial) keeping constant the global fraction of oceans.

The results of the first test are shown in Figure 10, where we plot the curves of habitability h versus a calculated at three values of ocean fraction: f_o = 0.25, 0.50, and 0.75 (dashed, dotted, and solid curves, respectively). Each panel shows the results obtained at a constant pressure p = 0.1, 0.3, 1.0, and 3.0 bar. Two effects are visible as the coverage of oceans increases: (1) the HZ tends to shift to larger distances from the star; (2) at low pressures (p ≲ 1 bar) the habitability tends to increase. Both effects are relatively small, at least for the combination of parameters considered. The increase in the thermal inertia of the climate system with increasing global ocean fraction is the key to interpreting these results. At low pressure the high thermal inertia of the oceans tends to compensate the inefficient surface distribution of the heat typical of low-pressure atmospheres. At high pressure, the combination of a high thermal inertia with an efficient atmospheric diffusion tends to give the same temperature at a somewhat smaller level of insolation.

Zoom In Zoom Out Reset image size

To investigate the effects of continental/ocean distribution, we considered the three model geographies proposed by WK97: (1) present-day Earth geography, (2) equatorial continent, and (3) polar continent. In practice, each model is specified by a set of ocean fractions, f_o, of each latitude zone (Table III in WK97). The case (2) represents a continent located at latitudes |φ| < 20° covering the full planet. The case (3) represents a polar continent at φ ≲ −30°. The global ocean coverage is approximately the same (〈f_o〉 ≃ 0.7) in the three cases. As a result, we find that these different types of model geographies introduce modest effects on the habitability curves. The habitability of present-day and equatorial continent geography are essentially identical at all pressures. The polar continent geography is slightly less habitable. This is probably due to the combination of two factors that tend to form a larger ice cap in the presence of a polar continent: (1) ice on land has a higher albedo than ice on water and (2) the thermal capacity of continents is lower than that of oceans. In any case, the differences in habitability are small and tend to disappear at p ≳ 3 bar due to the fast rise of the horizontal heat transport.

4.4.4. Albedo of the Continents

4.4.5. Surface Gradient of Latitudinal Heat Transport

4.4.6. Summary

4.5.1. The Edges of the Habitable Zone

4.5.2. Which Type of Life?

The term C represents the thermal inertia of the atmospheric and surface layers that contribute to the surface energy budget. For a layer with density ρ, specific heat capacity c_p, and depth Δℓ, the effective thermal capacity is C = ρ c_p Δℓ. In a planet with oceans, the surface oceanic layers mixed by the winds provide the main contribution to the thermal inertia. A wind-mixed ocean layer with Δℓ = 50 m has C_ml50 = 210 × 10⁶ J m⁻² K⁻¹ (WK97). Also the contribution of the atmosphere can be expressed in terms of an equivalent mixed ocean layer. For the Earth atmosphere, the equivalent ocean depth is 2.4 m (Pierrehumbert 2010), corresponding to C_{atm, ○} = 10.1 × 10⁶ J m⁻² K⁻¹. For other planetary atmospheres the thermal inertia will scale as

$$\begin{equation} \left({ C_{\rm atm} \over C_{{\rm atm},\circ } } \right) = \left({ c_{\rm p} \over c_{{\rm p},\circ } } \right) \, \left({ p \over p_\circ } \right), \end{equation} \tag{ A1 }$$

where c_p and p are the specific heat capacity and total pressure of the atmosphere, respectively (Pierrehumbert 2010). The effective heat capacity of the solid surface is generally negligible. As a representative value, a layer of rock or ice with Δℓ = 0.5 m yields a contribution C_solid ≃ 1 × 10⁶ J m⁻² K⁻¹. This term, even if small, is not negligible in planets without oceans and with very thin atmospheres. On the basis of the above considerations, we adopt

$$\begin{eqnarray} &&C_o= C_{\rm ml50} + C_{\rm atm} \nonumber \\ &&C_l= C_{\rm solid} + C_{\rm atm} \end{eqnarray} \tag{ A2 }$$

for the thermal inertia of the atmosphere over oceans and lands, respectively. For ices not undergoing liquid–solid transition, we adopt C_i = C_l. In the temperature range 263 < T < 273, we add an extra contribution to the thermal inertia of ice to take into account the latent heat of phase transition. In practice, based on similar recipes adopted by North et al. (1983) and WK97, we adopt C_i = C_l + C_ml50/5.

The mean effective thermal capacity of each latitude zone is calculated using as weighting factors the coverage of oceans, f_o, and continents, f_l = 1 − f_o. The expression that we use

$$\begin{equation} C = f_l [(1-f_{i}) \, C_l+f_{i} \, C_{i}] + f_o [ (1-f_{i}) C_o + f_{i} \, C_{i}] \end{equation} \tag{ A3 }$$

takes into account the ice coverage, f_i, estimated as explained in Appendix A.6.

All together, our values of thermal capacity are very similar to those adopted by WK97 and SMS08, with the exception of C_{atm, ○}, which is about twice in our case. Attempts to use higher values of the ocean capacity previously proposed by North et al. (1983) and representative of a 75 m deep wind-mixed ocean layer, yield a worse match between the model predictions and the Earth experimental data.

The horizontal transport of heat on the planet surface is mainly due to the general circulations and the related instabilities of the ocean and the atmosphere, which are governed by many factors. The most relevant are the physical and chemical properties of the atmosphere and oceans, the presence of Coriolis force induced by the planet rotation, and the topography of the planet surface. Since EBMs are zonally averaged, longitudinal variations of the heat transport are not considered. Even so, it is extremely difficult to model all the factors that govern the latitudinal transport. To keep low the computational cost, EBMs incorporate the efficiency of the transport into a single quantity, namely the diffusion coefficient D, even though representing latitudinal heat transport by a diffusive term is, in itself, a gross simplification. In previous EBM work, this term has been parameterized using a scaling law that takes into account the dependence of the diffusion coefficient on the main physical quantities involved in the latitudinal transport. The relation proposed by WK97,

$$\begin{equation} \left({D \over D_\circ } \right) = \left({ p \over p_\circ } \right) \, \left({ c_\mathrm{p} \over c_{{\rm p},\circ }} \right) \, \left({ m \over m_\circ } \right)^{\!\! -2} \, \left({ \Omega \over \Omega _\circ } \right)^{\!\! -2} , \end{equation} \tag{ A4 }$$

gives D as a function of atmospheric pressure, p, specific heat capacity, c_p, mean molecular weight, m, and angular velocity of planet rotation, Ω. All these quantities are scaled with respect to the reference Earth values. A discussion of the physics underlying this relation can be found in WK97. A number of recent applications of EBM to exoplanet climates have adopted the same relation (SMS08; Dressing et al. 2010). We refer to these papers for cautionary remarks regarding the limits of validity of the scaling relation D∝Ω⁻².

In our work, we address an intrinsic limitation of the above formulation, namely the fact that, once the planet parameters are specified, D is constant in latitude and time. This approximation is too crude because Earth's latitudinal transport is driven by atmospheric convective cells characterized by a latitudinal pattern with seasonal variations. As a consequence, in real planets the efficiency of the diffusion will vary in latitude and time. Here we propose a formalism that introduces latitudinal and seasonal variations of D related to the stellar zenith distance. This formalism can be applied to a generic planet.

Most of the latitudinal heat transport is carried out by the circulation that takes place in the atmospheric convective cells. On Earth, the most notable of these features are the two Hadley cells that depart from the intertropical convergence zone (ITCZ) and produce the net effect of a polarward transport. The position of the Hadley cells shifts during the year, influenced by the solar zenith distance. This is demonstrated by the seasonal shift of the ITCZ, which moves to higher latitudes in the summer hemisphere. Since the Hadley cells yield the largest contribution to the latitudinal heat transport, we conclude that the maximum efficiency of this type of transport is larger when the zenith distance Z_⋆ is smaller. Inspired by this behavior, and without making specific assumptions on the systems of atmospheric cells that may be present in a generic planet, we have introduced a dependence of the diffusion coefficient on μ = cos Z_⋆. More specifically, we assume that D scales with the mean diurnal value $\overline{\mu }$, derived from Equation (A23). From this equation and Equation (A25) one can see that $\overline{\mu }$ depends on the latitude, φ, and the obliquity of the axis of rotation, ε; it also depends on the time of the year, t, which determines the seasonal value of δ_⋆ and λ_⋆. Therefore, by assuming a dependence $D=D(\overline{\mu })$ we introduce, in practice, a dependence on φ, t and ε. To incorporate this effect in the EBM prescriptions, we adopt the following expression:

$$\begin{equation} \left({D \, \over D_\circ } \right) = \zeta _\varepsilon (\varphi ,t) \, \left({ p \, \over p_\circ } \right) \, \left({ c_\mathrm{p} \over c_{{\rm p},\circ }} \right) \, \left({ m \, \over m_\circ } \right)^{\!\! -2} \, \left({ \Omega \, \over \Omega _\circ } \right)^{\!\! -2} , \end{equation} \tag{ A5 }$$

where ζ_ε(φ, t) is a modulating factor that scales linearly with $\overline{\mu }$. This factor is normalized in such a way that the mean global annual value of D equals that given by the scaling law (A4). Thanks to this fact, we can make use of the previous studies of the diffusion coefficient, as far as the dependence on p, c_p, m, Ω is concerned, and of our own formalism, as far as the dependence on φ, t, and ε is concerned. We now show how the modulating factor ζ_ε(φ, t) can be expressed as a function of a single parameter, $\mathcal {R}$, that represents the ratio between the maximum and minimum values of the diffusion coefficient in the planet.

A.2.1. The Modulating Factor ζ

To obtain an expression for the modulating factor, ζ, we use the following defining conditions: (1) ζ scales linearly with $\overline{\mu }$; and (2) ζ is normalized in such a way that its mean global annual value is 〈ζ_ε(φ, t)〉 = 1. The first condition translates into

$$\begin{equation} \zeta _\varepsilon (\varphi ,t) = c_0 + c_1 \, \overline{\mu }_\varepsilon (\varphi ,t) , \end{equation} \tag{ A6 }$$

where c₀ and c₁ are two constants. The second condition becomes

$$\begin{equation} \langle \zeta _\varepsilon (\varphi ,t) \rangle = c_0 + c_1 \, \langle \overline{\mu }_\varepsilon (\varphi ,t) \rangle = 1 , \end{equation} \tag{ A7 }$$

where

$$\begin{equation} \langle \overline{\mu }_\varepsilon (\varphi ,t) \rangle = { \int _0^P dt \int _{-{\pi /2}}^{+{\pi /2}} d\varphi \cos \varphi \overline{\mu }_\varepsilon (\varphi ,t) \over \int _0^P dt \int _{-{\pi /2}}^{+{\pi /2}} d\varphi \cos \varphi }. \end{equation} \tag{ A8 }$$

To determine c₀ and c₁, we need an additional relation between these two constants. For this purpose, we follow an approach similar to that of North et al. (1983), who adjusted the diffusion coefficient to be three times as large at the equator as it is at the poles. In the framework of our formulation, this condition would be expressed in the form $\mathcal {R}=3$, where

$$\begin{equation} \mathcal {R}={ [\zeta _\varepsilon (\varphi ,t)]_{\max } \over [\zeta _\varepsilon (\varphi ,t)]_{\min } } \end{equation} \tag{ A9 }$$

is the ratio between the maximum and minimum values of the modulation factor in any latitude zone and at any orbital phase. The above condition is equivalent to

$$\begin{equation} \mathcal {R}={ c_0 + c_1 \, [\overline{\mu }_\varepsilon (\varphi ,t)]_{\max } \over c_0 + c_1 \, [\overline{\mu }_\varepsilon (\varphi ,t)]_{\min } } = { c_0 + c_1 \, [\overline{\mu }_\varepsilon (\varphi ,t)]_{\max } \over c_0 } , \end{equation} \tag{ A10 }$$

where we have used the fact that $[\overline{\mu }_\varepsilon (\varphi ,t)]_{\min }=0$ since $\overline{\mu }=0$ in the zones where the star lies below the horizon during a period of rotation. We treat $\mathcal {R}$ as an input parameter of the model used to estimate c₀ and c₁. By combining Equations (A7) and (A10), it is easy to show that

$$\begin{equation} c_1 = \left({ [\overline{\mu }_\varepsilon (\varphi ,t)]_{\max } \over \mathcal {R}-1 } + \langle \overline{\mu }_\varepsilon (\varphi ,t)\rangle \right)^{-1} \end{equation} \tag{ A11 }$$

and

$$\begin{equation} c_0 = c_1 \times { [\overline{\mu }_\varepsilon (\varphi ,t)]_{\max } \over \mathcal {R}-1 }. \end{equation} \tag{ A12 }$$

The quantities $[\overline{\mu }_\varepsilon (\varphi ,t)]_{\max }$ and $\langle \overline{\mu }_\varepsilon (\varphi ,t)\rangle$ are calculated given the obliquity ε. To compute the double integral in Equation (A8) we use the expression of $\overline{\mu }$ given in Equation (A23), with the condition $\overline{\mu }=0$ when H ⩽ 0, i.e., when the Sun is below the local horizon for a complete rotation of the planet.

A.2.2. Calibration of $\mathcal {R}$

In climate EBMs, the OLR is usually expressed as a function of the temperature. In our problem, we are interested in probing variations of total pressure and therefore we need an expression I = I(T, p). The pressure dependence of the OLR enters through the infrared optical depth of the atmosphere, τ_IR, which governs the intensity of the greenhouse effect. In principle, one may expect a strong dependence of τ_IR on p as a result of two effects. First, in a planetary atmosphere τ_IR = κ_IR p/g, where κ_IR is the total absorption coefficient and g is the gravitational acceleration. Second, in the range of planetary surface pressures, κ_IR is dominated by collisional broadening, which introduces a linear dependence on p of the widths of the absorption lines (Pierrehumbert 2010; Salby 2012). In the absence of line saturation, κ_IR should therefore increase linearly with p. All together, one might expect τ_IR = κ_IR p/g∝p². In practice, however, the rise of the optical depth is milder because (1) the absorptions of different species may overlap in wavelength, (2) lines are partly saturated, and (3) the amount of water vapor, which is an important contributor to κ_IR, is independent of p. Given the complex interplay of these factors, it is not possible to derive a simple analytical function τ_IR = τ_IR(T, p). Therefore, in order to obtain I = I(T, p), we performed a series of radiative calculations and tabulated the results as a function of T and p.

We used standard radiation models developed at the National Center for Atmospheric Research (NCAR), as part of the Community Climate Model (CCM) project. The calculations were performed in a column radiation model scheme for a cloud-free atmosphere (Kiehl & Briegleb 1992). We used the set of routines CliMT (Pierrehumbert 2010; Caballero 2012), adopting Earth's value of surface gravitational acceleration. We varied p while keeping the mixing ratios of non-condensable greenhouses gases (CO₂ and CH₄) equal to Earth's values. The reference values of total pressure and partial pressure of greenhouse gases are shown in Tables 1 and 2. The contribution of water vapor was parameterized through the relative humidity, rh, as we explain below. In the regime of high temperatures, we increased the resolution of the vertical strata (up to 10,000), in order to better track the water vapor in the higher atmospheric levels. To test the radiative calculations in the regime of very low pressure (i.e., negligible atmospheric greenhouse), we compared the predicted OLR with the blackbody radiation calculated at the planet surface. As expected, the results are identical as long as water vapor is negligible.

To calibrate the OLR for the case of Earth, we used as a reference the mean annual Earth OLR measured with the Earth Radiation Budget Satellite (ERBE) satellite (online material published by Pierrehumbert 2010). In practice, we combined two ERBE data sets—i.e., the OLR versus latitude and the surface temperature versus latitude—to build a set of OLR data versus temperature. These experimental data are shown with crosses in Figure 11. To test different types of models, we then used the CCM to calculate the clear-sky OLR as a function of temperature for different values of relative humidity (rh = 0.1, 0.5, and 1.0). We then subtracted 28 W m⁻² to the clear-sky results to take into account the mean global long-wavelength forcing of the clouds on Earth (Pierrehumbert 2010). The three CCM curves resulting for rh = 0.1, 0.5, and 1.0 are shown in Figure 11. The curve obtained for rh = 0.5 (solid line) gives the best match to the experimental data and is very similar to the OLR adopted by SMS08 and SMS09 (dashed line). Based on these results, we adopted rh = 0.5 in our model. For applications to planets with cloud coverage 〈f_c〉 different from that of Earth, 〈f_{c, ○}〉, we subtract 28 (〈f_c〉/〈f_{c, ○}〉) W m⁻² to the clear-sky OLR obtained from the radiative calculations.

Zoom In Zoom Out Reset image size

In the formulation of SMS08, the albedo is an analytical function of the temperature, which, in practice, considers two types of surfaces: with ice and without ice. In our work, we employ a formalism that can be applied to planets with any type of surface characteristics. For each latitude zone we calculate a mean albedo by weighting the contributions of continents, a_l, oceans, a_o, clouds, a_c, ice on continents, a_il, and ice on oceans, a_io. The weighting factors are the zonal coverage of oceans, f_o, and continents, f_l; within oceans and continents we separate the contribution of ice, f_i; the cloud coverage on water, land, and ice is specified by the parameters f_cw, f_cl, and f_ci, respectively. In this way, we obtain

$$\begin{eqnarray} a_s &=& f_o \lbrace (1-f_{i}) [ a_{o} (1-f_{{\rm cw}}) + a_{c}f_{{\rm cw}} ] \nonumber \\ &&+ f_{i} [ a_{{\rm io}} (1-f_{{\rm ci}}) + a_{c} f_{{\rm ci}} ] \rbrace \nonumber \\ &&+ f_l \lbrace (1-f_{i}) [ a_{l} (1-f_{{\rm cl}}) + a_{c} f_{{\rm cl}} ] \nonumber \\ &&+ f_{i} [ a_{{\rm il}} (1-f_{{\rm ci}}) + a_{c} f_{{\rm ci}}] \rbrace . \end{eqnarray} \tag{ A13 }$$

In this equation, we have omitted for simplicity the latitude dependence of the zonal coverage of oceans and lands. This formulation of the surface albedo is similar to that adopted by WK97. At variance with that work, we consider the ice on continents (not only on oceans) and variable cloud cover. Only a few parameters in the above expression are free. Most of the albedo parameters are estimated with specific prescriptions that we describe below. The zonal coverages of oceans, lands, ice, and clouds are described in Appendix A.6. In practice, only the zonal fraction of oceans, f_o, and the continental albedo, a_l, are free parameters.

A.4.1. Albedo Parameters

For the albedo of the continents, we adopt a_l = 0.2, the same value used by WK97, as a representative value of the Earth continents. For other planets we treat a_l as a free parameter that can typically vary between 0.1 and 0.35 (see Section  4.4.4).

The albedo of the oceans plays a crucial role in the estimate of a_s. The ocean reflection depends on the zenith distance of the star, Z_⋆, that the EBM provides self-consistently as a function of time and latitude. WK97 adopted the Fresnel formula (Kondratyev 1969) to model the ocean reflectivity. However, the Fresnel formula is an ideal approximation valid for smooth surfaces while the actual sea level is never ideally smooth. For this reason we opted for a function calibrated with experimental data (Briegleb et al. 1986; Enomoto 2007):

$$\begin{eqnarray} a_o &=& { 0.026 \over (1.1 \, \mu ^{1.7} + 0.065)} \nonumber \\ && + 0.15 (\mu -0.1) \, (\mu -0.5) \, (\mu -1.0) , \end{eqnarray} \tag{ A14 }$$

where μ = cos Z_⋆ is calculated with Equation (A20). Compared to the Fresnel formula, this expression gives very similar results at Z_⋆ ≲ 40°, but a less steep, more realistic, rise at high values of zenith distance (Enomoto 2007).

For the albedo of ice over lands and ocean we adopted a_il = 0.85 and a_io = 0.62, respectively. The existence of a difference between these two types of albedo has been documented in previous work (Kondratyev 1969; Pierrehumbert 2010). By adopting a high albedo for ice over continents, we improve the modelization of the Earth climate on Antarctica, which is a critical point of EBMs (see Figure 1 and Section 2.3).

The cloud albedo, a_c, depends on the microphysical properties and geometrical dimensions of the clouds, which determine their optical properties. For a cloud of given optical thickness, the albedo increases with increasing stellar zenith angle, Z_⋆, which elongates the slant optical depth (Salby 2012). A linear dependence of a_c on Z_⋆ was reported by Cess (1976) for a set of data representative of the Earth cloud albedos. On the basis of that result, WK97 adopted a_c = α + βZ_⋆ and tuned the parameters α and β in such a way so as to match the total Earth albedo with their model. Since this tuning yields α < 0, a problem with this formalism is that the cloud albedo vanishes at low zenith distances and becomes negative at Z_⋆ = 0. We expect, instead, a minimum value of cloud albedo to exist at Z_⋆ = 0, corresponding to the minimum vertical thickness of the cloud. For this reason, in our work we adopt the same parameterization, but with a minimum value of cloud albedo at low zenith distances, a_c0. In practice, we use the expression

$$\begin{equation} a_c = \max \lbrace a_{c0}, [ \alpha + \beta Z_\star ] \rbrace . \end{equation} \tag{ A15 }$$

We estimate α and β by requiring the relation α + βZ_⋆ to yield a good fit to Cess data. We then tune a_c0 in such a way so as to match the total Earth albedo with our model. In this way we obtained α = −0.07, β = 8 × 10⁻³(°)⁻¹, and a_c0 = 0.19.

The incoming stellar radiation is calculated as the diurnal average

$$\begin{equation} S= { \int _{0}^{2\pi } s(t_\star) \, dt_\star \over \int _{0}^{2\pi } dt_\star } \end{equation} \tag{ A16 }$$

of the stellar flux s(t_⋆) incident on the planet at latitude φ, where t_⋆ is the instantaneous hour angle of the star measured in angle units from the local meridian. To estimate s(t_⋆) we proceed as follows. At the time t, the planet is located at a distance r = r(t) from its star. The stellar flux at such distance is

$$\begin{equation} q(r) = { L_{\star } \over 4 \pi r^2 } , \end{equation} \tag{ A17 }$$

where L_⋆ is the bolometric luminosity of the star. If we call Z_⋆ the stellar zenith angle, the instantaneous flux at the planet surface is

$$\begin{equation} s(t_\star) = \left\lbrace \begin{array}{@{} ll}\tau _a \, q(r) \, \cos Z_\star & {\rm if} {|t_\star | < H} \\ 0 & {\rm if} {|t_\star | \ge H} \end{array}\right. \end{equation} \tag{ A18 }$$

where τ_a is the short-wavelength transmissivity of the atmosphere and H is the half-day length. By definition, −H ⩽ t_⋆ ⩽ +H represents the portion of rotation period during which the star stays above the horizon. The half-day length in radians can be estimated from the expression (WK97)

$$\begin{equation} \cos H = {-}\tan \varphi \tan \delta (0 < H < \pi). \end{equation} \tag{ A19 }$$

The flux s(t_⋆) is null when |t_⋆| ⩾ H because in this case the star is below the horizon (Z_⋆ ⩾ π/2). The zenith distance, Z_⋆, is related to the latitude, φ, the stellar declination, δ_⋆, and the hour angle, t_⋆, by means of the equation

$$\begin{equation} \cos Z_\star = \sin \varphi \, \sin \delta _\star + \cos \varphi \, \cos \delta _\star \, \cos t_\star. \end{equation} \tag{ A20 }$$

With the above relations, we can calculate the average (A16) along circles of constant latitude. In these circles, the terms sin φ and cos φ are constant in the integration. Also the declination δ_⋆ = δ_⋆(t) and the flux q(r) = q(r[t]) can be treated as constants if the rotation period is much smaller than the orbital period, i.e., if

$$\begin{equation} P_\mathrm{rot} \ll P_\mathrm{orb}. \end{equation} \tag{ A21 }$$

With these assumptions, if the atmosphere is transparent in the short wavelength range (τ_a = 1), it is easy to show that

$$\begin{equation} S= { q(r) \over \pi } \left(H\sin \varphi \, \sin \delta _\star + \cos \varphi \, \cos \delta _\star \, \sin H \right). \end{equation} \tag{ A22 }$$

In a similar fashion, by defining μ = cos Z_⋆, it is possible to calculate the mean diurnal value of μ when the star is above the horizon, $\overline{\mu }=\int _{-H}^{H} \mu \, dt_\star /\int _{-H}^{H} dt_\star$, and obtain the relation

$$\begin{equation} \overline{\mu } = \sin \varphi \, \sin \delta _\star + \cos \varphi \, \cos \delta _\star \, { \sin H \over H } \end{equation} \tag{ A23 }$$

that we use in our formulation of the diffusion coefficient.

Once we have obtained the diurnal average S, we calculate its temporal variation in the course of the orbital period, S = S(t). The time t specifies the planet position along the orbit and S = S(t) represents the seasonal evolution of the flux at a given latitude. The only quantities that depend on t in Equation (A22) are q(r) and δ_⋆. To calculate the seasonal evolution of q(r) and δ_⋆ we proceed as follows.

We specify the position of the planet on its orbit using a system of polar coordinates centered on the star and origin of the angles in the direction to the pericenter, as shown in Figure 12. The position of the planet at time t, P_t, is specified by r = r(t) and the true anomaly, ν = ν(t). The instant stellar flux is given by q(r) = (q₀/(r/a)²) where q₀ = L_⋆/(4πa²). The instantaneous value of r/a can be calculated by introducing an additional angular variable E (the eccentric anomaly) such that

$$\begin{equation} r/a = 1-e \cos E , \end{equation} \tag{ A24 }$$

where e is the eccentricity and E is interpreted geometrically in the left panel of Figure 12. We show below how to compute E.

Zoom In Zoom Out Reset image size

To calculate the stellar declination δ_⋆ we consider the planetocentric orbital longitude of the star, λ_⋆, shown in the right panel of Figure 12. The origin of this angle is the line of nodes, defined by the intersection of the orbital plane with the equatorial plane. Without loss of generality, we set t = 0 the instant in which the planet crosses the ascending node.⁵ With this choice, the instant value of the stellar declination is (Allison & McEwen 2000)

$$\begin{equation} \delta _\star = \arcsin \left(\sin \varepsilon \, \sin \lambda _\star \right) , \end{equation} \tag{ A25 }$$

where ε is the obliquity, i.e., the inclination of the planet's orbit to its equator. To calculate λ_⋆(t) we consider its relation with the true anomaly, ν. From the top panel of Figure 12, it is easy to see that

$$\begin{equation} \lambda _\star (t) = \nu (t)+\lambda _\mathrm{P} , \end{equation} \tag{ A26 }$$

where λ_P is the planetocentric longitude of the star at the moment in which the planet is at the pericenter (P₁ in the figure); in terms of orbital parameters, λ_P = ω, where ω is the argument of the pericenter. To solve Equation (A26), we use the expression (Bertotti & Farinella 1990, Equation (10.30), p. 211)

$$\begin{equation} \tan (\nu /2) = \left({ 1+e \over 1-e} \right)^{1/2} \tan (E/2) \end{equation} \tag{ A27 }$$

that can be inferred from the top panel of Figure 12.

At this point we are left with the calculation of the eccentric anomaly E. To this end we use Kepler's equation

$$\begin{equation} E-e \sin E = M , \end{equation} \tag{ A28 }$$

where M is the mean anomaly, defined to be a linear function of time which increases by 2π per revolution according to the expression

$$\begin{equation} M= \overline{n} \, t + M_\circ , \end{equation} \tag{ A29 }$$

where $\overline{n}=2\pi /P$. The constant M₀ is determined by the choice of the initial conditions. To solve the transcendent Kepler's equation (A28), we use Newton's iteration method.⁶

The zonal ocean fraction, f_o, is a free parameter that also determines the fraction of continents f_l = 1 − f_o. By assigning f_o and f_l to each latitude zone of the planet, the EBM takes into account the planet geography, although in a schematic way.

The zonal coverage of ices, f_i, is calculated as a function of the mean diurnal zonal temperature, T. This approach, adopted by WK97 and followed by SM08, allows to incorporate the climate feedback between temperature and ice albedo into the EBM. In their original formulation, WK97 used data from Thompson & Barron (1981) to calculate the ocean fraction covered by sea ice as

$$\begin{equation} f_i (T)= \max \lbrace 0, [ 1 - e^{ (T-273\,\mathrm{K})/ 10\,\mathrm{K} }]\rbrace . \end{equation} \tag{ A30 }$$

This expression provides a "freezing curve" that ideally describes the formation of ices as the temperature decreases. However, in this formulation the ice melts completely and instantaneously as soon as T > 273 K. This approximation is too crude for zones that are frozen for most of the time: in such regions we expect a buildup of permanent ices if the timescales for ice formation and ice melting are comparable. To avoid this problem, we compare the time intervals in which T < 273 K and T ⩾ 273 K in each zone. If T < 273 K during more than 50% of the orbital period, we adopt a constant ice coverage for the full orbit, $f_i=f_i (\overline{T})$, where $\overline{T}$ is the mean annual zonal temperature. In the other cases, we follow the original formulation (A30). In this way we obtain the formation of permanent ices, with a coverage that increases as the mean zonal annual temperature decreases. Our treatment significantly improves the description of the ice coverage on Earth, avoiding a sudden appearance and disappearance of the polar caps in the course of each orbit. To keep the model simple, we adopt for the continents the same ice coverage of the oceans. A more realistic treatment of the ice coverage would require a formulation of the formation and melting of ices as a function of time. The dynamical treatment of the ice formation and melting is beyond the purpose of this paper and will be the subject of subsequent work. We refer to Spiegel et al. (2010) for a detailed discussion of ice formation and melting using an EBM.

As far as the cloud coverage is concerned, we adopt different values for clouds over oceans, continents, and ices in our formulation of the albedo (Equation (A13)). The observational support for variations of the cloud coverage on different types of underlying surfaces comes from an analysis of Earth data performed by Sanromá & Pallé (2012). These authors used data collected by the International Satellite Cloud Climatology Project and presented the cloud coverage as a function of latitude for clouds on water, ice, desert, and vegetation. We averaged these latitude profiles to obtain a global representative cloud coverage for each type of surface. The value for the water was slightly adjusted to improve the agreement of our model with the experimental albedo–latitude profiles of Earth. The cloud coverage over continents was calculated by weighting deserts and vegetation with factors 2/3 and 1/3, respectively.⁷ As a result, we adopt f_cw = 0.67, f_cl = 0.50, and f_ci = 0.50 for clouds over water, continents, and ices, respectively. With this choice of parameters, we obtain a global cloud coverage 〈f_c〉 = 0.612 in our best model for Earth. This value is in agreement with the experimental value 〈f_c〉 = 0.603 obtained from the ERA Interim reanalysis for the years 1979–2010 (Dee et al. 2011). For comparison, WK97 adopted 〈f_c〉 = 0.5 in their formulation of the surface albedo of Earth. An interesting feature of our formalism is that the cloud coverage is automatically adjusted for planets with different fractions of continents, oceans, and ices.