Limit cycles can reduce the width of the habitable zone
The liquid water habitable zone (HZ) describes the orbital distance at which a terrestrial planet can maintain above-freezing conditions through regulation by the carbonate-silicate cycle. Recent calculations have suggested that planets in the outer regions of the habitable zone cannot maintain stable, warm climates, but rather should oscillate between long, globally glaciated states and shorter periods of climatic warmth. Such conditions, similar to ‘Snowball Earth’ episodes experienced on Earth, would be inimical to the development of complex land life, including intelligent life. Here, we build upon previous studies with an updated an energy balance climate model to calculate this ‘limit cycle’ region of the habitable zone where such cycling would occur. We argue that an abiotic Earth would have a greater CO partial pressure than today because plants and other biota help to enhance the storage of CO in soil. When we tune our abiotic model accordingly, we find that limit cycles can occur but that previous calculations have overestimated their importance. For G stars like the Sun, limit cycles occur only for planets with CO outgassing rates less than that on modern Earth. For K and M star planets, limit cycles should not occur; however, M-star planets may be inhospitable to life for other reasons. Planets orbiting late G-type and early K-type stars retain the greatest potential for maintaining warm, stable conditions. Our results suggest that host star type, planetary volcanic activity, and seafloor weathering are all important factors in determining whether planets will be prone to limit cycling.
Earth’s orbit falls within the boundaries of the habitable zone (HZ) where a rocky planet can maintain liquid water on its surface, given a CO-N-HO atmosphere and some mechanism (e.g., plate tectonics) for recycling these volatiles (Kasting et al., 1993; Abe et al., 2011; Kopparapu et al., 2013, 2014; Leconte et al., 2013; Wolf & Toon, 2014; Yang et al., 2014). The inner HZ edge is bounded by either the runaway greenhouse effect, in which liquid water evaporates entirely, or the moist greenhouse111We use the term ‘moist greenhouse’ throughout our paper, although we acknowledge that this may be imprecise terminology, as Earthâs present atmosphere could be described as a moist greenhouse. Other terms such as ‘moist stratosphere’ or ‘diffuse tropopause’ might be better descriptors of this phenomenon, but they suffer from sounding arcane or obscure. effect, in which liquid water persists on a planet’s surface but the stratosphere becomes wet and water is lost by photodissociation followed by escape of hydrogen to space. The moist greenhouse effect is difficult to simulate with one-dimensional models that assume a saturated troposphere (Kopparapu et al., 2013), but it does appear in some general circulation models (Wolf & Toon, 2014). Planets farther out in the HZ are expected to accumulate dense CO atmospheres because of the negative feedback between silicate weathering (the loss process for atmospheric CO) and surface temperature (Walker et al., 1981). But this feedback loop is ultimately limited by CO condensation and Rayleigh scattering, which combine to create an outer HZ boundary termed the ‘maximum greenhouse limit’ (Kasting et al., 1993). A proposed extension of the outer HZ boundary by formation of CO ice clouds (Forget & Pierrehumbert, 1997; Mischna et al., 2000; Colaprete & Toon, 2003; Forget et al., 2013) appears less likely when the ‘scattering greenhouse effect’ of these clouds is recomputed using more accurate radiative transfer models, (Kitzmann, 2016).
The conventional thinking regarding the HZ outer edge may be too optimistic, however, because it fails to account for mass transfer rates of CO. CO is released from volcanoes and is consumed by silicate weathering followed by deposition of carbonate sediments (Walker et al., 1981; Berner et al., 1983). These processes are in approximate balance on modern Earth, creating a climate that is stable and relatively warm, even if it is sometimes perturbed by glacial-interglacial cycles. Occasional ‘Snowball Earth’ episodes in which the planet is fully glaciated (Hoffman et al., 1998) have been attributed to a variety of complicating factors, including changes in atmospheric O and CH (Pavlov et al., 2003), as well as limit cycles involving atmospheric CO (Tajika, 2007; Mills et al., 2011). Limit cycles—oscillations between ice-free and globally glaciated states—occur in models of the early Earth in which volcanic outgassing rates are too low to sustain a CO-warmed climate (Tajika, 2007).
Several new papers have argued that such limit cycle behavior could be far more prevalent than previously thought (Kadoya & Tajika, 2014, 2015; Menou, 2015). Depending on the volcanic outgassing rate and host stellar type, some Earth-like planets are subject to climatic limit cycles when stellar insolation is low, as it is in the outer parts of a star’s HZ. All of these recent studies employed parameterized versions of energy-balance climate models (EBMs). Such models include radiative balance between incident stellar and outgoing infrared radiation, along with diffusional heat fluxes between different latitude bands. These models are useful tools for studying climate variations that occur on time scales of thousands to tens of millions of years.
Here we discuss the possibility that limit cycles could reduce the outer edge of the HZ. We use our own EBM, which implements an updated parameterization of radiative transfer based on 1-D radiative-convective HZ calculations (Kopparapu et al., 2013, 2014). Our model includes a representation of the carbonate-silicate cycle used by previous studies (Menou, 2015), and also expands on previous work to include the effect of CO condensation and the impact of seafloor weathering. These improvements allow us to determine, for the first time, the limit cycle boundaries relative to the conventional liquid water HZ for different types of stars.
2 Model Description
Latitudinal energy balance models (EBMs) are computationally efficient models that are well-suited to exploring glacial cycling and its influence on climate. Although EBMs can be useful with a highly idealized linearization of infrared absorption (North et al., 1981; Gaidos & Williams, 2004; Haqq-Misra, 2014), investigation of HZ limits requires an EBM with a more sophisticated radiative transfer parameterization (Williams & Kasting, 1997; Fairén et al., 2012; Vladilo et al., 2013; Kadoya & Tajika, 2014; Menou, 2015). However, no previously published EBM has yet been updated to include parameterizations based upon radiative-convective HZ model calculations (Kopparapu et al., 2013, 2014), which use coefficients derived from the HITRAN 2008 (Rothman et al., 2009) and HITEMP 2010 (Rothman et al., 2010) spectroscopic databases to account for additional absorption features of HO and CO compared to earlier models (Kasting et al., 1993). Existing EBMs that consider planetary habitability (Vladilo et al., 2013; Kadoya & Tajika, 2014, 2015; Menou, 2015) use an older radiation scheme (Williams & Kasting, 1997), which is based upon a polynomial parameterization of prior radiative-convective calculations (Kasting, 1991). These EBM investigations continue to be useful for interpreting exoplanet habitability and guiding more complicated modeling studies toward interesting regions of parameter space, but the lack of up-to-date radiative transfer parameterization in such models can limit their utility when comparing with other recent climate models or interpreting exoplanet observations.
We use a one-dimensional EBM that has been developed in previous studies (Williams & Kasting, 1997; Gaidos & Williams, 2004; Fairén et al., 2012; Haqq-Misra, 2014). This EBM calculates meridionally averaged temperature as a function of latitude and time according to the equation
Eq. (1) expresses the change in temperature as the sum of stellar heating, infrared cooling, and meridional diffusion. Diurnally averaged solar flux is the product of a constant solar flux and a function of latitude (which assumes a circular orbit) to yield seasonally varying insolation (Gaidos & Williams, 2004). The effective heat capacity of the surface and atmosphere, , depends on the fraction of ocean and ice coverage at a given latitude (Williams & Kasting, 1997; Fairén et al., 2012). We prescribe Earth-like conditions by setting 70% ocean coverage at all latitudes and allowing fractional ice coverage between 263 K to 273 K. Radiative fluxes are represented by the top-of-atmosphere albedo and the infrared outgoing radiative flux . Our FORTRAN implementation of this EBM is discretized into 18 equally spaced latitudinal zones with a initial uniform temperature profile of and stepped through 1000 (or more) complete orbits by increments of to numerically solve Eq. (1) with a forward finite differencing scheme222We can be confident that our choice of parameters will yield a converged solution by examining the CFL condition for numerical stability, , where for explicit time-marching solvers like ours. If we approximate our length interval as the radius of Earth divided by the total number of latitudinal bands, and if we assume a typical advective wind speed on Earth of m s, then . Only when the number of latitude bands approaches 75 does .. We initialize all simulations with snowball conditions () and CO partial pressure CO.
The diffusive parameter, , describes the efficiency of meridional energy transport and scales with changes in atmospheric pressure, heat capacity, atmospheric mass, and rotation rate (Williams & Kasting, 1997). This parameter accounts for the exchange of sensible and latent energy fluxes between the tropics and midlatitudes by adjusting the energy balance equilibrium temperature at each latitude. Williams & Kasting (1997) assumed that is proportional to the inverse square of rotation rate, so that a more rapidly (slowly) rotating planet will show a decreased (increased) tendency toward meridional energy distribution. This behavior mimics well-known results from general circulation models (Williams & Halloway, 1982; Showman et al., 2013) while still maintaining the computational efficiency of an EBM. We focus our study on a model planet with a rotation rate equal to present-day Earth, but this assumption does not significantly alter our results. A more slowly rotating planet would have increased meridional energy transport and would therefore tend toward uniform conditions at all latitudes; such a limit is analogous to the globally-averaged EBM used by Menou (2015). A more rapidly rotating planet would show more extreme variations between equatorial and polar temperatures, but this contrast does not significantly alter the onset of global glaciation and deglaciation events. Sensitivity studies indicate that either of these cases will still exhibit transitions into a limit cycle state at approximately the same value of solar constant, so our assumption of a fixed present-day rotation rate does not necessarily limit the scope of our calculations.
Surface albedo is calculated as a weighted sum of the albedos of unfrozen land, unfrozen ocean, and ice coverage at each latitude. We assume an albedo of 0.2 for unfrozen land, while the albedo of ocean varies as a function of solar zenith angle (Williams & Kasting, 1997; Fairén et al., 2012). We also implement a band-dependent ice albedo (Pollard & Kasting, 2005) that partitions the broadband ice albedo into a visible component, , and a near-infrared component, . This latter value applies to snow-covered ice, as we discuss further below. According to Shields et al. (2013), the near-IR albedo of bare ice is 0.3. Our model does not include a hydrologic cycle, and so it does not calculate the ratio of bare ice to snow-covered ice, as previous studies (Shields et al., 2013) did. For a given stellar spetrum, we define total ice albedo as the sum , where and are the respective percent contributions of visible ( 700 nm) and near-infrared ( 700 nm) radiation. Following a previously published method (Kopparapu et al., 2013, 2014), we use âBT-Settlâ model spectra (Allard et al., 2007) to calculate the percent contribution of visible and near-infrared radiation for F-type (7200 K, 67% visible), G-type (5800 K, 52% visible), K-type (4600 K, 32% visible), and M-type (3400 K, 10% visible) stars.
We parameterize the infrared outgoing radiative flux as a fourth-order polynomial function of the partial pressure of carbon dioxide, CO, and surface temperature, . Likewise, we parameterize top-of-atmosphere albedo as a third-order polynomial function dependent on CO, , surface albedo , and stellar zenith angle. Assuming a noncondensible background pressure of 1 bar N, we use a radiative-convective climate model developed for HZ calculations (Kopparapu et al., 2013, 2014) to obtain best fits of more than 50,000 calculations for and over a parameter space spanning , , and , across all zenith angles. We calculate separate radiative transfer parameterizations for F, G, K, and M stars, using model spectra (Allard et al., 2007), which we describe further in the Appendix.
We assume in our EBM that any condensing CO accumulates on the surface as dry ice. CO condensation occurs when CO exceeds the CO saturation vapor pressure at surface temperature within a latitudinal zone. We assume that dry ice will radiatively dominate over water ice, so we set total ice albedo for frozen carbon dioxide (Warren et al., 1990) when CO condensation occurs. We also keep an inventory of the thickness of CO ice that condenses or melts on the surface at a given latitude, and we adjust the radiative contribution of CO at each latitude, as well as the global value of CO, by a corresponding amount each iteration. We calculate the thickness of accumulating ice as , where is the partial pressure of CO that condenses into ice, , and is the density of dry ice. Assuming that ice thickness is limited only by geothermal heat flow, we express the maximum CO ice thickness as: (Pollard & Kasting, 2005), where is the thermal conductivity of solid CO (Kravchenko & Krupskii, 1986; Stewart & Nimmo, 2002), is the temperature difference between the atmosphere and seawater beneath the ice, and is an Earth-like geothermal heath flux. For a temperature difference characteristic of a globally glaciated planet (Pollard & Kasting, 2005), this gives a maximum CO ice thickness of . Any additional accumulation would result in basal melting of CO glaciers and the transport of liquid CO to lower latitudes, although none of our simulations in this study reach conditions where .
Our latitudinal model can only represent clouds through adjustments to surface albedo, and we assume in our calculations that water clouds cover half of the surface. We have explored the sensitivity of climate to this cloud fraction parameter in a previous study (Fairén et al., 2012), which showed that excessive cloud cover can cause a planet to plummet into global glaciation. We account for the absorption of infrared radiation by clouds by subtracting a fixed amount of 8.5 W m from at each latitude band, following Williams & Kasting (1997). This value was selected by requiring that the EBM should produce a present-day Earth temperature of 288 K when the model is initialized with above-freezing initial conditions at . Possible sources of additional warming include CO ice clouds, which could warm the surface by up to 15 K by providing additional downward-directed infrared radiation through a scattering greenhouse effect (Forget & Pierrehumbert, 1997; Mischna et al., 2000; Colaprete & Toon, 2003; Forget et al., 2013). But this is not by itself enough to deglaciate a planet (Forget et al., 2013), and Kitzmann (2016) has shown that these previous calculations may have significantly overestimated the warming from such clouds. By neglecting their radiative impact, we generate a somewhat pessimistic outer limit on planetary habitability.
2.1 Outgassing Rates
CO accumulates in our model atmosphere as a result of the carbonate-silicate cycle, which allows a frozen planet to eventually deglaciate. We represent the 0.5 Myr timescale of the carbonate-silicate cycle with the time variable, , to contrast this slower timescale from the faster time step, , used in Eq. (1). Following Menou (2015) we implement the time evolution of CO into our EBM according to
where represents the volcanic outgassing rate of CO, represents the uptake of CO by rock weathering, and represents the uptake of CO by seafloor weathering. We begin by assuming volcanic outgassing to be a constant , where is the present-day value.
We also re-examine prior assumptions regarding the total rate of CO outgassing. Menou (2015) assumed a value of , which, converted to geochemists’ units for an Earth-mass planet, corresponds to rate of . Menou (2015) appears to have obtained this value from Abe et al. (2011), who in turn obtained this value from Sasal et al. (2002) as the expected CO flux from the depleted mantle alone. However, best estimates of the total terrestrial CO outgassing rate suggest (Gerlach, 2011; Jarrard, 2003), which is about a factor of ten greater than the value used by Menou (2015). By comparison, Tajika (2007) assumed a modern outgassing rate of 8 Tmol/yr. We argue that this modern rate is the best choice for assessing the habitability of a planet with Earth-like tectonic activity. We use a value of in most of our calculations, but we also explore the dependence of on defining the limit cycle HZ boundary.
2.2 Weathering Rates
where for , is an activation energy, and is a runoff efficiency factor. The value of is a critical factor determining the rate of CO evolution and has been estimated to be in the range of 0.25 to 1.0 in the absence of vascular plants (Kump et al., 2000; Berner & Kothavala, 2001; Abbot et al., 2012). A lower value of is appropriate for environmental conditions with pH 5, and probably represents a minimum dependence of weathering rate on CO (Berner, 1992). Tajika (2007) assumed , following Walker et al. (1981). Their parameterization, in turn, was based on laboratory studies of silicate dissolution by Lagache (1976). It remains unclear whether or not the presence of widespread biological systems should affect , although the presence of life should increase the weathering rate (Cawley et al., 1969; Schwartzman & Volk, 1989; Berner, 1992). We use a value of = 0.5 for most of our calculations, which is appropriate for conditions under which the weathering rate is proportional to dissolved [H]. We also perform a more limited set of calculations for .
Seafloor weathering has been neglected in previous calculations of limit cycling, but that is probably an oversight, particularly for planets with smaller amounts of land area than Earth. We consider the effects of seafloor weathering by assuming that the seafloor weathering rate is independent of temperature according to the functional form
where is a seafloor weathering parameter analogous to in Eq. (3) and is the baseline seafloor weathering that occurs in the absence of any CO dependence. Increases in the value of or cause a reduction in the frequency of limit cycling until a state of permanent glaciation occurs. This parameter may be near unity, but weaker weathering dependence with or lower may also be consistent with some computed histories of early Earth (Sleep & Zahnle, 2001). We initially assume in our calculations to assess the limit cycle boundary of the HZ, and we then consider sensitivity to changes in the seafloor weathering rate.
2.3 Partial pressure of CO in soil
A critical factor in determining the onset of limit cycles in the HZ is the partial pressure of CO in soil. The parameter represents the long-term balance between atmospheric and soil CO for present-day Earth. Menou (2015) assumed a value of to represent the pre-industrial CO level. However, the long-term balance between weathering and volcanism should be based on the value of CO in soil, rather than in the atmosphere. This implies that an abiotic Earth should have a higher value of atmospheric pCO than today. To put it another way, if all of life were to suddenly vanish, then pCO should increase until the atmospheric and soil (regolith) partiial pressures were equal. If we wish to use our model to test the habitability of abiotic planets, or of inhabited planets that lack vascular land plants, then we should tune the model to an abiotic state.
For a biotic planet like present-day Earth, root respiration by vascular plants increases the value of soil CO by a factor of 10 to 100 (Kump et al., 2010). We assume here that the enhancement is a factor of 30, in which case soil CO should be approximately bar. For , this implies that land plants accelerate silicate weathering by a factor of . For , the acceleration would be . In either case, an abiotic present-day Earth would be warmer than today because land plants would no longer be pumping atmospheric CO into soil. We therefore choose a value of for our calculations of the limit cycle HZ boundary.
We first consider an Earth-like (but abiotic) planet orbiting a G-type star like the Sun. At present-day stellar flux (), our weathering model for abiotic Earth with balances at soil CO = bar and average surface temperature of 296 K, while present Earth (with life) at a temperature of 288 K has a higher value of soil CO = bar (Fig. 1). By contrast, Menou (2015) argued that the carbon cycle on an abiotic Earth should balance at 288 K and CO = bar. This would only be true if land plants had zero effect on silicate weathering. Our model agrees with Menou (2015) in predicting no limit cycles for present-day Earth, but our results differ toward the outer edge of the HZ. At , our model predicts stable warm climates above the freezing point, whereas Menou (2015) predicts limit cycles with prolonged glacial conditions. Further outward in the HZ at (the effective solar flux at Mars’ orbit), our model and Menou (2015) both predict that the intersection between the weathering rate curve and the greenhouse effect curve falls beneath the freezing point, suggesting that limit cycles should occur. However, when we assume a value of in our model, limit cycles do not occur at all (Fig. 1, dashed green curve). Elsewhere (Batalha et al., submitted), we have argued that limit cycles are even more likely to have occurred on early Mars when the solar flux was significantly lower than today. For early Mars, though, the greenhouse effect must be supplemented by some absorber other than HO or CO; otherwise, even brief recoveries to above-freezing surface temperatures are impossible without outside stimuli, e.g., impact events (Segura et al., 2002).
The primary difference between the results of our model and those of Menou (2015) is caused by our higher assumed volcanic outgassing rate and by our treatment of the CO partial pressure of in soil. The relevance of these factors in determining the frequency of limit cycle events is shown in Fig. 2, which shows the limit cycle frequency (in units of cycles per Gyr) as a function of and . The grey region of this figure represents warm, stable solutions that are not prone to limit cycling. When we select and bar following Menou (2015), Fig. 2 predicts that limit cycles should occur with a frequency of about ten cycles per Gyr. Even if we assume values similar to those used by Kadoya & Tajika (2014), with and bar, this still results in climates prone to limit cycling. However, when we select our preferred values of as the present volcanic outgassing rate and bar as the present-day value of pCO in soil, then Fig. 2 predicts that limit cycles should not occur. Our improvements to the radiative transfer and our consideration of the mass balance of CO ice, which was absent in previous studies, make our model more accurate in predicting the frequency of limit cycle events; however, our assumptions about and determine when these limit cycles occur.
When we use our model to investigate the possibility of limit cycles occurring within the Sun’s HZ, we find that no such limit cycle boundary exists for our choices of volcanic outgassing rate () and soil CO (). This model configuration allows an Earth-like planet to deglaciate from a snowball state at any point within the conventional HZ due to the accumulation of a dense CO atmosphere from the carbonate-silicate cycle. This result may initially appear inconsistent with Fig. 1, where we predict limit cycles for . However, the climate calculations in Fig. 1 are global averages from a one-dimensional model (Kopparapu et al., 2013, 2014), in which the surface must be either completely ice-free or completely ice-covered. By contrast, our calculations with a latitudinal EBM in Fig. 3 allow polar ice caps to form while still retaining ice-free conditions at lower latitudes, which permits stable climates to persist across the entire HZ. This illustrates the importance of using latitudinally-resolved models for situations where a planet’s surface is partly ice-covered but remains otherwise habitable. Thus, we still retain a classic picture of the HZ, where warm stable climates (even for abiotic planets) are possible all the way out to the maximum greenhouse effect. This lack of limit cycling throughout the HZ applies to planets orbiting F, G, K, and M stars, as the tuning of our model to soil CO results in warm stable climates regardless of stellar type. A decrease in would cause the weathering rate to be even less sensitive to changes in CO, which only expands the climatically stable parameter space where limit cycles do not occur. This suggests that the effect of limit cycling for planets orbiting near the outer edge of the HZ may have been overestimated by previous studies (Kadoya & Tajika, 2014, 2015; Menou, 2015).
We have argued that volcanic outgassing should balance weathering at the soil CO, rather than the atmospheric CO, which results in no limit cycles for Earth-like planets in the HZ. However, this result also depends on our assumed volcanic outgassing rate, which we have tuned to present-day Earth. Kadoya & Tajika (2014, 2015) have argued convincingly that the habitability of an Earth-like planet is highly sensitive to the CO degassing rate. (See also earlier papers by Tajika (2003, 2007).) Lower rates of volcanism will decrease atmospheric CO, which will increase the susceptibility to limit cycles, as we show in Fig. 2. Furthermore, increased fractional land area could accelerate loss of CO by silicate weathering and make a planet more prone to limit cycles, provided that sufficient rainfall is available. Thus, limit cycles could occur on planets that have a reduced volcanic outgassing rate or increased weathering rate compared to abiotic Earth.
We next consider the effect of volcanic outgassing rate on the occurrence of limit cycles in the HZ. Beginning with a planet orbiting a G-type star, we decrease the volcanic outgassing rate to , and we calculate the distance at which limit cycles begin by decreasing the relative solar flux until cycles start to develop. (Solar flux and orbital distance are related by the inverse square law.) This transition occurs at = 0.77, which corresponds to an orbital distance of about 1.2 AU. At that distance and beyond, the planet’s climate exhibits warm intervals of 5 Myr in duration with equatorial temperatures above the freezing point of water, separated by extended periods of global glaciation lasting 60 Myr (Fig. 3, top panel). During the warm intervals, increased weathering rates lead to a decrease in atmospheric CO, which eventually triggers the subsequent glaciation. Because continental weathering ceases entirely during glaciation, atmospheric CO accumulates from ongoing volcanic outgassing.
The effects of limit cycling are quite different for hot, blue, F-type stars than for cooler, red, K- and M-type stars, largely because of the way in which incident stellar radiation interacts with surface water ice. Water ice is highly reflective at visible wavelengths (700 nm), but becomes an increasingly efficient absorber at longer, near-infrared wavelengths. Ice-albedo feedback is therefore greater for planets around F stars than it is for planets around K and M stars because F stars emit a greater percentage of their radiation at visible wavelengths (Joshi & Haberle, 2012; Shields et al., 2013). For F stars, again assuming , the limit cycle transition occurs at = 1.07, which corresponds to an orbital distance of about 0.96 AU. This planet experiences warm periods lasting 1 Myr with prolonged glacial periods of 20 Myr (Fig. 3, bottom panel). Planets orbiting K- and M-type stars do not experience limit cycles, even at this reduced volcanic outgassing rate, but instead these systems can maintain stable warm conditions all the way out to the maximum greenhouse limit.
Seafloor weathering (Coogan & Gillis, 2013) is another potential sink for atmospheric CO, which can make a planet more susceptible to limit cycles or permanent glaciation. For a G star planet with = 0.7 and that is already experiencing limit cycles, an increase in the rate of seafloor weathering serves to decrease the frequency in cycles between climate states, which results in longer extended periods of glaciation between warm episodes (Fig. 4). The rate of seafloor weathering may also be sensitive to atmospheric CO through the parameter in Eq. (4), but increases in only accentuate the tendency of a planet in a limit cycle toward permanent glaciation. For planets beyond the limit cycle region, an increase in seafloor weathering will cause the limit cycle region of the HZ to expand, analogous to a decrease in volcanic outgassing of the same magnitude. The rest of our calculations will consider sensitivity to the volcanic outgassing rate, , but these results apply more broadly to the difference between volcanic outgassing and seafloor weathering, , that drives long-term changes in CO (Eq. (3)). Strong seafloor weathering could cause a planet to be prone to limit cycles even if the volcanic outgassing rate is at or above present-day values.
We summarize our findings by performing similar calculations of the limit cycle boundary for F, G, K, and M stars with and combining these calculations with the results for in Fig. 5. As discussed above, planets with a present-day volcanic outgassing rate avoid limit cycling altogether, while a rate of or lower can be sufficient to induce cycling for at least part of the HZ. Planets orbiting F stars exhibit cycling behavior throughout a large region of their HZs and are the most responsive to changes in , as a result of their increased sensitivity to ice-albedo feedback. A planet like Earth around a G star should have experienced limit cycles in its past only if volcanic outgassing were much lower than today. Modern Earth avoids this fate, but early Earth may have experienced repeated Snowball Earth episodes if volcanic outgassing rates were relatively low (Sleep & Zahnle, 2001; Tajika, 2007). The rock record is sparse or nonexistent during the first half of Earth’s history, so such behavior could have occurred but remained undetected.
By contrast, late K and M star planets can maintain stable climates without cycling until the maximum greenhouse limit is reached. Such planets are subject to tidal locking and may show only one face to the star, as the Moon does to Earth, which may or may not present a problem for complex life. However, planets orbiting these stars may lose their entire water inventory as a consequence of a runaway greenhouse during the extended, hot pre-main sequence evolution of the host star (Ramirez & Kaltenegger, 2014; Luger & Barnes, 2015; Tian & Ida, 2015). This suggests that such planets may be entirely uninhabitable unless they begin with an abundant water inventory or water is somehow resupplied after the star enters the main sequence.
Note that the HZ itself is also narrower in Fig. 5 than prior estimates because the requirement for a planet to be able to recover from global glaciation shifts the outer edge slightly inward in our model. This last result is sensitive to the fraction of the planet’s surface covered by (lower albedo) bare ice, as opposed to (higher albedo) snow-covered ice. In our model the ice is snow-covered at all latitudes, and so its albedo remains relatively high, even in the near-infrared. By contrast, Shields et al. (2013) found that the outer edge should be completely insensitive to surface albedo because its effect would be completely masked by the overlying dense CO-rich atmosphere. However, such dense CO atmospheres may not be long-lived when carbonate-silicate weathering is included.
The Rare Earth hypothesis (Ward & Brownlee, 1999) suggests that complex life may be uncommon in the universe even if simple life is widespread. Although many arguments have been raised against this idea (Kasting, 2001), limit cycling in parts of the conventional HZ, for planets with relatively low outgassing rates, presents problems for both simple and complex life. Photosynthetic algae and cyanobacteria would go extinct on a ‘hard Snowball’ planet with sea-ice thicknesses of a kilometer or more, unless illuminated refugia were available. Both types of organisms survived the Neoproterozoic Snowball Earth episodes (Hoffman et al., 1998), so some types of refugia must have existed. Models with either thin ice (Pollard & Kasting, 2005) or open water (Abbot et al., 2011) near the equator may provide the explanation. But multicellular land life would be highly challenged by this type of climatic behavior which, fortunately, has not occurred since the late Precambrian. It follows that animal life, and thus intelligent life, may not be able to evolve on planets with low incident stellar flux and a low volcanic outgassing rate, even if they are within the conventional HZ.
The presence of a limit cycle boundary depends critically upon the assumed volcanic outgassing rate, and planets with a CO outgassing rate similar to today may not experience limit cycles at all. Kadoya & Tajika (2015) suggest that low outgassing rates may be expected for Earth-like planets orbiting old stars, which would make such planets more prone to limit cycles. One should exercise caution in accepting this conclusion, as it is based on their assumption that the CO outgassing rate declines monotonically with time on an Earth-like planet as its interior cools. Other authors (Holland, 2009) have suggested that Earth’s CO outgassing rate actually increased with time during the first half of its history as the growth of continents allowed greater storage of carbonate rocks and greater recycling of CO through weathering, carbonate deposition, and sediment subduction.
If the CO outgassing rate of present-day Earth is anomalously large compared to typical terrestrial planets, then Earth might be uncommon in its ability to sustain a stable warm climate. Conversely, planets more massive than Earth (known as super-Earths) may exhibit higher rates of volcanism than Earth today, although the dependence of plate tectonics on planetary mass remains unclear (Valencia et al., 2007; Kite et al., 2009; Korenaga, 2010; van Heck & Tackley, 2011; Haghighipour, 2013). So, it is at least conceivable that super-Earths in the outer parts of the HZ would be better abodes for complex life than would true Earth analogs.
The actual dependence of the silicate weathering rate on CO is unknown for abiotic planets. We have assumed , which matches the behavior of the H concentration in rainwater; however, the actual exponent could range from 0.25 to 1 (Berner, 1994). The weathering rate on modern Earth is sometimes argued to have zero direct dependence on atmospheric CO (Berner et al., 1983), because CO in soils is decoupled from atmospheric CO by the presence of vascular plants, which pump up soil CO by way of root respiration. Plants also generate humic acids which accelerate weathering, again without any direct relation to atmospheric CO. Consequently, Menou (2015) argued that the emergence of land life on a planet should stabilize its climate against limit cycles. But this inference is incorrect. Land plants accelerate weathering (Berner, 1992) by anywhere from a factor of 2-3 (Cawley et al., 1969) to a factor of 10-100 or more (Schwartzman & Volk, 1989) compared to an abiotic environment. (We assume a weathering acceleration factor in the range of 2.8-5.5, as dicussed above.) Lichens, algae, and other microorganisms also secrete acids that accelerate weathering (Berner, 1992), so the emergence of life on a planet in the outer, limit-cycling region of the HZ should only help to pull down atmospheric CO, making the planet even more subject to global glaciation. Global glaciation would kill any plants, allowing atmospheric CO to again accumulate, and so cycling should re-initiate at a rate that would depend on whether the plants themselves were able to regenerate. Life (as we know it) would not stabilize a planet’s climate against limit cycling, but it might create a more complex, biologically mediated form of limit cycling.
What do these arguments imply about the prevalence of animal life and the possible evolution of intelligent life? The limit cycle region of the HZ depends upon the assumed, and somewhat uncertain, behavior of the volcanic outgassing and seafloor weathering rates for abiotic planets. For planets around K and M stars, this does not appear to pose a problem, but for some types of stars, the outlook is less optimistic. According to our results in Fig. 5, the HZ around F stars can be nearly eliminated if the volcanic outgassing rate is or less. F stars also have relatively short main sequence lifetimes, and they brighten quickly as they age, which limits the time available for biological evolution. Meanwhile, although planets around K and M stars may avoid this problem, they may have other issues that could preclude their habitability. If Ramirez & Kaltenegger (2014), Luger & Barnes (2015), and Tian & Ida (2015) are correct, nearly all planets around late K and M stars should experience drastic pre-main-sequence water loss. Early K stars and late G stars avoid both of these problems. So, there are still many stars that could host planets with complex life. But any search for such life should be concentrated on planets around late G- to early K-type stars, which are only a subset of the planets that might support simple life.
We have calculated the limit cycle boundary of the habitable zone as a function of stellar type and CO outgassing rate. Earth-like planets with volcanic outgassing rates similar to today are able to maintain stable climates across the entire range of the HZ, regardless of stellar type. But planets with lower volcanic outgassing rates or significant seafloor weathering rates should experience limit cycles, with punctuated episodes of warm conditions followed by extended glacial periods. F star planets are the most prone to this behavior as a result of increased susceptibility to ice-albedo feedback. Planets orbiting late K and M stars avoid limit cycles because of reduced ice-albedo feedback, but they may suffer from water loss during their formation. Thus, systems with the greatest potential for habitability are those around late G and early K type stars.
The net outgassing rate of CO and partial pressure of in soil are key parameters in understanding the habitability of an Earth-like planet. If Earth has maintained a net CO outgassing rate at or above its current value for its entire history, then it may never have been prone to limit cycles at any point in time. By extension, if Earth’s outgassing rate is typical for other terrestrial planets, then limit cycles may not pose a problem for habitability at all. Likewise, we expect that an abiotic planet would accumulate more atmospheric CO than its inhabited counterpart, which would lead us to expect that few planets, if any, should reside in limit cycles. However, we should expect a diversity of exoplanet environments to exist, and we cannot rule out the possibility that Earth’s CO inventory is atypical. In the search for habitable worlds, we should at least consider the possibility that some Earth-like planets may exhibit lower net outgassing rates than today, which could develop limit cycles and preclude the development of complex life.
Appendix A Polynomial fits to OLR and planetary albedo for F, G, K, and M stars
We parameterized top of the atmosphere (TOA) albedo, , and the outgoing IR flux, , as polynomials with the following variables: surface temperature (K) used as ; CO, where CO is the partial pressure of CO (bars); where is the zenith angle; and surface albedo . The parameterizations were derived by running a 1-D radiative convective (RC) model (Kopparapu et al., 2013, 2014) over a range of values of the above parameters with a 1 bar N noncondensible background for each stellar type. The fits are valid in the range K K, bar CO bar, and . For planetary albedo, we made separate fits above and below 250 K.
Solar zenith angle is explicitly calculated at each latitude in the EBM as a function of solar declination and solar hour angle, averaging over a complete rotation to obtain the insolation-weighted zenith angle (Williams & Kasting, 1997; Cronin, 2014). In general, we expect that TOA should increase as a function of . The usual configuration of our 1-D RC model (Kasting, 1991; Kasting et al., 1993; Kopparapu et al., 2013, 2014) assumes that stratospheric temperature is equal to the ‘skin temperature’ of a gray atmosphere (i.e., the temperature at an optical depth of zero), given by
where is the Stefan-Boltzmann constant. However, Eq. (A1) is only appropriate for global average conditions (i.e., for ) and cannot be applied to calculate and for the range of solar zenith angles across all latitude bands represented in our EBM. In order to circumvent this problem, we followed Williams & Kasting (1997) and parameterized as a function of :
where is the absorbed fraction of incident solar flux, which was calculated for a variety of zenith angles between and using the radiative-convective model. The value of () was obtained using Eq. (A1) above. This configuration of our RC model allows the parameterizations of TOA below to accurately decrease with as expected.
a.1 Parametric expressions for OLR and Planetary Albedo
The coefficients are provided as downloadable ASCII data tables in online supplementary information. The rows in the data table correspond to the order of the coefficients in the following expressions. The OLR coefficients are the same for all different stars, as it is a planet specific property and variation in stellar type does not change it. For planetary albedo coefficients, each stellar type has two sets of data tables, one for surface temperatures between K to K and another for K to K.
These paramtetric fits provide a rapid means of obtaining OLR and planetary albedo from our 1-D RC climate model calculations. The error between the climate model data and the parametric fits for OLR does not exceed 2%, and the error distributions of planetary albedo for different stellar spectral types show that a majority of parametric fits have less than 20% uncertainty. Error plots for OLR and planetary albedo are provided as online supplementary information.
- Abbot et al. (2011) Abbot, D., Voigt, A., & Koll, D. 2011, J. Geophys. Res., 116, D18103
- Abbot et al. (2012) Abbot, D. S., Cowan, N. B., & Ciesla, F. J. 2012, ApJ, 756, 178
- Abe et al. (2011) Abe, Y., Abe-Ouchi, A., Sleep, N. H., & Zahnle, K. J. 2011, Astrobiol., 11, 443
- Allard et al. (2007) Allard, F., Allard, N. F., Homeier, D., et al., 2007, A&A, 474, L21-L24
- Batalha et al. (submitted) Batalha, N. E., Kopparapu, R. K., Haqq-Misra, J., & Kasting, J. F., submitted, Earth Planet. Sci. Lett.
- Berner et al. (1983) Berner, R. A., Lasaga, A. C., & Garrels, R. M. 1983, Am. J. Sci., 283, 641
- Berner (1992) Berner, R. A., 1992, Geochim. Cosmochim. Acta, 56, 3225
- Berner (1994) Berner, R. A., 1994, Am. J. Sci., 294, 56
- Berner & Kothavala (2001) Berner, R. A. & Kothavala, Z., 2001, Am. J. Sci., 301, 182
- Cawley et al. (1969) Cawley, J. L., Burruss, R. C., Holland, H. D., 1969, Science, 165, 391
- Colaprete & Toon (2003) Colaprete, A. & Toon, O. B., 2003, J. Geophys. Res., 108, 5025
- Coogan & Gillis (2013) Coogan, L. A. & Gillis, K. M., 2013, Geochem. Geophys. Geosys., 14, 1771
- Cronin (2014) Cronin, T. W., 2014, J. Atmos. Sci., 71, 2994
- Fairén et al. (2012) Fairén, A. G., Haqq-Misra, J. D., & McKay, C. P. 2012, A&A, 540, A13
- Forgan (2013) Forgan, D., 2013, MNRAS, 437, 1352
- Forget & Pierrehumbert (1997) Forget, F. & Pierrehumbert, R., 1997, Science, 278, 1273
- Forget et al. (2013) Forget, F., Wordsworth, R., Millour, E., et al., 2013, Icarus, 222, 81
- Gaidos & Williams (2004) Gaidos, E. & Williams, D. M., 2004, New Astron., 10, 67
- Gerlach (2011) Gerlach, T., 2011, Eos, 92, 201
- Haghighipour (2013) Haghighipour, N., 2013, Annu. Rev. Earth Planet. Sci., 41, 469
- Haqq-Misra (2014) Haqq-Misra, J., 2014, J. Adv. Model. Earth Syst., 6, 950
- Holland (2009) Holland, H. D., 2009, Geochim. Cosmochim. Acta 73, 5241
- Hoffman et al. (1998) Hoffman, P. F., Kaufman, A. J., Kalverson, G. P., & Schrag, D. P., 1998, Science, 281, 1342
- Jarrard (2003) Jarrard, R. D., 2003, Geochem. Geophys. Geosyst., 4, 8905
- Joshi & Haberle (2012) Joshi, M. M. & Haberle, R. M., 2012, Astrobiol., 12, 3
- Kadoya & Tajika (2014) Kadoya, S. & Tajika, E., 2014, ApJ, 790, 107
- Kadoya & Tajika (2015) Kadoya, S. & Tajika, E., 2015, ApJ, 815, L7
- Kaltenegger & Sasselov (2011) Kaltenegger, L. & Sasselov, D., 2011, ApJ, 736, L25
- Kane & Gelino (2012) Kane, S. R. & Gelino, D. M., 2012, Astrobiol., 12, 940
- Kane & Hinkel (2013) Kane, S. R. & Hinkel, N. R., 2013, ApJ, 762, 7
- Kasting (1991) Kasting, J., 1991, Icarus, 94, 1
- Kasting (2001) Kasting, J., 2001, Persp. Biol. Med., 44, 117
- Kasting et al. (1993) Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108
- Kite et al. (2009) Kite, E. S., Manga, M., & Gaidos, E., ApJ, 700, 1732
- Kitzmann (2016) Kitzmann, D., 2016, ApJ, 817, L18
- Kopparapu et al. (2013) Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al., 2013, ApJ, 765, 131
- Kopparapu et al. (2014) Kopparapu, R. K., Ramirez, R., SchottelKotte, J., et al., 2014, ApJ, 787, L29
- Korenaga (2010) Korenaga, J., 2010, ApJ, 725, L43
- Kravchenko & Krupskii (1986) Kravchenko, Y. & Krupskii, I., 1986, Sov. J. Low Temp. Phys, 12, 46
- Kump et al. (2000) Kump, L. R., Brantley, S. L., & Arthur, M. A. 2000, Annu. Rev. Earth Planet. Sci., 28, 611
- Kump et al. (2010) Kump, L. R., Kasting, J. F., & Crane, R. G. 2010, The Earth System, Upper Saddle River, NJ
- Lagache (1976) Lagache, M., 1976, Geochim. Cosmochim. Acta, 40, 157
- Leconte et al. (2013) Leconte, J., Forget, F., Charnay, B., et al., 2013, A&A, 554, A69
- Luger & Barnes (2015) Luger, R. & Barnes, R., 2015, Astrobiol., 15, 119
- Menou (2015) Menou, K., 2015, Earth Planet. Sci. Lett., 429, 20
- Mills et al. (2011) Mills, B., Watson, A. J., Goldblatt, C., et al., 2011, Nat. Geosci., 4, 861
- Mischna et al. (2000) Mischna, M. A., Kasting, J. F., Pavlov, A., & Freedman, R., 2000, Icarus, 145, 546
- North et al. (1981) North, G. R., Cahalan, R. F., & Coakley, J. A. 1981, Rev. Geophys., 19, 91
- Pavlov et al. (2003) Pavlov, A. A., Hurtgen, M. T., Kasting, J. F., & Arthur, M. A., 2003, Geology, 31, 87
- Pollard & Kasting (2005) Pollard, D. & Kasting, J. F., 2005, J. Geophys. Res., 110, C07010
- Ramirez & Kaltenegger (2014) Ramirez, R. & Kaltenegger, L., 2014, ApJ, 797, L25
- Rothman et al. (2009) Rothman, L. S., Gordon, I. E., Barbe, A., et al., 2009, J. Quant. Spectrosc. Radiat. Transfer, 110, 533
- Rothman et al. (2010) Rothman, L. S., Gordon, I. E., Barbe, A., et al., 2010, J. Quant. Spectrosc. Radiat. Transfer, 111, 2139
- Sasal et al. (2002) Saal, A. E., Hauri, E. H., Langmuir, C. H, & Perfit, M. R., 2002, Nature, 419, 419
- Schwartzman & Volk (1989) Schwartzman, D. W. & Volk, T., 1989, Nature, 340, 457
- Segura et al. (2002) Segura, T. L., Toon, O. B., Colaprete, A., Zahnle, K., 2002, Science 298, 1977
- Shields et al. (2013) Shields, A. L., Meadows, V. S., Bitz, C. M., et al., 2013, Astrobiol., 13, 715
- Showman et al. (2013) Showman, A. P., Wordsworth, R. D., Merlis, T. M., et al., 2013, Comparative Climatology of Terrestrial Planets, 277
- Sleep & Zahnle (2001) Sleep, N. H. & Zahnle, K., 2001, J. Geophys. Res., 106, 1373
- Stewart & Nimmo (2002) Stewart, S. T. & Nimmo, F., 2002, J. Geophys. Res., 107, 5069
- Tajika (2003) Tajika, E., 2003, Earth Planet. Sci. Lett., 214, 443
- Tajika (2007) Tajika, E., 2007, Earth Planet. Space, 59, 293
- Tian & Ida (2015) Tian, F. & Ida, S., 2015, Nat. Geosci., 8, 177
- Valencia et al. (2007) Valencia, D., O’Connell, R. J., & Sasselov, D. D., ApJ, 670, L45
- van Heck & Tackley (2011) van Heck, H. J. & Tackley, P. J., 2011, Earth Planet. Sci. Lett., 310, 252
- Vladilo et al. (2013) Vladilo, G., Murante, G., Silva, L., et al., 2013, ApJ, 767, 65
- von Braun et al. (2011) von Braun, K., Boyajian, T. S., Kane, S. R., et al., 2011, ApJ, 729, L26
- Walker et al. (1981) Walker, J. C. G., Hays, P. B., & Kasting, J. F., J. Geophys. Res., 86, 9776
- Ward & Brownlee (1999) Ward, P. D. & Brownlee, D., 1999, Rare Earth, Springer-Verlag, New York
- Warren et al. (1990) Warren, S. G., Wiscombe, W. J., & Firestone, J. F. 1999, J. Geophys. Res., 95, 14717
- Williams & Halloway (1982) Williams, G. P. & Halloway, J., 1982, Nature, 297, 295
- Williams & Kasting (1997) Williams, D. M. & Kasting, J. F., 1997, Icarus, 129, 254
- Wolf & Toon (2014) Wolf, E. T. & Toon, O. B., 2014, Geophys. Res. Lett., 41, 167
- Yang et al. (2014) Yang, J., Boué, G., Fabrycky, D. C., & Abbot, D. S., 2014, ApJ, 787, L2