GCM Simulations of Unstable Climates in the Habitable Zone

TL;DR

This study confirms that Earth-like planets in the outer habitable zone can experience unstable climate cycles, using advanced 3D climate models and a novel hybrid integrator to validate previous low-dimensional model results.

Contribution

It introduces a hybrid 0D-3D climate modeling approach to verify climate cycle phenomenology in complex GCM simulations, extending prior idealized models.

Findings

GCM climate cycles align with low-dimensional model results

Outer habitable zone planets may often be in transient glaciated states

Hybrid integrator effectively simulates long-term climate evolution

Abstract

It has recently been proposed that Earth-like planets in the outer regions of the habitable zone experience unstable climates, repeatedly cycling between glaciated and deglaciated climatic states (Menou 2015). While this result has been confirmed and also extended to explain early Mars climate records (Haqq-Misra et al. 2016; Batalha et al. 2016), all existing work relies on highly idealized low-dimensional climate models. Here, we confirm that the phenomenology of climate cycles remains in 3D Earth climate models with considerably more degrees of freedom. To circumvent the computational barrier of integrating climate on Gyr timescales, we design a hybrid 0D-3D integrator which uses a general circulation model (GCM) as a short relaxation step along a long evolutionary climate sequence. We find that GCM climate cycles are qualitatively consistent with reported low-dimensional results. This establishes on a firmer ground the notion that outer habitable zone planets may be preferentially found in transiently glaciated states.

Equations

$$\frac{W}{W_{\oplus }}={\left(\frac{{\text{pCO}}_2}{{\text{pCO}}_{2,\oplus }}\right)}^{\beta }{e}^{[k_{\text{act}}(T_{\text{s}}-288)]}[1+k_{\text{run}}(T_{\text{s}}-288){]}^{0.65},$$

$$\frac{d}{dt}pC{O}_2=V-W(T_{\text{s}},pC{O}_2),$$

$$\Delta pC{O}_2=\text{sign}(V-W)\Delta T\frac{\partial pC{O}_2}{\partial T},$$

$$\Delta t=\frac{V-W}{\Delta pC{O}_2}.$$

$$\frac{V}{V_{\oplus }}=\kappa e^{-\alpha \left(\frac{F}{F_{\oplus }}\right)},$$

$$\frac{W_{sf}}{W_{\oplus }}=\frac{{\varphi }^n+\nu }{1+\nu }$$

$$W={\beta }_0{W}_{\text{c}}+(1-{\beta }_0)W_{sf},$$

Full text

planets and satellites: terrestrial planets, planets and satellites: atmospheres, methods: numerical, astrobiology

††journal: The Astrophysical Journal††software: PlaSim (Fraedrich et al., 2005), matplotlib (Hunter, 2007)

\censorruleheight

=.1ex

I Introduction

Recent estimates of the habitable zone for Earth-like planets with modern Earth’s rotation rate place an inner edge at 0.99 AU and an outer edge at 1.7 AU (Kopparapu et al., 2013). However, these estimates come from 1D climate models which only consider equilibria on short timescales. That short-term equilibrium comes from an energy balance requirement—the heat released into space must balance the net solar radiation (insolation) received from the sun North et al. (1981). We refer to this as thermal equilibrium. Planets can also demonstrate stability or instability on longer timescales. Changes in glaciation or greenhouse gases over time can alter the planet’s albedo and outgoing thermal emission, affecting the thermal equilibrium temperature (North et al., 1981; Pierrehumbert, 2005). It may be that a more realistic habitable zone that takes into account long-term stability will be significantly more constrained. The emergence of complex life depends as much on the long-term stability of a habitable climate as it does on its immediate habitability—an informed search for Earth analogs thus depends strongly on understanding climate variation throughout the habitable zone.

A planet’s climate stability can be tied to the carbon-silicate cycle. Atmospheric carbon dioxide dissolves into rainwater as carbonic acid, is delivered to the surface, reacts with the rock, and is eventually carried to the oceans, where it can be recycled into the mantle via subduction (Berner, 2004; Pierrehumbert, 2010). On planets containing vascular land plants, the direct dependence of weathering on the CO2 partial pressure (hereafter pCO2) is weakened by plant activity affecting carbon levels in the soil (Berner, 2004; Pierrehumbert, 2010; Schwartzman & Volk, 1989). However, on planets lacking land plants, the carbon cycle is more strongly-coupled and acts as a negative feedback on the planet’s surface temperature—increased temperatures or pCO2 surpluses lead to increased weathering and CO2 sequestration, while a deficit reduces weathering, allowing volcanic outgassing to replenish CO2 (Pierrehumbert, 2010; Walker et al., 1981; Schwartzman & Volk, 1989; Williams & Kasting, 1997; Kump et al., 2000). The climate on such a planet tends toward an equilibrium point. At low insolation, the equilibrium point may be cold enough for total glaciation (Menou, 2015; Pierrehumbert, 2005; Hoffman et al., 1998), called a “snowball” state. Budyko (1969) showed that such states can exist at the same insolations and CO2 partial pressures as temperate climates, and therefore require either high insolation or significant greenhouse gases for deglaciation. In this scenario of a snowball state triggered by weathering, weathering ceases in the absence of liquid water, implying that no weathering equilibrium can be reached. Volcanic outgassing continues, and at a certain point pCO2 will be high enough for the ice to melt, resulting briefly in warm wet conditions, followed by rapid weathering and cooling. Thus, on geological timescales the climate may cycle between long snowball events and brief warm wet periods (Menou, 2015). Evidence of glacial deposits at then-tropical latitudes on Earth has been linked to this hypothesis (Tajika, 2007), and previous studies (e.g. Pierrehumbert, 2005; Lucarini et al., 2010; Voigt et al., 2011) have shown that such climates are not ruled out by 3D climate models. We do not however aim to determine whether such episodes on Proterozoic Earth were triggered by weathering and later deglaciated by this process as opposed to other external triggers; there are sufficient uncertainties in both our model and Earth’s conditions at the time that more sophisticated tools would be necessary.

Menou (2015) used a zero-dimensional reduction of 1D models to approximate energy-balanced climates of Earth-like planets at varying orbital distances throughout the habitable zone. He found that climate cycles, wherein the planet oscillates between brief warm periods and long frozen periods, should dominate much of the outer habitable zone. However, energy-balance models (EBMs) necessarily miss or must approximate inherently 3D phenomena and feedbacks, and so are of limited utility for predicting new climate phenomena. In the case of climate cycles, glaciation and deglaciation are dependent on factors such as sea ice distribution, precipitation patterns, heat transport by wind streams, local heating and cooling by clouds, seasonal melting, etc. Stronger or weaker regional temperature gradients can mean that the onset of snowball glaciation happens at global average surface temperatures that deviate significantly from 273.15 K. Since weathering is dependent on surface temperature, cooling or warming of a region due to changes in precipitation, cloud activity, and air streams as the climate changes can have a significant on weathering. It is therefore necessary to verify the existence of climate cycles using more sophisticated 3D climate models that can capture a diversity of climate feedbacks.

II Methods

II.A Simulating Long Timescales

Since climate cycles, as observed in 1D models by Menou (2015) and Haqq-Misra et al. (2016), have periods of at minimum millions of years, the primary obstacle to simulating them in 3D models is the computational cost of millions of years or more of climate evolution in existing 3D models. PlaSim, the GCM used in this study, takes on the order of minutes of wall-time to compute one simulated year. Computing millions of years directly is thus unfeasible. A more detailed discussion of the computational cost of direct computation vs. our method is presented in Appendix A. An approach is needed that allows a model to skip through millions of years of evolution, only computing a couple years at certain key points.

We achieve this by noting that the climate exhibits different types of equilibria, which operate on different timescales. Thermal equilibrium, when the planet exhibits a relatively unchanging seasonally-averaged global temperature, is achieved on the order of years to decades following a perturbation. Weathering equilibrium, on the other hand, when global outgassing and weathering rates cancel each other out and the global pCO2 level is relatively constant, is achieved (if at all) on much longer geological timescales of millions of years. This latter timescale is many orders of magnitude longer than the thermal timescale, which means that weathering effects can typically be disregarded on thermal timescales, and vice versa. The CO2 partial pressure is a global quantity, which means that only global, zero-dimensional terms are needed to compute its evolution. Thermal equilibrium, however, must be computed using a full 3D treatment to capture various mechanisms such as sea ice feedback, wind circulation, cloud formation, etc. The global weathering rate also depends on 3D features such as landmasses and ice caps. Therefore, we construct a hybrid 3D-0D model in which a full 3D model computes the evolution of a new thermal equilibrium following a 0D perturbation in pCO2. We compute the size and duration of that perturbation by computing the global weathering rate from the results of the 3D model. This means that millions of years of evolution of slow change can instead be approximated by allowing the CO2 partial pressure to change by an equivalent amount, and then allowing the climate to relax back to thermal equilibrium. We neglect the intermediate timescale on which the atmospheric and oceanic carbon reservoirs equilibriate, assuming that partitioning is instantaneous. This reduces the computational cost, at the expense of model accuracy. In this way, we treat the GCM not as the primary driver of model evolution, but as a relaxation step in a simpler 0D model that controls the long-term model evolution.

This is very similar to the timescale-coupling approach used in Donnadieu et al. (2006) and Le Hir et al. (2009) in their studies of the Pangea breakup and the Snowball Earth episode’s aftermath, respectively. However, their approach involved sampling the expected evolutionary track and then interpolating between points using lookup tables to compute the time-dependence of the climate’s evolution. Because our aim is to determine the very nature of the evolutionary track and to establish its time-dependence, we take a slightly different approach by extrapolating forward in time without a priori knowledge of what the future climate states will look like, computing timescales using evaluated rather than interpolated weathering fluxes to determine timescales.

For the 3D relaxation step, we use the Planet Simulator (PlaSim) developed at the Universität Hamburg, a 3D Earth-system general circulation model (GCM) developed for fast and easy study of Earthlike climate systems. PlaSim solves primitive equations for vorticity, temperature, divergence, and pressure using the spectral transform method. It includes a 50-meter mixed-layer slab ocean, thermodynamic sea ice, a bucket model of soil hydrology, moist and dry convection, interactive clouds, large-scale precipitation, and a parameterized longwave and shortwave radiation model (Fraedrich et al., 2005). It uses a latitude-longitude grid with an atmosphere discretized into multiple vertical levels. For a more complete description of the model, readers are referred to the PlaSim Reference Manual (Lunkeit et al., 2012). PlaSim has been used before to study snowball climates and to map climate parameter spaces of Earth-like planets (Lucarini et al., 2010; Boschi et al., 2013; Lucarini et al., 2013). PlaSim succeeds at reproducing the snowball hysteresis described in Budyko (1969) and Pierrehumbert (2005) (Lucarini et al., 2010), and yields results consistent with the findings in Pierrehumbert (2005) that snowball climates have greater seasonal variability, particularly in the troposphere lapse rate and tropopause height, and demonstrate on average weaker lapse rates and a lower tropopause. In PlaSim, snowball deglaciation is associated with the onset of melting sea ice in the tropics. We note that Abbot et al. (2010) has called this deglaciation mechanism into question as a possible artifact of simplistic sea ice schemes.

PlaSim however is a model of intermediate complexity, and while it qualitatively reproduces the snowball climates described by Pierrehumbert (2005) and others, it necessarily makes approximations, parameterizations, and omissions as a tradeoff for a large advantage in speed. The aim of this paper therefore is not to thoroughly investigate the dynamics of a snowball climate or determine precise quantitative limits on glaciation and deglaciation, as slower, more sophisticated models would be better-suited to that task, but rather to take advantage of PlaSim’s speed and simulate long timescales across wide swaths of parameter space, to verify whether or not climate cycles persist when investigated with models of increasing complexity and to identify overall trends in the habitable zone parameter space.

Our weathering implementation is based on a parameterization of Earth’s silicate weathering in the absence of vascular land plants. We use the same form used in Menou (2015), which is a form of the parameterization introduced in Berner (1994) and expanded on in Berner (2001). That parameterization is based on observations of cation concentrations in river runoff, laboratory measurements, and seafloor drill cores. Our functional weathering form is thus

[TABLE]

where $T_{s}$ is the surface temperature, pCO2,⊕ is the pre-industrial pCO2 of 330 $\mu$bar, and $W_{\oplus}$ is the weathering rate at $T_{\text{s}}$ = 288 K. $k_{\text{act}}$ is related to the chemical activation energy of the weathering reaction and is set to 0.09, and $k_{\text{run}}$ is a runoff efficiency factor set to 0.045, which helps account for changes in precipitation. $\beta$ is the dependence on pCO2, and is set to 0.5 (Menou, 2015; Pierrehumbert, 2010; Kump et al., 2000). We follow a similar procedure to Edson et al. (2012) and Le Hir et al. (2009) to compute this locally across the planet and then integrate to compute the global weathering rate. When the planet is not in a snowball state, implying global mean temperatures above approximately 260 K, we evaluate Equation 1 at each cell if the cell has land surface and has surface temperatures above 273.15 K, 4 times per simulated day in 6-hour intervals. The weathering rate is averaged over a full year in each cell, to compute a localized annual mean weathering rate. This is then spatially averaged over the globe to compute the annual mean global weathering rate.

We run PlaSim in the T21 configuration, with 64 latitudes and 32 longitudes. We use 10 vertical grid levels, turn on the mixed-layer ocean, and use a timestep of 45 minutes. We allow pCO2 to change year-after-year according to the annual weathering and outgassing rates. This change is usually extremely small, but it costs nothing if we are already computing the weathering rate, and yields slightly more accuracy, particularly during episodes of very fast weathering when a planet has just emerged from a snowball state. We run PlaSim for a minimum of 11 years, at which point we use NumPy to perform a linear fit to the annual mean surface temperature. We average the slopes of fits for 9-year histories ending in each of the three most recent years to give a sense of the general trend of the past decade and to try to minimize the effect of outliers. If this average slope has an absolute value greater than 0.05 K/year, then we consider PlaSim to not have reached thermal equilibrium, and allow it to continue to run for one more year. This check is considered at the end of each simulated year. This procedure has been determined empirically to achieve approximate thermal equilibrium, and other methods for determining thermal equilibrium should work as well, achieving similar results.

The zero-dimensional pCO2-evolution component of the model is computed by solving

[TABLE]

as in Menou (2015). We first compute the annual temperature change relative to the change in pCO2 over the last two years of the simulation, $\partial T/\partial\text{pCO}_{2}$. Since the surface temperature is the quantity most directly related to sea ice extent, which determines the qualitative model state, we aim to control the temperature evolution, rather than simply controlling the pCO2 evolution. This is because sea ice extent and surface temperature become increasingly sensitive to perturbations in pCO2 near the global freeze/thaw point. We therefore compute the pCO2 step size as

[TABLE]

where $\Delta T$ is the absolute value of the desired temperature change, $V$ is the volcanic CO2 outgassing rate, and $W$ is the global CO2 weathering rate. The weathering rate is only needed at this point to determine the sign of the change. In order to gain higher sampling resolution of the transitions between snowball and warm states, we decrease the temperature step size $\Delta T$ near the global freeze/thaw point. For temperatures between 240 and 268 K, we aim for a change of 0.5 K per step. For temperatures between 235–240 K and 268–273 K, we aim for 1 K per step. And for temperatures below 235 K or above 273 K, we aim for 2 K per step. We use finer sampling near the transition because a model with too-coarse sampling could artificially introduce climate cycles in models with weathering equilibrium temperatures near the global freeze/thaw point.

We can combine this pCO2 step size with the pCO2 rate of change at the end of the last model year to estimate the duration of the pCO2 adjustment:

[TABLE]

Because this estimate uses only the pCO2 rate of change at the end of the previous relaxation step, this method is likely not very accurate. However, for the purposes of identifying climate cycles and the order of magnitude of their timescales, this method is likely sufficient. A more robust and sophisticated technique could perhaps be used in its place to yield more accurate estimates of the period of a given cycle, but is beyond the scope of this paper.

II.B Establishing Limits on Climate Stability

We further attempt to determine the limits of the temperate habitable zone, beyond which limit cycles should occur. We accomplish this with a parameter sweep, running PlaSim to thermal equilibrium for 31 different insolations, ranging from 1.0 F⊕ to 0.75 F⊕, where F⊕ is the incident solar flux received by the present-day Earth, taken here to be 1367 W m*-2*. These insolations correspond to orbital radii between 1 and 1.15 AU, separated by 0.005 AU. For each insolation, we prescribe a CO2 partial pressure ranging from 0.1 $\mu$bars to 6 bars, separated by 0.1 dex, and run the model to thermal equilibrium for that given pCO2 and insolation. At the transition points between snowball and temperate climates, we reduce the sample spacing to 0.01 dex. We do this for both warm, Earth-like initial conditions with only limited glaciation, and snowball initial conditions. We use Earth-like land distributions and orography. We compute the global weathering rate averaged over two decades for each model at thermal equilibrium. These values give an estimate of the edges of the temperate habitable zone in terms of outgassing rates and insolation.

III Results

The parameter sweep, whose results are shown in Figure 1, show a range of temperatures between roughly 250 and 260 K for which no equilibriated models are found. This is because in that temperature range, sensitive albedo feedbacks with sea ice lead to either rapid glaciation into a snowball state or a rapid thaw into a temperate state. The weathering rates at the lower limit of the partially-frozen models are therefore the minimum outgassing rates required to maintain stable temperate climates. These rates are shown in Figure 2. The minimum outgassing rate is exponentially related to the insolation, such that

[TABLE]

where $\kappa$ and $\alpha$ are constants, and in our fit equal to approximately $3.6\times 10^{7}$ and 23 respectively. This has a similar functional form to the limit of the temperate habitable zone derived by Abbot (2016). Haqq-Misra et al. (2016) argued that Menou (2015) used an inappropriate estimate of Earth’s outgassing rate, citing a more appropriate value of approximately 63 bars Gyr*-1*. We test the model’s dependence on this by running the same models with 63 bars Gyr*-1* as the Earth outgassing rate. The results are shown in Figure 2 as a second set of datapoints. Since the weathering parameterization, Equation 1, uses values normalized to Earth, changing the Earth outgassing rate doesn’t affect the onset of limit cycling relative to Earth outgassing, instead affecting the timescales of the cycles. The result is that climate cycles occur at actual outgassing rates approximately 10 times higher than suggested when assuming Earth outgasses 7 bars Gyr*-1*, and at timescales approximately 10 times shorter.

Our hybrid 0D-3D models were run at $2\times 10^{-3}$ V⊕, or 14 mbar Gyr*-1*, at 1350 W m*-2*, 1325 W m*-2*, and 1300 W m*-2* insolation. We assume modern Earth obliquity and seasons, orbital eccentricity, land fraction, and continental distribution and orography. As expected, we find that these models have unstable climates characterized by limit cycling, as shown in Figure 3. The cycles have periods ranging from hundreds of millions to approximately 1 billion years. We phase-fold the warm periods to improve their resolution and increase the accuracy of our estimate of their average duration, finding that they last for approximately 1.2 Myr, 930 kyr, and 720 kyr respectively, as shown in Figure 4. The warm periods are therefore of order 0.1% of the total cycle duration at these insolations and outgassing rates. It is also of note that the mean temperature and pCO2 curves throughout the cycle have the same form as those identified in Menou (2015). Inconsistencies between different cycle periods are likely the result of imperfections in the pCO2-stepping algorithm and estimation of the time intervals between data points. With an improved model and algorithm for spanning these geological timescales, we expect greater consistency between cycles. Using the outgassing estimate in Haqq-Misra et al. (2016), we find that cycle periods decrease by roughly a factor of 10, consistent with a nearly factor of 10 increase in the speed at which atmospheric pCO2 is replenished. These are consistent with the cycle periods found by Abbot (2016).

To probe the sensitivity to continental distribution, we run a set of models at Earth insolation and 21 mbar Gyr*-1* outgassing, or 0.3% V⊕. Each has the same land fraction as modern Earth, but the land is distributed either in a rectangular supercontinent centered on the equator, a circular supercontinent centered on a pole, two circular continents centered on the poles, and three models with most land at the poles but with a small continent of varying size at the equator. Each model assumes flat orography for simplicity. Of these models, only the equatorial supercontinent demonstrates an unstable climate at this insolation and outgassing, suggesting a strong sensitivity to continental distribution, consistent with the findings in Donnadieu et al. (2006) that weathering rates are strongly-affected by changes in land distribution, with tropical continents leading to higher weathering rates.

IV Simulation Caveats

IV.A Model Uncertainties

While PlaSim represents a model of much greater complexity than the simplified model used in Menou (2015), it does have limitations. The radiation model used in PlaSim is highly simplified, largely ignoring aerosols and finer spectral variations in absorptivities, relying on broadband emmisivities, polynomial prescriptions, and constant albedos, and insufficiently treating the stratosphere, resulting in biases (Lunkeit et al., 2012). For example, different models can predict significantly different cooling rates for a given $T$-pCO2 combination, resulting in different abundances and durations of snowball states. Williams & Kasting (1997, hereafter WK97) ran a suite of 1D models for various pCO2 levels, and then fit polynomials to the radiative results of the models. They claim an rms error in outgoing longwave radiation (OLR) of less than 5 W/m2 in the 190–370 K temperature range and from 10*-5*–10 bars pressure range. Almost all of our models fall within this pCO2 range. We compute a global annual mean of outgoing longwave (thermal) radiation by recording the outgoing longwave radiation (OLR) at each grid point 4 times per simulated day. We then compare this with the OLR predicted by the WK97 polynomial fit assuming the same global annual surface temperature and pCO2, computing a single polynomial with the global values as inputs. The comparison is shown in Figure 7. PlaSim consistently predicts more cooling, by 10% or less during warm phases, but by up to 70% in the frozen phase when pCO2 levels are higher. Even in a best-case scenario of Earth insolation, using 21 mbar Gyr*-1* outgassing to trigger limit cycles, PlaSim predicts more cooling. This shortcoming of the model almost certainly affects the quantiative accuracy of the results. We however consider this loss of accuracy to be one of the tradeoffs to be made in exchange for improved model speed that allows us to simulate long timescales and large swaths of parameter space. There is a balance to be struck, and so the reader should be wary of quantitative estimates pertaining to the high-pCO2 regimes. Thus, without repeating this study with different GCMs, our results generally cannot be considered quantitative predictions for where climate cycles might occur beyond identifying qualitative trends—we expect limit cycles to dominate the outer habitable zone, with a boundary roughly described by an exponential trend, as in Equation 5 and as predicted by Abbot (2016).

PlaSim also neglects glacial dynamics—instead, snow simply accumulates on land. This means that glacial flow is ignored, and any differences in albedo or melting properties between large accumulations of snow and glaciers are ignored. And while PlaSim includes the effect on albedo of snow-covered sea ice as opposed to bare ice, the effect of dust accumulation in snow is not. This could produce qualitatively different results in how the planet melts and freezes. We also ignore any interaction between glaciation and volcanism. There likely is some important interaction, as some studies suggest that rapid changes in lithostatic pressure due to rapid melting could trigger volcanic activity (Sternai et al., 2016), and others suggest that increased lithostatic pressure on continents due to glacial loading could also trigger volcanic activity (Kyle et al., 1981). The amount of volcanic activity will also have an effect on snow and ice albedo by changing the amount of soot and other contaminants in the snow and ice (Gow & Williamson, 1971).

Our models also do not take into account changes in sea level due to melting or freezing. As weathering rate depends on land fraction, this may have a significant impact. We also do not fully-treat large-scale ocean circulation, which tends to discourage the advance of sea ice, as currents transport heat away from the equator (Poulsen et al., 2001). We also ignore variations in fundamental orbital properties such as stellar luminosity and secular variations in the orbit itself, both of which can be expected to vary on relevant timescales (Greenberg & Van Laerhoven, 2011; Gómez-Leal et al., 2016).

As a result of these approximations and inaccuracies, it is likely that there is significant uncertainty and inaccuracy in the conditions under which PlaSim predicts entry into a snowball state and deglaciation out of it. Since our aim is not to precisely quantify the conditions under which snowballs and limit cycles occur, but rather to verify that they persist as a phenomenon in higher-dimensional models and to roughly outline the parameter space in which they can be expected, we believe this is an acceptable tradeoff for PlaSim’s computational speed. However, it is worthwhile to consider how changes in snowball entry and exit conditions would affect our quantitative results. If glaciation into snowball conditions happens at colder temperatures (lower pCO2), then temperate climates should be stable at lower outgassing rates than those predicted by our results. This would generally diminish the amount of parameter space occupied by limit cycles. However, since we are concerned with a logarithmic outgassing parameter space, we do not expect a reduction in glaciation temperature of a few degrees to imply the near-absence of limit cycles from the parameter space. Conversely, if the entry into snowball happens at warmer temperatures, temperate climates would cease being stable at higher outgassing rates, increasing the prevalence of limit cycles. Changes in the deglaciation temperature should have less of an effect on limit cycle prevalence. If deglaciation happens at warmer temperatures (higher pCO2), then the post-deglaciation climate would be significantly warmer with higher weathering rates. Similarly, colder deglaciation would imply a cooler post-deglaciation climate. However, since the dominant factor in whether limit cycles occur is whether a cooling climate is able to reach a weathering equilibrium before glaciation into snowball, changing the deglaciation point leaves the parameter space unchanged. However, changes to either transition temperature could significantly affect the duration of each part of the cycle. Colder glaciation would increase the length of each phase due to the longer weathering and outgassing times needed for the larger shifts in pCO2, as would warmer deglaciation. Conversely, warmer glaciation and cooler deglaciation would reduce the weathering and outgassing times needed.

IV.B Weathering Uncertainties

The carbon-silicate cycle implementation is itself a simplification, limited by parameterizations of Earth weathering, and ignoring mantle CO2 processes. One important simplication is in precipitation and runoff. Equation 1 includes a term that depends on runoff efficiency and temperature. This is intended to account for changes in precipitation and runoff and the resulting effect on weathering (Abbot et al., 2012). However, this may miss complexities arising from glacial melt and changes in the hydrologic cycle due to large changes in ice fraction. As ice sheets advance to the equator, the amount of precipitable water present in the atmosphere is expected to dramatically decline, leading to a reduction in weathering due to precipitation and runoff. Our implementation also involves computing what was originally a global, 0D parameterization at each grid point and then taking the average, rather than taking into account such variables as the precipitation specific to that gridpoint. It is important to note that we are effectively taking the average of a polynomial, rather than computing a polynomial with an average input. A more sophisticated weathering model might yield quantitatively different results when used with our long-time evolution technique. Nor do we consider the impact of CO2 condensation on weathering feedbacks. Turbet et al. (2017) showed that snowball states at low insolation may result in large amounts of CO2 deposition as dry ice at the planet’s pole. We ignore the effect that could have on limit cycles.

We ignore seafloor weathering in our parameterization, as our interest is in the role that continental weathering can play in long-term climate stability, and the role that seafloor weathering plays as a climatic feedback is highly uncertain, as it is not well-known how much carbon is sequestered through that process, nor whether the seafloor weathering rate depends on heat flux from the mantle or the atmospheric temperature (Caldeira, 1995; Abbot et al., 2012). We can however estimate what its impact on our results might be by assuming a parameterization of the same form as that used in Abbot et al. (2012), adopting their notation of $\phi$ as the CO2 partial pressure relative to modern Earth, $pCO_{2}/pCO_{2,\oplus}$:

[TABLE]

where $W_{sf}$ is the seafloor weathering, $\nu$ is related to the seafloor spreading rate and is taken to be 12.6, and $n=\frac{1}{3}$ (Abbot et al., 2012). The total contribution of seafloor and continental weathering rates is estimated to be a linear combination,

[TABLE]

where $\beta_{0}$ is the relative modern contribution of continental weathering, and $W_{c}$ is the continental weathering rate. This is generally not going to be constant across a wide range of climates, but in the absence of good constraints, we assume it is, and set it to be 0.75 (Abbot et al., 2012). The results are shown in Figure 8, and suggest that seafloor weathering could provide a weathering floor at high insolations, such that at outgassing rates less than the seafloor weathering rate planets could become trapped in snowball conditions. This regime of parameter space however already assumes very low outgassing rates, and this particular seafloor weathering parameterization does not significantly affect the shape of the stability parameter space at lower insolations and higher outgassing rates, where it is more likely to be relevant for observed Earth-like exoplanets.

The result of these caveats is that fully modeling a planet with any predictive accuracy is difficult, requires a great deal of physics and computation, and suggests both breadth and diversity in the habitable zone exoplanet parameter space. However, by limiting the investigation to a particular slice of that parameter space, we may investigate the effects of a particular mechanism on climate. Our results suggest that the carbon-silicate feedback mechanism permits climate cycles at low insolations and moderate-to-low outgassing rates, in general agreement with the results of Tajika (2007), Menou (2015), and Haqq-Misra et al. (2016).

V Conclusion

We present a technique for simulating the evolution of climate systems on geological timescales by using a 3D GCM as a relaxation step in a 0D carbon evolution model. We use this to present the first 3D GCM simulations of the climate cycles predicted by 1D models in Menou (2015) and Haqq-Misra et al. (2016). This shows that climate cycles are not limited to simple EBMs, and strengthens the argument for their relevance to real planetary systems. We find that Earth analogs undergo climate cycles at moderate outgassing levels and low insolations comparable to early Earth. Further work includes refining the model assumptions, accounting for other physical processes such as glacial dynamics, seafloor weathering (Haqq-Misra et al., 2016), and variable outgassing rates in response to glacial loading. We note that our technique is not specific to PlaSim, and any GCM could be used in its place using the same technique. We also note that we have not mapped the entirety of the parameter space proposed by our models to exhibit climate cycles; we have only probed a few points in that space to find climate cycles. Future work should map the cycle frequencies and characteristics in this region and compare with the results presented in Menou (2015), Haqq-Misra et al. (2016), and Abbot (2016). It may be additionally useful to explore more sophisticated weathering models which capture more of the relevant physical processes, rather than merely using a 3D temperature field to compute a global average weathering rate.

VI Acknowledgments

The authors would like to thank Professor Ray Pierrehumbert for his very helpful comments, insights, and suggestions, particularly pertaining to glaciation and deglaciation conditions and prior work on coupling GCMs to processes that happen on longer timescales. The authors would also like to thank Michael Way at the NASA Goddard Institute for Space Studies for his helpful comments on the ROCKE-3D GCM and the dependence of the habitable zone on rotation period. AP is supported by a Centre for Planetary Sciences Graduate Fellowship at the University of Toronto, Scarborough, and by the Department of Astronomy & Astrophysics at the University of Toronto, St. George. KM is supported by the Natural Sciences and Engineering Research Council of Canada. Computing time was provided by the Canadian Institute for Theoretical Astrophysics at the University of Toronto, St. George. Many thanks are given to my peers for their support.

Appendix A Computational Cost

Earth-system general circulation models are often very computationally expensive. PlaSim is one of the faster GCMs available. On a compute node with two Intel®  Xeon®  E5310 quad-core CPUs, clocked to 1.60 GHz, running PlaSim in parallel with 8 threads takes 162 seconds of wall-time on average per simulated year, where wall-time is the real time that a user must wait for a job to finish. On a compute node with 2 octa-core Intel®  Xeon®  E5-2640 CPUs, clocked at 2.00 GHz, running PlaSim in parallel with 16 threads takes 62 seconds of wall-time on average per simulated year. Overall, PlaSim takes approximately 1300 CPU-seconds per simulated year on the E5310 cores and 1000 CPU-seconds per simulated year on the E5-2640 cores. As an example of a more sophisticated GCM, the ROCKE-3D Earth GCM developed at the NASA Goddard Institute for Space Studies can compute a simulated year in as little as 11 CPU-hours (38,000 CPU-seconds) on their Xeon®  E5-2697 nodes using 44 cores (Way et al., 2017). Climate cycles, as proposed by Menou (2015), occur on the order of hundreds of millions of years to billions of years. Taking as an optimal computing case, if we assume a climate cycle could be found which lasts only 100 Myr, then on the computing hardware available to us, computing the entire cycle directly year-by-year with PlaSim on 16 cores would take $6.2\times 10^{9}$ seconds of wall-time, or approximately 196 years. Since submitting a computing job with a run-time of two centuries is not practical or desirable, directly simulating the carbon-silicate cycle with a GCM is impossible without considerable computing resources. Our technique, however, allows us to avoid directly computing each year of climate evolution, at the cost of computational accuracy. Instead, we compute only a few decades intermittently throughout the long-term evolution of the planet. Computing two periods of the cycles observed in this paper involved directly simulating approximately 1800 years, which is 6 orders of magnitude less than directly computing the entire cycle, and allows full computation in only 31 hours of wall-time on our hardware.