The “Snowball Earth” hypothesis, proposed to explain the Neoproterozoic glacial episodes in the period 750–580 million years ago, suggested that the earth was globally covered by ice/snow during these events. This study addresses the problem of the forcings required for the earth to enter such a state of complete glaciation using the Community Climate System Model, version 3 (CCSM3). All of the simulations performed to address this issue employ the geography and topography of the present-day earth and are employed to explore the combination of factors consisting of total solar luminosity, CO2 concentration, and sea ice/snow albedo parameterization that would be required for such an event to occur. The analyses demonstrate that the critical conditions beyond which runaway ice–albedo feedback will lead to global freezing include 1) a 10%–10.5% reduction in solar radiation with preindustrial greenhouse gas concentrations; 2) a 6% reduction in solar radiation with 17.5 ppmv CO2; or 3) 6% less solar radiation and 286 ppmv CO2 if sea ice albedo is equal to or greater than 0.60 with a snow albedo of 0.78, or if sea ice albedo is 0.58 with a snow albedo equal to or greater than 0.80. These bifurcation points are very sensitive to the sea ice and snow albedo parameterizations. Moreover, “soft Snowball” solutions are found in which tropical open water oceans stably coexist with year-round snow-covered low-latitude continents, implying that tropical continental ice sheets would actually be present. The authors conclude that a “soft Snowball” is entirely plausible, in which the global sea ice fraction may reach as high as 76% and sea ice margins may extend to 10°S(N) latitudes.
During the Neoproterozoic era, it is generally argued that two or more global-scale glaciations occurred, one that began at ~716 Ma (Ma ≡ 106 yr ago: the Sturtian glaciation) and a second that began at ~635 Ma (the Marinoan glaciation) (e.g., Macdonald et al. 2010), during which continental ice sheets reached low latitudes and the ocean was completely covered by thick sea ice. This will be referred to as the “Snowball” or “hard Snowball” hypothesis (Harland and Bidgood 1959; Kirschvink 1992; Hoffman et al. 1998; Hoffman and Schrag 2002). An alternative to this hard Snowball hypothesis was first proposed by Hyde et al. (2000) and further developed in Peltier et al. (2004, 2007) and is referred to in the literature as the “soft Snowball” or “Slushball” hypothesis. In the latter model, although the continents are covered by thick continental ice sheets, stable states are suggested to exist in which open water in the tropics may coexist with low-latitude continental glaciers. A further alternative hypothesis is that suggested by Pollard and Kasting (2005), which involves a so-called “thin ice” solution, in which the global ocean is covered by sea ice but remains thin enough to permit sunlight to reach the underlying ocean so that photosynthetic biota are able to survive the otherwise extreme “ice house” conditions. Analyses by Warren and Brandt (2006) suggest that the “thin ice” solution would be stable only in an extremely restricted region of parameter space. A further alternative solution to the “Neoproterozoic enigma” that has been suggested is the “high obliquity” hypothesis, in which the tilt of the spin axis with respect to the plane of the ecliptic is suggested to have been in excess of 54° (Williams et al. 1998; Williams 2008). In this state, the present-day “tropical” surface would be colder than the poles so that the existence of low-latitude continental ice sheets would not imply global glaciation at all. This suggestion is, however, in conflict with the results of geodynamic analyses that strongly suggest the earth’s obliquity has remained close to its present-day value (e.g., Levrard and Laskar 2003). We will therefore focus in this present paper on the assessment of the relative merits of the hard Snowball and soft Snowball hypotheses.
The issue of the forcings that would be required to enable the earth to enter a state of global glaciation is clearly one of the most critical issues in the Snowball Earth debate. The models that have been employed to investigate this question have ranged from simple quasi-analytic energy balance models (EBMs) to more realistic models of the coupled atmosphere, ocean, sea ice, and continental ice sheet system. Based on the application of simple one-dimensional EBMs, including the influence of ice–albedo feedback, Budyko (1969) and Sellers (1969) showed that a modest decrease of solar radiation might be sufficient to induce ice age conditions; a 10% reduction in the solar constant was suggested to possibly lead to catastrophic descent into a completely ice-covered earth state, as summarized in the review by North and Stevens (2006). However, the results of energy balance models are very sensitive to the specified thermal diffusion coefficient required to represent meridional heat transport (North 1975; Ikeda and Tajika 1999).
It is furthermore the case that such EBMs are unable to estimate the influence of clouds, of lapse rate feedback, of atmospheric and oceanic heat transports, of sea ice dynamics, and of ice sheet dynamics, which suggests that much more realistic and complex models must be employed to study this problem. Using a 2D energy balance model coupled with an ice sheet model (EBM–ice sheet model), Hyde et al. (2000), Crowley et al. (2001), Peltier et al. (2004), and Liu and Peltier (2010) found that thick low-latitude continental ice sheets may in fact coexist with a substantial area of open water in the tropical region and the open ocean fraction may be 50% or greater. By inserting both ice sheet extent and topography from the EBM–ice sheet coupled model into an atmospheric general circulation model coupled with a slab ocean, Hyde et al. (2000) and Baum and Crowley (2001) found that, although the region of open water on the equator is reduced, it remains stable in the presence of the continental ice sheets when the CO2 concentration is ~2.5 times the present concentration.
Moreover, by employing an atmospheric general circulation model coupled to a slab ocean, various authors have reported that the earth would enter into a completely frozen state when the models are subjected to extreme forcings; however, the required forcings were considered unrealistic for the Neoproterozoic era. For example, a 13%–14% reduction in solar radiation has been found to induce a hard Snowball state (Stone and Yao 2004), but the solar constant during the Neoproterozoic era was only 6%–7% less than present (Gough 1981; Crowley and Baum 1993). In a similarly extreme example, in which the CO2 concentration was effectively reduced to 1 ppmv in the general circulation model Portable University Model of the Atmosphere (PUMA), the tropical ocean was found to remain free of sea ice (Romanova et al. 2006). These results have begun to be placed in content through the application of state-of-the-art coupled atmosphere–ocean general circulation models (AOGCMs). The most recent analysis is that which has been performed using the ECHAM5/Max-Planck-Institute for Meteorology Hamburg Primitive Equation Ocean Model (MPI-OM) in which a 6%–9% reduction in the solar constant was found to initiate a “modern” hard Snowball Earth with the present geography and topography (Voigt and Marotzke 2009) and a 4%–4.5% reduction in solar constant was found to be sufficient to initiate a Marinoan hard Snowball Earth (Voigt et al. 2011). Since these conditions do suggest the plausibility of hard Snowball states during the epoch of time, their results certainly warrant close scrutiny and detailed checks using other equally capable modern models of the coupled system. This is especially so that this result contrasts with that obtained using the AOGCM Fast Ocean Atmosphere Model (FOAM), in which a 9% reduction in solar radiation and 140 ppmv CO2 was required to trigger entry into the hard Snowball state (Poulsen and Jacob 2004).
A further peculiarity in the literature on this subject is the finding that there are no stable states with global sea ice cover exceeding 50%, corresponding to sea ice margins equatorward of ~30° latitude (Poulsen and Jacob 2004; Voigt and Marotzke 2009; Voigt et al. 2011), consistent with the results of energy balance models (Budyko 1969; Sellers 1969; North 1975), earth climate system models of intermediate complexity (EMICs) (Bendtsen 2002; Donnadieu et al. 2004; Lewis et al. 2007), and some GCMs (Jenkins and Smith 1999; Romanova et al. 2006). However, in other GCMs, various authors have found that a narrow strip of sea-ice-free tropical ocean may remain stable even as the global sea ice cover reaches as high as 80% (Chandler and Sohl 2000; Micheels and Montenari 2008; Pierrehumbert et al. 2011; Abbot et al. 2011). Based on the significant discrepancies between the results obtained using these different models described above, this issue also warrants further scrutiny and is another of the several purposes of the present work.
Using a fully coupled modern ocean–atmosphere model, the National Center of Atmospheric Research Community Climate System Model, version 3 (CCSM3), we will reestimate the critical conditions in which a hard Snowball Earth would be forced to develop. Our analyses are formulated so as to provide a direct comparison with the work of Voigt and Marotzke (2009) based on the ECHAM5/MPI-OM. We design a series of experiments to systematically examine the response of the climate system to different forcings, including solar luminosity, carbon dioxide concentration, initial conditions, and the insolation variations associated with the changing geometry of Earth’s orbit around the sun. We will focus on the dependence of the results on the details of the sea ice/snow parameterization, including sea ice and snow albedos, surface melting of snow, and of the underlying sea ice, the vertical resolution of the sea ice (number of sea ice layers), the sea ice thickness distribution, and sea ice dynamics.
The paper is organized as follows. Section 2 describes the climate model, and section 3 describes the design of the series of simulations we will perform. Section 4 reports our results. We will focus on the critical point at which a hard Snowball Earth occurs and on the range of conditions over which the existence of a soft Snowball solution is expected. In these analyses, when the sea ice margins are outside of the tropics and the sea ice coverage is greater than 0% and less than 50%, this state will be referred to as a partial-ice-covered solution. When the sea ice margins enter the tropics, the sea ice coverage is greater than 50% but less than 100%, and, while a significant part or all of the tropical continents are perennially snow covered, such a state will be referred to as a soft Snowball solution. When the global ocean is covered by sea ice/snow, it will be referred to as a hard Snowball solution. Section 5 presents a detailed discussion and interpretation of the results of section 4. Section 6 offers a summary of our conclusions. In Yang et al. (2012, hereafter Part II), we proceed to discuss the feedbacks that are engendered in the model when either the solar constant or the CO2 concentration is reduced, feedbacks that act to further decrease (or increase) the surface temperature and thus promote (or prevent) the equatorward expansion of sea ice, including the feedbacks associated with clouds, water vapor, sea ice dynamics, the Hadley cells, and oceanic wind-driven and thermohaline circulations.
2. Model description
The Community Climate System Model version 3 is a state-of-the-art coupled model for simulating past, present, and future climates. CCSM3 has four major components (atmosphere, ocean, land, and cryosphere) linked via a flux coupler, where no corrections are applied to the fluxes between components (Collins et al. 2006b). The atmospheric part, the Community Atmosphere Model version 3 (CAM3) is a three-dimensional primitive equation model (Collins et al. 2004). The land component, the Community Land Model version 3 (CLM3) (Oleson et al. 2004), treats five surface types for each grid cell, namely, glacier, lake, wetland, urban, and dynamic vegetation. The ocean circulation model is based on the Parallel Ocean Program (POP), a three-dimensional primitive-equation-based model using the hydrostatic and Boussinesq approximations, originally developed at the Los Alamos National Laboratory (Smith et al. 2010). The sea ice component is the Community Sea Ice Model version 5 (CSIM5) (Briegleb et al. 2004), which includes an energy-conserving thermodynamics model (Bitz and Lipscomb 1999) and an elastic–viscous–plastic representation for sea ice dynamics (Hunke and Dukowicz 1997, 2002).
The model is designed to enable a wide range of spatial resolutions. Here, we use the most economical low-resolution, T31 grid (an approximately 3.75° × 3.75° resolution transform grid in which the atmosphere is described on 26 vertical levels); the ocean is described with a spatial discretization, labeled “gx3,” which has 3.6° resolution in the zonal direction but variable resolution in the meridional direction with a spacing of approximately 0.6° near the equator, and 25 vertical levels. The land model employs the same horizontal grid as the atmospheric component, and each grid cell of the land model incorporates a subgrid hierarchy, including multiple land unit types, snow/soil columns, and plant functional types (PFTs). The horizontal grid of the sea ice model is the same as that of the ocean and describes the sea ice distribution in five categories with four vertical layers for each thickness category on the basis of which the model computes the vertical heat conduction through the ice. The time steps for atmosphere, land, sea ice, and ocean are 30, 30, and 60 min, and 2 h, respectively, and the coupling interval is 1 h between the atmosphere, land, and sea ice but 1 day between the ocean and the other three components. When applied to simulate the present-day climate, many aspects of the prediction of the low-resolution model compare quite favorably with observations. For example, the strength of the Atlantic overturning circulation is approximately 16 Sv (Sv ≡ 106 m3 s−1) (Yeager et al. 2006), which is well within the error bars of observations (Cunningham et al. 2007).
The albedo of glacier over land is set to be 0.80 in the visible band and 0.55 in the near-infrared band, while frozen wetland and lake albedo is assumed to be 0.60 in the visible and 0.40 in the near infrared. The albedo of soil varies with soil color and ranges from 0.05 to 0.24 in the visible and from 0.10 to 0.48 in the near infrared. The albedo of snow over land depends on the snow depth, snow aging, the cosine of the solar zenith angle, and soot and dust accumulation. The typical albedo of new snow is 0.95 in the visible and 0.65 in the near infrared, which is respectively 0.04 and 0.02 higher than those of snow over sea ice because snow over land contains relatively low moisture. The upper limit of snow depth over land is fixed at 1.0 m in liquid water equivalent. For a more detailed description of the land surface albedo parameterization, see Dickinson et al. (1993) and Oleson et al. (2004).
For Snowball Earth simulations, the treatment of sea ice and snow is extremely important (Lewis et al. 2006; Pierrehumbert et al. 2011), so a more detailed discussion of the sea ice module of CCSM3 is warranted. The growth of sea ice is determined by a combination of thermodynamic and dynamic contributions. The energy balance at the top and bottom surfaces determines the thermodynamic contribution, while the dynamic contribution depends on the Coriolis force, surface wind stress, the stress applied to the underside of the sea ice by ocean currents, the internal stress within the ice, and the stress due to sea surface tilt of the geoid (Briegleb et al. 2004). Furthermore, the thermodynamic contribution includes six subcomponents: snow–ice conversion (snoice), surface ice melt (meltt), basal ice melt (meltb), lateral ice melt (meltl), congelation ice growth (growb), and frazil ice growth (frazil). Any ice that forms at depth is assumed to float upward to the surface, where it may either partially or completely melt in its upper layers where temperatures may be above the freezing point. The freezing point of seawater is set to be −1.8°C, which is assumed to be independent of oceanic salinity. The sea ice salinity is assumed to be 4 psu. When ice forms or melts, potential temperature and salinity are both adjusted so as to ensure energy and salt conservation (Smith et al. 2010).
In several respects, the sea ice thermodynamic and dynamic modules of CCSM3 are more accurate and realistic than those in ECHAM5/MPI-OM. For the sea ice thermodynamics, the zero-layer Semtner model used in ECHAM5/MPI-OM (Semtner 1976) overestimates the amplitude of the annual cycle of the sea ice compared with the energy-conserving model in CCSM3 (Bitz and Lipscomb 1999). For sea ice dynamics, the elastic–viscous–plastic (EVP) model in CCSM3 obtains more accurate ice transient behavior than the viscous–plastic (VP) model (Hibler 1979) used in ECHAM5/MPI-OM (Hunke and Dukowicz 1997).
Despite the accurate thermodynamic and dynamic treatments of sea ice evolution, the albedo parameterization for sea ice and snow is relatively crude. In CCSM3, the albedos depend on spectral band (visible, wavelength <0.7 μm, and near infrared, wavelength >0.7 μm), snow thickness hs, ice thickness hi, and surface temperature TS (Briegleb et al. 2004; Hunke and Lipscomb 2010). When snow depth is above 0.02-m liquid water equivalent and surface temperature is below −1°C, snow albedo is 0.78 [0.91 in the visible and 0.63 in the near infrared, assuming the partition between the visible band and the near-infrared band is 53%/47% (Briegleb et al. 2004)], which will be referred as “cold snow albedo” or “snow albedo” Hereafter. When sea ice is thicker than 0.5 m and surface temperature below −1°C, the sea ice albedo is set to be 0.50 (0.68 in the visible and 0.30 in the near infrared), which will be referred as the “cold sea ice albedo” or “sea ice albedo” hereafter.
As surface temperature increases above the −1°C (“temperature criterion”), the sea ice albedo and the snow albedo decrease such that by TS = 0°C, a “warm sea ice albedo” and a “warm snow albedo” are assumed to be 0.425 and 0.656, respectively. This temperature dependence of sea ice and snow albedo is employed as a means of accounting for the potentially important influence of snow wetness and the development of melt ponds on the surface of the sea ice. This is illustrated in Fig. 1, not only for CCSM3 but also for other models: ECHAM5/MPI-OM, FOAM, and the Geophysical Fluid Dynamics Laboratory Climate Model version 2.1 (GFDL CM2.1). It is notable that different models use different sea ice and snow albedos, and the temperature dependence of the parameterizations also differs significantly. The temperature criterion is −10°C in FOAM and GFDL CM2.1, but it is −1°C in CCSM3 and ECHAM5/MPI-OM. Both of these assumed critical temperatures for the onset of melt ponds are somewhat arbitrary, the accuracy of which is not yet well understood (Briegleb and Light 2007). Unfortunately, the temperature criterion is a critical factor. As we will demonstrate in section 4e, near the ice margins from May to October, the surface temperatures are always between −4° and −10°C; therefore, a −1°C temperature criterion will have no influence on the surface albedo but a −10°C temperature criterion will significantly decrease the surface albedo during these months. We will test the influence of the temperature criterion on the results, which indirectly relates to the influence of melt ponds and leads.
As ice thickness decreases from 0.5 m (“sea ice thickness scale”) to 0 m, the ice albedo decreases smoothly to the ocean albedo 0.06, and, as snow depth decreases from 0.02-m liquid water equivalent (“snow depth scale”) to 0 m, the snow albedo decreases smoothly to zero. The sea ice thickness scale of 0.5 m is comparable to observations (Briegleb and Light 2007), but the snow depth scale of 0.02 m is much larger than observations. As reviewed by Warren et al. (2002), 5–10 mm of fresh snow can increase the sea ice albedo from 0.49 to 0.81. In the extremely cold conditions of the Snowball Earth, the snow, when it exists, is always thicker than 0.1 m (see section 4b). We have designed an experiment in which the snow depth scale is set to be half of the default value and find that it has essentially no effect on the results. Therefore, it may well be that this crude assumption of a fixed snow depth scale (0.02-m liquid water equivalent) will not significantly affect the results.
The albedos in CCSM3 are within the uncertainty of observations in the Arctic and Antarctic (Perovich et al. 2002; Brandt et al. 2005) but, when new cold dry snow covers the sea ice, the albedo is ~0.07 lower than observations (Collins et al. 2006b). This diminished albedo is tuned so as to compensate the ~50 W m−2 less incoming shortwave radiation due to other biases in the CCSM3 model, such as the poor Arctic low cloud cover in summer (Gent et al. 2011), and, with this simple adjustment, the model correctly simulates the timing of the onset of sea ice/snow melting (Collins et al. 2006b). Moreover, it is clear that the Intergovernmental Panel on Climate Change Fourth Assessment Report (IPCC AR4) climate models systematically overestimate the sea ice/snow albedo by 0.05 in summer over the northern latitudes of the Western Hemisphere (Wang et al. 2006).
As shown in Fig. 1, the sea ice albedo of CCSM3 is lower than that in the other three models. For example, in the ECHAM5/MPI-OM experiments of Voigt and Marotzke (2009) and Voigt et al. (2011), the cold sea ice albedo is assumed to be 0.75—which is 0.25 higher than that in CCSM3; much higher than observations of sea ice, 0.47–0.52, and sea glacier, 0.55–0.66 (Perovich 1996; Warren et al. 2002; Warren and Brandt 2006); and equal to the observed maximum value for bare subeutectic sea ice (TS = −37°C) containing hydrohalite crystals in brine inclusions (Warren et al. 2002). The sea ice albedo of 0.75 may be considered unrealistic, especially near sea ice fronts where the surface temperatures are always above −37°C (see Fig. 12b). We have tested different albedo settings within the CCSM3 structure and find that it is a very important factor for Snowball Earth initiation, similarly important to the issue of tropical sea ice thickness for photosynthesis under the hard Snowball Earth conditions (Pollard and Kasting 2005; Warren and Brandt 2006) and to the issue of the conditions required for deglaciation of the hard Snowball Earth (Pierrehumbert 2004, 2005; Lewis et al. 2006; Abbot and Pierrehumbert 2010; Hu et al. 2011).
Furthermore, CCSM3 is one of only two models that can successfully forecast the decline of the Arctic sea ice extent in September from 1979 to 2006 (Meehl et al. 2006; Stroeve et al. 2007). The other similarly successful model is the Hadley Centre Global Environmental Model version 1 (HadGEM1) (McLaren et al. 2006). Notably, both of these models incorporate relatively realistic sea ice models, for example, by including a subgrid-scale sea ice thickness distribution (ITD), which is important for sea ice albedo feedback (Stroeve et al. 2007; Holland et al. 2006). We will test the influence of the ITD on the results for the Snowball Earth development in section 4f.
3. Design of the snowball transition experiments
In all simulations to be described, the land–sea distribution, land surface type, surface topography, and ocean bathymetry are maintained at their present-day conditions as in the experiments with the ECHAM5/MPI-OM (Voigt and Marotzke 2009). A zonally symmetric latitude–height distribution of ozone is specified, corresponding to the monthly mean climatology used in the second phase of the Atmospheric Model Intercomparison Project. A complete list of the design criteria applied to each of the simulations in this study is provided in Table 1. A detailed description of the applied forcings and the sea ice/snow albedo sets for each of the simulations to be performed are presented.
a. Solar luminosity
In the time interval 750–580 Ma, solar luminosity was 6%–7% less than present (Gough 1981; Crowley and Baum 1993). The present total solar irradiance is 1367 W m−2, termed “TSI0” in what follows. A sequence of simulations has been carried out with the solar constant set to the sequence of values: 0.99TSI0, 0.98TSI0, 0.97TSI0, 0.96TSI0, 0.95TSI0, 0.94TSI0, 0.92TSI0, 0.91TSI0, 0.90TSI0, 0.895TSI0, 0.89TSI0, 0.88TSI0, 0.85TSI0 and 0.83TSI0; the number in front of TSI0 represents the fraction of the present-day solar irradiance. For each 1% decrease in solar luminosity, the reduction in net solar radiation is approximately 2.4 W m−2 at the top of the model [1367/4 × 0.01 × (1 − α); assume a value of 0.3 for planetary albedo α; the number 4 is the ratio of the area of a sphere to a disk]. All of these experiments are initialized from the present-day climate by abruptly decreasing the solar constant, and most of the simulations are integrated for as long as 3000 model years. In this group of experiments, the greenhouse gases are fixed at preindustrial levels.
b. Greenhouse gases
According to the Snowball Earth hypothesis, onset of these glacial episodes may have been triggered by extremely weak greenhouse forcing due to a low atmospheric CO2 concentration (Kirschvink 1992), but what this level may have been is unknown. The control experiment (1.0TSI0) for this sequence of analyses is integrated under the “present day” greenhouse gas conditions (assumed to be 355 ppmv CO2). For the other experiments, the sequence of CO2 levels employed in the set of experiments is 286, 140, 70, 50, 35, 30, 25, 20, or 17.5 ppmv. Recent estimates show that the radiative forcing (longwave plus shortwave) due to doubled CO2 is ~2.66 W m−2 at the top of the model (Collins et al. 2006a). Other atmospheric greenhouse gases are fixed at preindustrial levels, CH4 = 805.6 ppbv, N2O = 276.7 ppbv, and no chlorofluorocarbons (CFCs). However, during the Proterozoic Eon (from 2500 Ma to 542 Ma), the CH4 concentration may have been higher than those specified here (Poavlov et al. 2003). In this group of experiments, the solar radiation is fixed at a value that is 6% less than the present level and, therefore, closer to the Marinoan level (~635 Ma) and slightly higher than the Sturtian level (~716 Ma). All of these experiments are initialized from the present-day climate.
c. Sea ice/snow albedo
In CCSM3, by default, the sea ice albedo is assumed to be equal to 0.50 and the snow albedo to 0.78; hereafter, this albedo set is denoted as (0.50, 0.78). We have performed a series of experiments with different sea ice and snow albedos as follows: (0.55, 0.80), (0.58, 0.78), (0.58, 0.80), (0.60, 0.80), (0.65, 0.80), and (0.75, 0.80). In all of these simulations, as the surface temperature increases from −1° to 0°C, the sea ice albedo and the snow albedo decrease by 0.075 and 0.124, respectively. In this group of experiments, solar luminosity is fixed at 94% of the present level, and the carbon dioxide amount is set at 286, 572 (2 times), and 1144 ppmv (4 times). These experiments are initialized from the present-day climate, from a quasi-equilibrium state of the experiment with 94% solar radiation and 35 ppmv CO2 (approximately 65% sea ice coverage), or from a quasi-equilibrium state of the experiment with 94% solar radiation and 70 ppmv CO2 (approximately 40% sea ice coverage).
d. Melting snow and melting sea ice
Observations show that melt ponds can significantly decrease the summer surface albedo (Perovich et al. 2002; Eicken et al. 2004), and model simulations have provided further support for the importance of this phenomenon (Ebert and Curry 1993; Curry et al. 1995; Lüthje et al. 2006; Pedersen et al. 2009). To perform a simple test of the importance of melting snow and melting sea ice in CCSM3, we have designed three experiments—with fixed sea ice albedo of 0.65 (no ice melting, Fig. 1a), fixed snow albedo of 0.83 (no snow melting, Fig. 1b), and a linear temperature-dependent sea ice albedo and snow albedo (ice/snow melting, Fig. 1c) as for sea ice
and for snow
where αice is sea ice albedo, αsnow is snow albedo, and TS is surface air temperature at 2 m. In this group of experiments, solar luminosity is fixed at 94% of the present level and carbon dioxide is fixed at 286 ppmv. These experiments are initialized from a quasi-equilibrium state of the experiment, with 94% solar radiation and 35 ppmv CO2 (approximately 65% sea ice coverage).
e. Sea ice vertical resolution and sea ice thickness distribution
By default, the vertical distribution of sea ice in CCSM3 is described for four layers, three layers in FOAM and a single layer in ECHAM5/MPI-OM; all of these models poorly resolve the penetration depth of the diurnal cycle. A vertical grid spacing of less than 0.3 m is needed to accurately simulate the diurnal cycle of sea ice surface temperature (Abbot et al. 2010). Here, we increase the number of the sea ice layers to 12 or 60 to test the influence of this increased resolution on the result. We have also tested the influence of the sea ice thickness distribution. By default, the sea ice is represented in five categories and each category has a specified upper limit on thickness. To find the optimal number of categories, we decrease the number to 3 and increase it to 10. In this group of experiments, solar luminosity is fixed at 94% of the present level and carbon dioxide amount at 35 ppmv. All of these experiments are initialized from the present-day climate, and the default sea ice albedo and snow albedo are employed, 0.50 and 0.78, respectively.
f. Sea ice dynamics
To test the role of sea ice dynamics, we have performed two types of experiments, one with both dynamic and thermodynamic sea ice, and the other with only thermodynamic sea ice. The solar radiation is fixed at 89% or 90% of the present level, and the carbon dioxide concentration is fixed at 286 ppmv. These experiments are initialized from the present-day climate or from a quasi-equilibrium state of the experiment with 94% solar radiation and 25 ppmv CO2 (approximately 68% sea ice coverage). In this group of experiments, the default sea ice/snow albedos are employed.
g. Different initial conditions
Because the climate system is strongly nonlinear, the final equilibrium state depends on the initial conditions. We have tested the strength of this dependence by employing different initial sea ice fractions, accompanied by different atmospheric and oceanic temperatures. Five different initial sea ice fractions have been tested: 21%, 65%, 69%, 76%, and 81%. These initial conditions are taken from different transient states of the run with 89% solar radiation and 286 ppmv CO2, which finally enters a globally ice-covered state. Solar luminosity is fixed at 90% or 94% of the present level, and carbon dioxide concentration is fixed at 286 ppmv. In this group of experiments, the default sea ice/snow albedos are employed.
h. Orbital geometry
In all of the previously described experiments, the eccentricity, obliquity, and precession of Earth’s orbit around the sun were fixed to the modern values 0.0167, 23.44°, and 102.72°, respectively. During the period of deep glaciation, however, Earth’s orbit could have been significantly different, which would have had some effect on the extent of sea ice, especially for the soft Snowball if such a state did in fact exist.
Here, we test the impact of the above parameters for three significantly different orbital regimes, the mid-Holocene warm period at 6000 yr before present (BP) (0.0187, 24.11°, and 0.87°), the Last Glacial Maximum (LGM) at ~21 000 yr BP (0.0190, 22.95°, and 114.43°), and that corresponding to the end of Eemian interglacial (PEGI) at 116 000 yr BP (0.0414, 22.49°, and 94.17°) (Peltier 2007). For the purpose of these experiments, we have not changed the solar constant and have kept it fixed to the value of 1285 W m−2, that is, 6% less than the present-day level. All of these experiments are initialized with a soft Snowball solution in which the sea ice margins are at ~20°S (N) latitude and the carbon dioxide concentration is fixed at 35 ppmv. In this group of experiments, the default sea ice/snow albedos are employed.
Although the rate of the earth’s rotation decreases over time due to lunar tidal stresses and the length of the solar day was ~21.9 h at 650 Ma (Williams 1990), previous results have shown that the effect of rotation rate is minor in comparison to CO2 concentration and surface albedo (Jenkins and Smith 1999; Pierrehumbert et al. 2011). Jenkins and Smith found that, when increasing the rotation rate from one corresponding to an 24-h day to one corresponding to an 18-h day (at 900 Ma), the annual and global mean surface temperature decreases by only 0.6 K. We have therefore elected not to directly investigate this additional influence.
a. Reducing solar luminosity
Sea ice evolution in the reduced solar radiation experiments is shown in Fig. 2. In the control experiment with the present solar constant (1367 W m−2) and greenhouse gases, the global-mean TS of 287.1 K and sea ice cover of ~6% are close to modern values (Peixoto and Oort 1992). After a 1% reduction in solar radiation and 69 ppmv decrease in CO2 (from 355 to 286 ppmv CO2), the global-mean TS decreases by 3.2 K and the sea ice cover increases from approximately 6% to 11% (Fig. 2). In further experiments, we simply reduce the solar radiation in a series of successive steps while keeping the CO2 concentration fixed at 286 ppmv. Our analyses demonstrate that runaway albedo feedback occurs after a 10.5% reduction in solar luminosity. From 0.90TSI0 to 0.895TSI0, the TS decreases by 23 K and the global-average sea ice fraction jumps from 49% to 100%.
The bifurcation point (above which runaway albedo feedback occurs), namely, a 10%–10.5% reduction in solar radiation in CCSM3, is close to that previously predicted using energy balance models (10%, North and Stevens 2006) and FOAM (9% for an ideal, tropical supercontinent, Poulsen and Jacob 2004) but much higher than those in ECHAM/MPI-OM [6%–9% for the modern geography (Voigt and Marotzke 2009); 4%–4.5% for the Marinoan geography (Voigt et al. (2011)].
b. Reducing CO2 concentration
By fixing solar radiation to be 94% of the present level, the sea ice cover response to carbon dioxide reduction is shown in Fig. 3. When reducing the CO2 concentration from 286 to 140 ppmv or from 140 to 70 ppmv, the global-mean sea ice fraction (ICEFRAC) increases by ~5%, whereas the reduction from 70 to 50 ppmv leads to an increase in ICEFRAC from approximatey 40% to 60%. This enhanced sea ice expansion is due to the influence of the Hadley cells. As the sea ice enters the domain of the Hadley cells (corresponding to a sea ice coverage of ~40% under modern geography), the sea ice flows are greatly enhanced by strong stresses associated with the trade winds (the near-surface branch of the Hadley cells, see Part II); the trade winds transport relatively cold air from the subtropics to the sea ice fronts (Bendtsen 2002) that also enhances the equatorward expansion of the sea ice. In other words, it is likely that there are no stable states between ~40% and 60% sea ice coverage during the initiation of the Snowball Earth; this phenomenon is further confirmed in the simulations with CCSM4 (Yang and Peltier 2012).
Upon further decrease of the atmospheric CO2 concentration from 35 to 17.5 ppmv, the ocean becomes totally ice covered at model year 450. This bifurcation point at 17.5 ppmv CO2 and 94% of the modern solar radiation in CCSM3 is lower than those reported for the Goddard Institute for Space Studies GCM (<40 ppmv CO2, Chandler and Sohl 2000) and in the EMIC Climate and Biosphere Group model, version 2 (CLIMBER-2) (89 or 140 ppmv, Donnadieu et al. 2004), but higher than that in the PUMA GCM (less than 1 ppmv, Romanova et al. 2006). Importantly, when the CO2 concentration is between 50 and 20 ppmv, the maximum stable sea ice cover (above which runaway albedo feedback occurs) may be as high as 76% and the sea ice margins may extend close to 10°S and 10°N latitude in CCSM3 (Fig. 4). In GCMs that do not incorporate active ocean dynamics but are run with specified ocean heat transports, various authors have also found similar states (Baum and Crowley 2001; Chandler and Sohl 2000; Micheels and Montenari 2008).
As clearly established by the sea ice margin time series shown in Fig. 4, when the sea ice approaches the equator, the amplitude of the seasonal cycle of the sea ice edge becomes increasingly large. This is the key process in determining the seasonal advance and retreat of the sea ice fronts and thus the strength of ice albedo feedback, as will be discussed further in section 4e.
The spatial patterns of snow, sea ice, and open water are illustrated in Fig. 5. Two points are notable in this sequence of graphics. 1) As the sea ice front extends into the region of 30°S–30°N, tropical continents become permanently snow covered (an indication that continental ice sheets would develop in these regions), while substantial open water can still persist in the tropical oceans (Fig. 5e). This state is consistent with the spatial distribution of Neoproterozoic glacial deposits (Hoffman et al. 1998; Evans 2000), although the positioning of the continents is not in accord with the reconstructions of either the Sturtian or the Marinoan ice ages (Li et al. 2008). 2) Over the surfaces of the sea ice and the continents at similar latitudes, there always exists thick snow (~1.0 m) (Figs. 5a–d). However, when the sea ice fronts extend to within ~18°S–18°N, the surface of the sea ice remains snow free, especially in the Southern Hemisphere (Fig. 5e), which serves to significantly limit the strength of the albedo feedback.
As shown in Fig. 6, the surface temperatures over the open water regions of the oceans lie between 0° and 7°C, close to the results obtained using the GCM GENESIS version 2 with specified ice sheets (Hyde et al. 2000; Baum and Crowley 2001). In the same circumstances, the surface air temperatures over the tropical snow-covered continents lie between 0° and −10°C, much higher than those predicted using the GENESIS-2 GCM. In the GENESIS-2 model, the surface temperatures over the continents were predicted to lie between −20° and −40°C, mostly due to the assumed high-altitude topography of the ice sheets (Hyde et al. 2000; Baum and Crowley 2001). Since the present analyses are based upon the assumption of modern geography and topography, this influence is not included.
The contributions of various processes to determine sea ice thickness in the final steady state of the soft Snowball Earth are depicted in Fig. 7. The thickness of sea ice has a strong meridional gradient, ranging from ~2 m near the ice edge to tens of meters at high latitudes. The sea ice distribution is controlled by two primary physical influences—thermodynamics and dynamics. The contributions of these influences are effectively equal but of opposite sign at all latitudes, and the net sea ice growth rate is close to zero at each latitude (Fig. 7a). The dynamical contribution is negative at mid- and high latitudes and positive in the tropics, which simply indicates that the sea ice is transported equatorward. Furthermore, nearly all of the sea ice transported into the tropics is melted by the incoming solar radiation over the surface, the cloud longwave forcing, and the poleward heat flux associated with the wind-driven oceanic circulations (for additional detail, see Part II).
The thermodynamic contribution includes six subcomponents: basal ice melt, surface ice melt, lateral ice melt, congelation ice growth, snow-to-ice conversion, and frazil ice growth (Fig. 7b). The values of these contributions within the tropics (30°S–30°N) are −0.18, −0.09, −0.06, +0.10, +0.02, and +0.0006 cm day−1, respectively. On the basis of the data shown in Fig. 7, one may infer that sea ice near the ice margins is controlled by two primary contributions, namely, dynamical transport from mid- and high latitudes and basal ice freezing (congelation ice growth) in local regions. The results also indicate that sea ice dynamics has important influences on the initiation of the Snowball Earth, as further addressed in section 4g.
c. The influences of reduced solar radiation as compared to reduced CO2 concentration
In general, the two types of experiments of reduced solar radiation and reduced CO2 concentration are much the same because both trigger very similar feedbacks (Hansen et al. 1984). The dissimilarity is that the solar radiative forcing is latitude dependent, whereas the greenhouse gas effect is nearly uniform (such as for CO2, CH4, and N2O, but not for O3). As demonstrated in Fig. 8b, in the two experiments with 94% solar radiation with 17.5 ppmv CO2 and 89.5% solar radiation with 286 ppmv CO2, both finally enter a hard Snowball Earth state but the sea ice evolution processes are quite different between them. Compared to the impact upon the sea ice evolution due to a reduction in solar radiation, sea ice in the reduced CO2 experiment first expands more rapidly and then more slowly when the sea ice fronts enter tropical latitudes. Although in terms of the global mean the radiative forcing is very close in the two experiments, −9.9 W m−2 for the experiment with 89.5% solar radiation with 286 ppmv CO2 and −10.6 W m−2 for the experiment with 94% solar radiation with 17.5 ppmv CO2 (relative to the experiment with 94% solar radiation and 286 ppmv CO2), in the tropics the reduction in the forcing associated with reduced solar radiation is larger than that associated with reduced CO2 concentration, while this situation reverses at high latitudes (Fig. 8a). In the calculation for Fig. 8a, the solar radiative forcing is defined as the incoming solar radiation at the top of the model multiplied by the planetary albedo, and the CO2 forcing (longwave plus shortwave) at the top of the model is −2.66 W m−2 for each halving its concentration (Collins et al. 2006a). The different influences of the solar radiation and CO2 forcing on the initiation of Snowball Earth can also be well understood in the context of a 1D energy balance model (see the appendix). In the energy balance model, it is clear that during the initiation of the Snowball Earth, the ice line is more sensitive in the early stage of the evolution and then less sensitive in the later stage (as sea ice line enters the tropics) to reduction in greenhouse gas forcing than to reduction in solar radiation (Fig. 9a).
Furthermore, as shown in Figs. 2 and 3, the bifurcation point (above which runaway glaciation occurs) in the experiments with reduction in CO2 concentration is higher than in the experiments with a reduction in solar forcing, 76% versus 49% sea ice coverage. A brief discussion based on the EBM is as follows. In the 1D EBM, the feedback factor is
where C is a parameter connected to the divergence of the net poleward energy flux due to the combined influences of the atmosphere and ocean, αp is the global mean albedo at the top of the atmosphere (TOA), B is a parameter related to the outgoing longwave radiation at the TOA (OLR), S is a normalized latitudinal distribution of solar insolation at TOA, αi is the albedo at the ice line, and xi is the location of the ice line (see the appendix; Roe and Baker 2010; Abbot et al. 2011). The value of f increases with decreasing latitude, and runaway glaciation occurs as f = 1 for both the reduction in the solar radiation and reduction in the greenhouse gas forcing (see the appendix). The higher the value of C, the more efficiently energy is redistributed from the equator to the poles, producing a greater cooling at the low latitudes and thereby promoting the formation of a globally ice-covered state, whereas a higher value of B represents a lower sensitivity of climate to perturbations (Roe and Baker 2010). The value of C depends on the atmospheric and oceanic circulations. In Part II of this study, it is found that there is essentially no difference in the tendency of poleward energy transport between the experiments with a reduction in CO2 concentration and with a reduction in solar radiation. The value of B depends on the tropospheric lapse rate, cloud cover, cloud-top temperature, water vapor amount, and greenhouse gases concentration (North et al. 1981). Possibly in the experiment with a reduction in CO2 concentration, the value of B will be higher than that in the experiment with a reduction in solar radiation since the atmosphere with less CO2 will allow more infrared radiation to be emitted to space if other feedbacks due to cloudiness, lapse rate, and water vapor variations may be assumed to be the same in these two types of experiments. A higher value of B enables f to reach 1 at a relatively low latitude, corresponding to a relatively high maximum stable ice coverage. As shown in Fig. 9b, an increase in B from 1.5 to 3.0 W m−2 °C allows the position of the maximum stable ice line to extend from approximately 35° to 20° latitude. Moreover, if earth’s obliquity is smaller than the present-day value (corresponding to a higher value of ∂S/∂xi at low latitudes), or if the albedo at the ice line αi is specified to have a smaller value (Abbot et al. 2011), the maximum stable ice line also will be closer to the equator. Therefore, even in an EBM, a soft Snowball solution exists. Further study is required to identify the appropriate values of the parameters B, C, and αi for Neoproterozoic climate.
d. Sea ice/snow albedo
For the purpose of testing the sensitivity of the results to surface albedo variations, we have performed a series of experiments with different sea ice and/or snow albedos. The results of these analyses are shown in Figs. 10 and 11. It is clear on physical grounds that increasing the sea ice and/or snow albedo will enhance the albedo feedback and promote the formation of a hard Snowball Earth. As shown in Fig. 3, the hard Snowball Earth occurs with 17.5–20 ppmv CO2 if sea ice albedo is set to 0.50 and snow albedo is set to 0.78. However, as the sea ice/snow albedo increases, global glaciation occurs with a higher value of the CO2 level, 286 ppmv for a sea ice albedo of 0.58–0.60 or 572 ppmv for a sea ice albedo of 0.65 (Figs. 10 and 11).
In comparing our results with the results obtained using ECHAM5/MPI-OM, it is important to note that their assumed sea ice and snow albedos are 0.75 and 0.80, respectively. When we assume the same albedos in CCSM3, the sea ice advances rapidly to produce a hard Snowball Earth (black dotted line in Fig. 10). It is therefore apparent that the primary difference between the CCSM3 and ECHAM5/MPI-OM predictions concerning the issue of Snowball Earth initiation may simply be associated with the assumed values for the albedos of snow and sea ice. However, it is also worth noting that differences in cloud forcing and atmospheric and oceanic heat transports among different models would also significantly contribute to the differences in forcings required to trigger Snowball Earth formation (Pierrehumbert et al. 2011; Voigt et al. 2011).
e. The influence of surface melting of snow and sea ice
To test the direct influence of melting snow on the surface of sea ice and of the direct surface melting of the sea ice itself, we have designed three simple experiments, as mentioned in section 3d. The results obtained in this set of experiments are displayed in Figs. 12 and 13.
As shown in Fig. 12a, the surface albedo has a strong seasonal cycle near the ice edge. From January to April the surface albedo is high (0.77–0.80); in May, the surface albedo begins to decrease; the albedo keeps declining until the end of August, and the average albedo for the month of August is 0.51; from fall to winter, the surface albedo increases again such that in December it reaches a value of 0.78. These variations correspond well to observations and to the results obtained in other simulations (Nazintsev 1964; Scharfen et al. 1987; Ebert and Curry 1993; Lindsay 1998; Perovich et al. 2002, 2003; Lüthje et al. 2006; Pedersen et al. 2009). For example, in the Arctic, Perovich et al. found that there were five distinct phases in the seasonal evolution of the albedo: for dry snow, the observed range is (~0.8–0.9); for melting snow, the albedo decreases from 0.8 to 0.7; for melt pond formation, a further decrease occurs from 0.7 to 0.5; as melt ponds expand, a further decrease occurs from 0.5 to 0.4; with fall freeze up, the albedo increases from 0.4 to 0.8 (Fig. 12a).
On the present earth, the length of the melting season is approximately 2 months in CCSM3 (gray lines in Figs. 12b–e), which is comparable to the observations; in the Arctic, the duration is ~60–80 days (Rigor et al. 2000; Perovich et al. 2002). As the sea ice enters into the tropics, however, the melt season near the ice edge may be as long as approximately four months (black lines in Figs. 12b–e). This is mostly because, as the sea ice moves to lower latitudes, the incoming shortwave radiation (ISW) at the surface is extremely high during most times of the year and never drops to below 120 W m−2, whereas in the Arctic, during half of the year, the ISW is very low, between 0 and 100 W m−2 (Fig. 12f). Moreover, near the ice edge, the lead area opening rate is comparable to that in the Arctic (not shown), which may also serve to lower the surface albedo.
If we neglect the temperature dependence of the albedo, corresponding to neglecting the effect of melting snow/ice, the earth would rapidly enter into a hard Snowball state (Fig. 13). Omitting the effect of sea ice melting would increase the surface albedo by ~0.15 near the ice edge in both hemispheres. In comparison, omitting the effect of snow melting will increase the surface albedo by ~0.20 in the Northern Hemisphere, but it nearly has no influence on the Southern Hemisphere because near the ice edge it is snow free (Figs. 13b and 13c). Both of these modifications serve to enhance ice/snow albedo feedback and thus promote the equatorward advance of the sea ice fronts.
We view these results as a demonstration of the plausibility of the idea that the formation of melt ponds and/or leads over the surface of sea ice during the melting season may be primarily responsible for preventing the runway albedo feedback from occurring, thereby inhibiting the formation of a hard Snowball. We intend to further test this idea with CCSM4, in which a more realistic parameterization for snow aging and melt ponds and a more accurate treatment for solar radiation have been implemented (Briegleb and Light 2007; Gent et al. 2011).
f. Sea ice vertical resolution and sea ice thickness distribution
As discussed in Abbot et al. (2010), present models overestimate the magnitude of the diurnal cycle of surface temperature by a factor of 2–3 owing to their low vertical resolution, which would underestimate the CO2 threshold for the deglaciation of the hard Snowball. Here, we test its effect on the initiation of the Snowball Earth by increasing the number of sea ice layers to either 12 or 60. It is found that increasing the number of layers enhances the albedo feedback (Fig. 14a), mostly through decreasing the surface temperature in daytime (Abbot et al. 2010). We have also found that increasing (decreasing) the number of sea ice thickness categories weakens (enhances) the surface albedo feedback; at least five categories being necessary to ensure the accuracy of the simulations (Fig. 14b), which is consistent with the results of Bitz et al. (2001).
g. The importance of sea ice dynamics
Driven by wind stresses and oceanic currents, the sea ice moves equatorward, cooling the atmosphere and surface ocean and/or increases the surface albedo at the margins. It is found that, when sea ice dynamics is excluded from the model, the sea ice advances significantly more slowly and the transition time from the present-day control climate to a state of global sea ice cover is postponed for about 105 years (Fig. 15a). Moreover, as illustrated in Fig. 15b, when solar radiation is 90% of the present level (0.90TSI0), the effect of sea ice dynamics is significant, increasing the global sea ice cover from approximately 66% to 73%. Therefore, sea ice dynamics has a significant impact on the initiation of the Snowball Earth, consistent with the results of Lewis et al. (2007).
Another notable result of these analyses is that, in the absence of sea ice dynamics, sea ice becomes much thicker than it does in the presence of this influence (~20 vs ~2 m), especially at mid and high latitudes, which is qualitatively consistent with the result of Lewis et al. (2007). However, this sea ice is still not thick enough to flow under its own weight in the manner proposed by Goodman and Pierrehumbert (2003). It should be noted that all of the experiments have been performed by “abruptly” decreasing the carbon dioxide concentration, which is clearly unrealistic since the time scale of the carbon cycle should be very long, during which the sea ice might grow to O(103 m) in thickness and might then be thick enough to “flow” (Goodman and Pierrehumbert 2003).
h. The influence of initial conditions
As suggested on the basis of EBM and GCM analyses (e.g., North 1975; Romanova et al. 2006), final equilibrium states of the climate system may be strongly dependent on the initial conditions from which the state is approached: sea ice cover, ocean temperature, air temperature, etc. We find this dependence upon initial conditions to be characteristic of the CCSM3 model in the analyses, based upon the assumption of different initial sea ice fractions (5%, 21%, 65%, 69%, 78% and 81%), accompanied by different air and ocean temperatures (Fig. 16). Under the same external forcing, 94% solar radiation and 286 ppmv CO2, the final sea ice cover can be 30% or 42%, depending on the initial ice fractions (Fig. 16a). Moreover, under 90% solar radiation and 286 ppmv CO2, the final sea ice fraction can be 50% or 72% (Fig. 16b). These results indicate that in the coupled ocean–atmosphere model CCSM3 there are multiple equilibria, consistent with the results of EBM-based analyses (e.g., North 1975).
Interestingly, even if the initial sea ice fraction is as high as 81%, the sea ice will retreat to midlatitudes [~40°S (N)], which is in some degree consistent with the results of the ECHAM5/MPI-OM of Voigt and Marotzke (2009). It should be mentioned that the experiments of both CCSM3 and ECHAM5/MPI-OM are initialized from various transient states rather than stable states. Actually, Lewis et al. (2004) found that as the University of Victoria EMIC model is initialized from a stable state with 80% sea ice, it may easily enter a hard Snowball Earth state.
i. The influence of variations in orbital geometry
It has been argued that the duration of the Neoproterozoic glaciations may have been as long as several million years (e.g., Hoffman et al. 1998; Macdonald et al. 2010) during which the earth’s eccentricity, obliquity, and precession would have evolved significantly over time. If the sea ice edge is within the tropical region, this evolution of the orbit may itself induce the sea ice to extend either equatorward or poleward and thus influence the robustness of the soft Snowball solution. Here, we test three different groups of parameters, namely, those corresponding to the mid-Holocene warm period at 6000 yr BP, to LGM at 21 000 yr BP, and to PEGI at 116 000 yr BP subsequent to which the onset of the most recent Late Quaternary glacial cycle began (Peltier 2007). The insolation at the top of the atmosphere for the modern earth and for the differences between these three periods and the modern earth are illustrated in Fig. 17. Our analyses of the impact of these differences in orbital configuration show that sea ice coverage is only slightly affected, less than 1% (Fig. 18). Our soft Snowball Earth solutions are therefore expected to remain stable in the presence of such perturbations.
j. Phase space plots
Based on the results of the above experiments, the dependence of the nature of the climate equilibrium on solar radiation, CO2 concentration, and sea ice albedo are summarized in Fig. 19. Apparently, between the partial-ice solution (present-day-like climate) and the globally ice-covered solution (hard Snowball Earth), there are intermediate and stable states (soft Snowball Earth) in which sea ice enters the tropics, snow covers low-latitude continents, and ice-free oceans persist. As the sea ice albedo is 0.50, soft Snowball solutions exist for 50–20 ppmv CO2 in addition to a 6% reduction in solar radiation (Fig. 3). While, as the sea ice albedo is increased to 0.55 (0.60), there remains a soft Snowball solution but with a higher CO2 level of 286 (572) ppmv.
As suggested by Pierrehumbert et al. (2011), a sufficiently low albedo is extremely important for the existence of a soft Snowball state; however, the sea ice albedo does not need to be as low as the value of 0.45 assumed in the EBM–ice sheet coupled model (Hyde et al. 2000; Peltier et al. 2007; Liu and Peltier 2010, 2011). Indeed, the sea ice albedo may be as high as 0.60–0.65 (Figs. 11 and 19), which is in reasonable accord with direct observations (Perovich 1996; Persson et al. 2002; Brandt et al. 2005), although it requires more greenhouse gases to maintain the state.
The suite of coupled ocean–atmosphere models previously employed to investigate Snowball Earth initiation prior to the current study included FOAM (Poulsen and Jacob 2004) and ECHAM5/MPI-OM (Voigt and Marotzke 2009; Voigt et al. 2011). In FOAM, the critical point for the hard Snowball onset corresponds to a 9% reduction in solar radiation if the CO2 concentration is 140 ppmv under boundary conditions consisting of an ideal, low-latitude supercontinent (Poulsen and Jacob 2004). This is in contrast to the result of ECHAM5/MPI-OM, which suggests that the critical point corresponds to a reduction of solar insolation by 6%–9% under modern geography and topography (Voigt and Marotzke 2009) or a 4%–4.5% reduction in solar radiation under the Marinoan geography (Voigt et al. 2011) with approximately 280 ppmv CO2. In our analyses based upon the CCSM3 model, the critical point corresponds to a 10%–10.5% reduction in solar constant with 286 ppmv CO2 under the modern geography and topography, which is close to that in FOAM but is a much more extreme requirement than that found to be necessary in ECHAM5/MPI-OM. As we have demonstrated in the aforementioned analyses, the major reason for the differences may be that the standard sea ice/snow albedo of CCSM3 (0.50, 0.78) is similar to that in FOAM [(0.55, 0.72), see Pierrehumbert et al. 2011] but is much smaller than that in ECHAM5/MPI-OM (0.75, 0.80). Differences in cloud radiative forcings and oceanic and atmospheric circulations also influence the results (Pierrehumbert et al. 2011; Part II).
To further reduce the uncertainty in these results, the detailed constraints on sea ice/snow albedo will require further investigation, especially near the sea ice edge, which is critical for determining whether the ice albedo instability will actually occur. On the hard Snowball Earth, the tropical ocean would be covered by sea glacier, and for this type of ice its albedo is found to be in the range of 0.55–0.66 (Warren and Brandt 2006; Dadic et al. 2010). On the soft Snowball Earth, the situation is more complex. The nature of the ice near the tropical sea ice front is a mixture of sea ice growing locally in place and sea ice exported from mid and high latitudes. For ice of this kind, a detailed understanding of its albedo has yet to be adequately investigated.
Importantly, the formation of melt ponds and leads during the melting seasons may significantly decrease surface albedo. With a simple temperature-dependent parameterization of melt ponds, we find that the sea ice/snow albedo near the ice edge will decrease from a high winter–spring value of 0.80 to a low summer–autumn value of ~0.51, and the melting season is significantly extended compared to that for high-latitude sea ice on the present earth. We infer that, when sea ice enters the tropics, the formation of melt ponds on the ice surface may be such that they are deeper and greater in areal extent and longer lived than those which form in the Arctic under the present-day conditions. More detailed snow/ice albedo parameterizations will be required to better understand the influence of melt ponds and leads.
There are, of course, many additional factors, which may either increase or decrease the sea ice/snow albedo, that have yet to be explicitly included in CCSM3 and other models. 1) The solar zenith angle (SZA) dependence of albedos for snow and sea ice has not been included, yet its influence is not expected to be negligible. As solar zenith angle increases, the snow and sea ice albedos will increase, and vice versa. For example, if the SZA increases from 40° to 80°, the snow albedo will increase by a few percent in the visible but by as much as 20% in the near infrared and the all-wavelength averaged albedo increases by about 0.08 (Wiscombe and Warren 1980). Also, SZA variations can explain up to 50% of the seasonal variability of the snow albedo in Greenland (Wang and Zender 2010). 2) As snow ages, with or without melting, the grain radius increases and the albedo is thereby reduced. The rate of decrease depends on the grain size and the shape of snow crystals. 3) Cloud cover can affect the sea ice/snow albedo by changing the spectral distribution and the effective zenith angle of shortwave radiation (Wiscombe and Warren 1980). 4) Synoptic weather events, such as rainfall, can temporarily decrease the albedo by ~0.05–0.15 (Perovich et al. 2002); general circulation models cannot accurately capture this influence owing to their low resolution. 5) Very low concentrations of highly absorbing impurities, such as black carbon or dust over sea ice/snow, may significantly lower the sea ice/snow albedo (Stephen and Warren 1980; Le-Hir et al. 2010; Abbot and Pierrehumbert 2010). Also, dust aerosol in the air may cause a significant warming over the sea ice/snow surface (Abbot and Halevy 2010). 6) The accurate parameterization of the solar spectrum is potentially very important. Using a broadband model, McKay (2000) suggested that tropical sea ice on a hard Snowball Earth would be on the order of 10 m thick, whereas using a spectral model with 60 wavelengths, Warren et al. (2002) found that the broadband model had overestimated the absorption depth of sunlight and the sea ice would actually be a few hundred meters thick. 7) The albedo of new snow is in the range of 0.75–0.87 and the albedo of melting snow ranges from 0.40 to 0.60; the albedo of sea glacier is in the range 0.55–0.66, and the bare sea ice albedo is in the range 0.47–0.52; the albedo of melt ponds varies from 0.15 for an old melt pond to 0.40 for a refrozen melt pond. Although the considerable variability of surface albedo as a function of surface type is readily apparent, it remains a challenge to determine the appropriate albedo to employ over a large area, which is inevitably >104 km2 in general circulation models (Perovich 1996; Warren et al. 2002; Warren and Brandt 2006).
By employing the fully coupled atmosphere–ocean–sea ice model CCSM3 to examine the conditions required for the initiation of either a soft Snowball or a hard Snowball state for the earth with modern geography and topography, we are led to the following conclusions.
In CCSM3, the forcings required for triggering a hard Snowball Earth is a 10%–10.5% reduction in solar radiation with 286 ppmv CO2 or a 6% reduction in solar radiation with 17.5–20 ppmv CO2. These thresholds are much more extreme than that in the ECHAM5/MPI-OM (Voigt and Marotzke 2009; Voigt et al. 2011). The explanation for these dramatic differences in required forcings is mainly the excessively high sea ice albedo employed in the ECHAM5/MPI-OM integrations (0.75 vs 0.50). Sea ice/snow albedo is critical to the results obtained in the Snowball Earth initiation experiments. Increased sea ice/snow albedo enhances albedo feedback and makes it easier to enter a globally ice-covered state.
In the range of CO2 concentration between 50 and 20 ppmv with the CCSM3 default sea ice albedo (0.50), we find soft Snowball Earth solutions. As the sea ice albedo is increased to 0.55–0.60, soft Snowball Earth solutions remain, but a CO2 level between 286 and 572 ppmv is required to maintain these states. In the soft Snowball solution, a belt of open water on the equator coexists with equatorial snow-covered continents and the sea ice margins may remain stable at 10°S and 10°N latitude, consistent with recent observations, which imply that some oceans should be ice free accompanying an active hydrological cycle during the Neoproterozoic glaciation events, as reviewed by Allen and Etienne (2008). We therefore conclude that the soft Snowball Earth solution proposed by Hyde et al. (2000) and Peltier et al. (2007) is entirely plausible.
The influences of solar radiation and CO2 forcing on the initiation of the Snowball Earth are not equivalent because the magnitude of solar radiation is latitude dependent, whereas the CO2 forcing is uniform. The sensitivity of the ice line to reduction in CO2 concentration is at first larger and then less (as sea ice enters the tropics) than that to reduction in solar radiation. Furthermore, the maximum stable sea ice coverage in the reduced CO2 experiments is greater than that in the reduced solar radiation experiments, 76% versus 49%.
The maximum stable sea ice coverage also depends on the sea ice/snow albedo. When solar radiation is reduced to 94% of the modern value, the maximum stable sea ice coverage is 76% for a sea ice albedo of 0.50. For an increased value of the sea ice albedo of 0.65, however, the maximum stable sea ice coverage declines to a value less than 65%.
The development of melt ponds and leads may significantly decrease the sea ice/snow albedo during the melting seasons and contributes significantly to the prevention of runaway albedo feedback, especially when the sea ice has extended into the tropics.
Increasing the vertical resolution of sea ice enhances the ice albedo feedback and thus promotes the equatorward extension of sea ice. However, increasing the sea ice thickness categories weakens the ice albedo feedback and thus inhibits the growth of sea ice; we find that at least five categories for the sea ice thickness distribution are necessary for the accuracy of the simulations.
Sea ice dynamics has an important effect on the initiation of the Snowball Earth. The sea ice dynamics accelerates the equatorward advance of sea ice fronts, and turning off the sea ice dynamics makes it harder to enter a hard Snowball state.
Multiple equilibria have been found to exist in the 3D coupled ocean–atmosphere–land–sea ice model CCSM3, similar to the characteristic of energy balance models.
Milankovitch-type variations in Earth’s orbit have very small influences on the equatorward extension of sea ice fronts. In particular, the soft Snowball solutions may not be destabilized by perturbations of this kind.
In Part II of the analyses, we examine further characteristic of the solutions delivered by this sequence of Neoproterozoic calculations, including details of various feedbacks, and find that the most important feedbacks during the initiation of the Snowball Earth are the albedo feedback and water vapor feedback. The most significant aspect of Neoproterozoic glacial climate that these papers do not address is the potentially important influence of continental configuration, which has recently been addressed in an EBM–ice sheet coupled model (Liu and Peltier 2010) and in a coupled ocean–atmosphere model (Voigt et al. 2011). Also important, of course, is the issue of the negative feedback of the carbon cycle, as has been discussed in Peltier et al. (2007) and Liu and Peltier (2011). These issues will be further addressed elsewhere.
We thank A. Voigt and D. S. Abbot for helpful and thoughtful reviews. We also thank Yonggang Liu for deep discussions. J. Yang and Y. Hu are supported by the National Basic Research Program of China (973 Program, 2010CB428606) and by NSF of China under Grant 40875042 and 41025018. J. Yang is partly supported by the Over-sea Study Program for Graduate Students of the China Scholarship Council. W. R. Peltier is supported by the Canadian Foundation for Climate and Atmospheric Sciences and a consortium of Canadian universities and by NSERC Discovery Grant A9627. The required computations were performed on the SciNet facility at the University of Toronto, which is a component of the Compute Canada HPC platform.
The Influences of Reduced Solar Radiation as Compared to Reduced Greenhouse Gas Concentration in a Budyko-Type EBM
where Q is the solar constant; x is the sine of latitude, x ∈ [0, 1]; S(x) describes the meridional distribution of solar insolation at the top of the atmosphere; T(x), , and α[T(x)] are the local temperature, the global-mean temperature, and the albedo as a function of temperature, respectively; A + BT(x) is a linearization of the outgoing longwave radiation; is a simple parameterization for the divergence of poleward heat transport by the atmosphere and ocean; A, B, and C are parameters (positive), and the temperature is in degrees Celsius.
The functional dependence of the albedo α[T(x)] is
where α1 and α2 are the albedo over ice-free and ice-covered surfaces (α1 < α2), respectively; αi is the albedo at the ice line (x = xi), simply taken to be (α1 + α2)/2; Ti is the temperature at the latitude of the ice line, generally taken to be −10°C. The standard approximation for S(x) is
where s2 = −0.482. Integrating Eq. (A1) from the equator to pole, one obtains the results for the global-mean energy balance:
where αp is the global-mean albedo,
and global-mean net solar radiation is
Considering Eq. (A1) at the ice line (x = xi) gives
First, we consider the sensitivity of the ice line to variations in the solar constant. In this calculation, A, B, and C are viewed as constants. Following Roe and Baker (2010), a first-order Taylor series expansion of Eq. (A8) gives
therefore, the sensitivity of the ice line to the solar constant can be written as
The sensitivity of the ice line to global-average net solar radiation can be written as
Next, we consider the sensitivity of the ice line to variations in greenhouse gas concentration. An increase in the value of A represents an increase in outgoing longwave radiation at fixed surface temperature, corresponding to a decrease in greenhouse gas concentration and an increase in ice coverage. In this calculation, the variation in A is taken to be the decrease from its modern value (ΔA = A0 − A) (Abbot et al. 2011), and the B, C, and Q are viewed as constants. Then, Eq. (A8) can be written as
Following Abbot et al., the sensitivity of the ice line to variations in greenhouse gas forcing can be written as
From Eqs. (A13) and (A17), it follows that the ratio of the sensitivity of the ice line to the reduction in global-mean net solar radiation (δxi/δQavg) to the sensitivity of the ice line to the reduction in greenhouse gas forcing [δxi/∂(ΔA)] is
when the ice line is at the pole, xi = 1, S(xi) = 0.518, and αp = α1 < αI; thus, β(xi) > 1, showing that the ice line is more sensitive to the reduction in greenhouse gas forcing than in solar radiation. When the ice line is at the equator, xi = 0, S(xi) = 1.241, αp = α2 > αI; thus, β(xi) < 1, meaning that the ice line is more sensitive to the reduction in solar radiation than in greenhouse gas forcing. The β(xi) as a function of the ice line (xi) is shown in Fig. 9a. It is clear that β(xi) > 1.0 when the ice line is at poleward of ~25°S(N) and β(xi) < 1.0 when the ice line enters the region between 25°S and 25°N. This suggests that, during the initiation of Snowball Earth, the sea ice advance rate in the experiment with a reduction in CO2 forcing would be faster in the early stage and slower in the later stage (as sea ice enters the tropics) than in the experiment with a reduction in solar radiation, consistent with the results in CCSM3 (Fig. 8b).
However, the value of the parameters B and C in these two types of experiments may be different (see section 4c for further discussion).