The initiation of Neoproterozoic “snowball” climates in CCSM3: the influence of paleocontinental configuration

We identify the “hard snowball” bifurcation point at which total sea-ice cover of the oceans is expected by employing the comprehensive coupled climate model CCSM3 (Community Climate System Model version 3) for two realistic Neoproterozoic continental configurations, namely a low-latitude configuration appropriate for the 720 Ma Sturtian glaciation and a higher southern latitude configuration reconstructed for 570 Ma but which has often been employed in the past to study the later 635 Ma Marinoan glaciation. Contrary to previous suggestions, we find that for the same total solar insolation (TSI) and atmospheric CO 2 concentration (pCO2), the 570 Ma continental configuration is characterized by colder climate than the 720 Ma continental configuration and enters the hard snowball state more easily on account of the following three factors: the higher effective albedo of the snow-covered land compared to that of sea ice, the more negative net cloud forcing near the ice front in the Northern Hemisphere (NH), and, more importantly, the more efficient sea-ice transport towards the Equator in the NH due to the absence of blockage by continents. Beside the paleogeography, we also find the optical depth of aerosol to have a significant influence on this important bifurcation point. When the high value (recommended by CCSM3 but demonstrated to be a significant overestimate) is employed, the critical values ofpCO2, beyond which a hard snowball will be realized, are between 80 and 90 ppmv (sea-ice fraction 55 %) and between 140 and 150 ppmv (sea-ice fraction 50 %) for the Sturtian and Marinoan continental configurations, respectively. However, if a lower value is employed that enables the model to approximately reproduce the present-day climate, then the critical values of pCO2 become 50–60 ppmv (sea-ice fraction 57 %) and 100–110 ppmv (sea-ice fraction 48 %) for the two continental configurations, respectively. All of these values are higher than previously obtained for the present-day geography (17–35 ppmv) using the same model, primarily due to the absence of vegetation, which increases the surface albedo, but are much lower than that obtained previously for the Marinoan continental configuration using the ECHAM5/MPI-OM model in its standard configuration (∼ 500 ppmv). However, when the sea-ice albedo in that model was reduced from 0.75 to a more appropriate value of 0.45, the criticalpCO2 becomes∼ 204 ppmv, closer to the values obtained here. Our results are similar to those obtained with the present-day geography (70–100 ppmv) when the most recent version of the NCAR model, CCSM4, was employed.

These states may be classified as either "hard" snowball (e.g., Kirschvink, 1992;Hoffman et al., 1998), "soft" snowball or slushball (e.g., Hyde et al., 2000;Peltier et al., 2007), or "thin-ice" snowball (e.g., Pollard and Kasting, 2005;Tziperman et al., 2012), respectively. In the hard snowball hypothesis, thick sea ice is assumed to have covered the entire surface of the oceans simultaneously with the continents being covered by thick ice sheets during the glaciation events ("Sturtian" at ∼ 720 Ma, "Marinoan" at ∼ 635 Ma and perhaps also the "Gaskiers" glaciation at ∼ 580 Ma). In the soft snowball/slushball hypothesis, open ocean remains in the low latitude equatorial region so that both hydrological and photosynthetic processes would have remained active during such episodes of glaciation. The thin-ice snowball hypothesis, as is suggested by its name, argues that thin (sea) ice might have prevailed over the low-latitude oceans, perhaps especially within marginal seas or isolated lakes during the glaciation events so that photosynthetic processes could persist, while the other parts of the planet were hard-snowball-like. In this sense, the thin-ice snowball is not a distinctly different climate state from the hard snowball, but in terms of biological activity it is identical to the slushball state. The fact that both photosynthetic and hydrological activity are strongly suspected to have continued throughout the Neoproterozoic interval has provided strong stimulus for the continuing effort to develop an alternative explanation to the hard snowball hypothesis for these deep glaciation events (Runnegar, 2000;Corsetti et al., 2006;Allen and Etienne, 2008). Please note that in all these three hypotheses, land ice is assumed to have covered all the continents excepting perhaps some isolated islands (e.g., Hyde et al., 2000;Peltier et al., 2007) since all three hypotheses have been advanced to explain the inferred very low-latitude glaciogenic deposits at sea level (Hoffman et al., 1998;Schrag, 2000, 2002) during the Neoproterozoic interval of time. Without this precondition, definition of a soft snowball becomes somewhat arbitrary if it is purely based on the fraction of total sea-ice area to ocean area or the position of sea-ice front.
Since the geological and geochemical evidence is not sufficient to determine the sea-ice state during any of the Neoproterozoic glaciations, the three hypotheses have been supported primarily on the basis of the different climate models employed by the different groups of proponents of the respective hypotheses. For example, Kirschvink (1992) proposed and Hoffman et al. (1998) have supported the hard snowball hypothesis based on the suggestion that the end state would have been governed by the action of runaway albedo feedback that obtains in the energy balance models (EBMs) of Budyko (1969) and Sellers (1969), which suggest that the Earth would rapidly enter a hard snowball state once the sea-ice edge extends below ∼ 30 • latitude. However, later analyses based upon the use of the EBM of North et al. (1983) coupled to the continental ice-sheet model of Peltier (1990, 1991) demonstrated that the dynamical action of the continental ice sheet supported the ex-istence of steady-state slushball solutions (Hyde et al., 2000;Crowley et al., 2001) as did models which coupled the influence of the continental ice sheets to an atmospheric general circulation model (AGCM) coupled in turn to a mixedlayer ocean model (Hyde et al., 2000). Subsequently, Pollard and Kasting (2005), and Tziperman et al. (2012) have suggested using an alternative sea-ice model (embodying the idea of sea glaciers) in which, even in a overall hard snowball state, thin ice could exist in some isolated regions. How stable such thin ice might be in such isolated regions is clearly a concern. A Jormungand state in which the sea ice enters the tropics (∼ 5 • latitude) without the Earth entering a hard snowball state has also been proposed based on AGCM simulations . The pCO 2 required to melt back from this state was found to be much higher than that required to form it; i.e., there was found to be a strong hysteresis associated with such a state, which was imagined to exempt it from suffering what is considered to be one of the major weaknesses of the soft snowball/slushball hypothesis by hard snowball Earth proponents (Schrag and Hoffman, 2001). However, such strong hysteresis is not found to exist in a fully coupled AOGCM model (e.g., Voigt and Abbot, 2012).
In order to settle the debate between the three primary competing hypotheses from a physical point of view, a model simulation which properly considers all of the relevant physical processes in the atmosphere, ocean, sea ice and land ice is required. The ideal model on the basis of which one might seek support for the soft snowball hypothesis, for example, would be an atmosphere-ocean general circulation model (AOGCM) coupled with an active land ice-sheet model, if this were able to adequately describe the initiation and development of global continental ice cover without leading to the ocean becoming completely ice covered. Unfortunately, such a model is still unavailable although both sophisticated AOGCMs and land ice sheet models do clearly exist. Nevertheless the very long timescale over which such a coupled structure would have to be integrated (∼ 100 kyr) associated with the formation of a snowball Earth (Donnadieu et al., 2003;Liu andPeltier, 2010, 2011) is not possible to contemplate on account of the computational resources that would be required. To adequately test the thin-sea-ice hypothesis, a sea-glacier model such as that developed in Tzipperman et al. (2012) may be required to replace the sea-ice module in the current AOGCMs as well as coupling to an explicit model of land ice-sheet evolution.
In the current state of research in this area, we are still several steps away from being able to reach the goal described above. Rather than applying a fully land-ice coupled AOGCM, only stand-alone AOGCMs have been applied (Poulsen et al., 2001(Poulsen et al., , 2002Donnadieu et al., 2004;Peltier et al., 2004;Poulsen and Jacob, 2004;Voigt and Marotzke, 2010;Voigt et al., 2011;Voigt and Abbot, 2012;Yang et al., 2012a, b, c) to search for the bifurcation point (in the space of e.g. sea-ice fraction versus radiative forcing; Voigt and Y. Liu et al.: The initiation of Neoproterozoic "snowball" climates in CCSM3 2557 , beyond which a hard snowball Earth, characterized solely in terms of total coverage of the oceans by sea ice, forms. Such models clearly cannot really be employed as a tool to argue definitively for or against the soft snowball hypothesis since a land ice sheet does not develop, although for specific continental configurations it may be used to argue whether it is possible for an ice sheet to develop before entering a hard snowball state by considering the existence of perennial snow cover over the continents . However, such models still represent an improvement over the AGCMs previously employed (e.g., Jenkins and Frakes, 1998;Jenkins and Smith, 1999;Chandler and Sohl, 2000) since they fully incorporate the influence of ocean dynamics.
The first relatively systematic study of the hard snowball bifurcation point using a reasonably modern AOGCM appears to have been that by Poulsen et al. (2001Poulsen et al. ( , 2002. The later resumption of such study using a more sophisticated AOGCM began with the work of Voigt and Marotzke (2010), who employed the present-day continental configuration as the basis for their analysis. Their results suggested that the Earth would enter a hard snowball state if the total solar irradiance (TSI) were reduced to somewhere between 91 and 94 % of the present-day value (TSI 0 , which is 1367 W m −2 ) and greenhouse gas concentrations were kept the same as their pre-industrial values. The model they employed was the ECHAM5/MPI-OM, and with this particular model, no nearsnowball slushball state was found to be stable; i.e., once sea-ice cover exceeded ∼ 57 % of the total ocean area, the climate quickly entered the hard snowball regime. However, with the same model, but for their more realistic Neoproterozoic continental configuration, Voigt et al. (2011) found that the Earth would enter a hard snowball state when TSI is reduced to between 0.955 TSI 0 and 0.96 TSI 0 , or if the TSI is fixed at 0.94TSI 0 (which is appropriate for the late Neoproterozoic; Gough, 1981), and the critical value of the atmospheric CO 2 concentration was between 500 and 556 ppmv, i.e., well above the modern level (see Voigt and Abbot, 2012). Again, no stable state was found in which the sea-ice front could enter within 25 • of the Equator and the system remain stable in a slushball state.
By employing a different coupled AOGCM, namely the National Center for Atmospheric Research (NCAR) Community Climate System Model version 3 (CCSM3), Yang et al. (2012a, b) demonstrated that for the present-day continental configuration, the critical values for both TSI and atmospheric CO 2 concentration (pCO 2 ) were much lower than those obtained using the ECHAM5/MPI-OM model, i.e., that it was much more difficult to form a hard snowball using this model. In particular it was found that the Earth would enter a hard snowball state only if the TSI were reduced to between 0.895 TSI 0 and 0.90 TSI 0 when pCO 2 is fixed to the pre-industrial value, or when pCO 2 is reduced to between 17.5 and 20 ppmv when TSI is fixed to 0.94 TSI 0 . As important, however, concerning their results is the fact that this model delivered stable slushball solutions in which the sea-ice front lay considerably equatorward of 25-30 • latitude. When a more advanced version of the model, CCSM4, was employed, these critical values at the bifurcation point were found to increase moderately to 0.91-0.92TSI 0 and 70-100 ppmv, respectively (Yang et al., 2012c). Therefore, the conditions required for the initiation of a hard snowball Earth have been found to be highly model dependent, the transition being much easier to achieve in the ECHAM5/MPI-OM than in either the CCSM3 or CCSM4 models. Yang et al. (2012a, b, c) identified the primary cause for the difference between the critical points predicted by the ECHAM5/MPI-OM and NCAR AOGCMs to be the seaice and snow albedos assumed in these models, with their values being highest in the ECHAM5/MPI-OM, slightly lower in CCSM4 and much lower in CCSM3. The most recent evaluation by Rogers et al. (2013) of AOGCM model performance in reproducing pan-Arctic sea-ice retreat between 1980 and 2008 has demonstrated that CCSM3 does better than ECHAM5. Follow-on analyses by Voigt and Abbot (2012) in which the sea-ice and snow albedos in ECHAM5/MPI-OM were reduced to values as low as those in CCSM3 led to the conclusion that it was still easier to initiate a hard snowball Earth than in either CCSM3 or CCSM4, while disabling sea-ice dynamics makes the initiation much more difficult to achieve. These authors therefore concluded that the sea-ice dynamics represented in ECHAM5/MPI-OM might be the major reason for easier initiation of hard snowball Earth than in either of the NCAR models. Pierrehumbert et al. (2011) have provided a useful review of the issues raised in these primary sources.
In the interpretation of the meaningfulness of these computations of the hard snowball bifurcation point, however, it will be important to keep in mind the fact that an additional negative feedback which acts through the agency of the ocean carbon cycle could entirely eliminate the existence of such a bifurcation point. Peltier et al. (2007) have suggested that as the climate system begins to cool towards a transition of this kind a strong drawdown of oxygen into the ocean would occur as a consequence of its increasing solubility as temperature decreases. Since the Neoproterozoic ocean is expected to be especially rich in dissolved organic carbon, the increased ventilation of the oceans would lead to the remineralization of this organic reservoir, the result being rapid production of carbon dioxide. Since part of this additional carbon dioxide would be exhausted into the atmosphere, the result would be a negative feedback through an enhancement of the greenhouse effect that would act to arrest the cooling towards the snowball bifurcation point and perhaps entirely eliminate the possibility of hard snowball Earth occurrence.
Our purpose herein is to fully test the influence of incorporating more realistic Neoproterozoic continental configurations for both the Sturtian and Marinoan glaciations on the bifurcation point for initiation of a hard snowball Earth climate state in CCSM3. This will provide an additional point of intercomparison with the results for the ECHAM5/MPI-OM 2558 Y. Liu et al.: The initiation of Neoproterozoic "snowball" climates in CCSM3 since Voigt et al. (2011) and Voigt and Abbot (2012) have obtained the bifurcation point for a continental configuration which they call Marinoan but which is more similar to the Sturtian continental configuration to be employed herein. This is a further step towards realistic snowball Earth simulations using the CCSM family of models. The goal is once more to find the greenhouse level (primarily CO 2 ) required to induce a snowball Earth when the TSI is fixed to be that appropriate for the Neoproterozoic, i.e., 0.94TSI 0 . The influence of aerosol on the results is also analyzed in the present paper for the first time and is also found to have nonnegligible influence on the outcome. Note that since we do not have explicit land ice (either prescribed or simulated) in the model employed in the present study, our results cannot be employed to comment in any definitive way on the plausibility of a soft snowball/slushball Earth.
For reasons that will be described in detail in Sect. 2.2, the continental configuration employed for the Sturtian glaciation in this paper is more similar to that employed by Voigt et al. (2011) for the Marinoan glaciation. This is unsurprising since according to the paleomagnetic reconstruction of Li et al. (2008), the continental configurations for the two glaciations are similar. However, because of limited data availability, some earlier studies employed the continental configurations constructed for 570 Ma for the Marinoan (e.g., Hyde et al., 2000;Peltier et al., 2004;Liu and Peltier, 2010, 2013a. We will continue to employ this continental configuration for the Marinoan glaciation in the analyses to follow because of its greater contrast regarding the latitudinal distribution of the supercontinental components with the 720 Ma continental configuration than the 630 Ma continental configuration (Li et al., 2008).
Before describing the modeling results for the more realistic Neoproterozoic continental configurations, we will first describe the results of sensitivity tests for the presentday geography. These tests will demonstrate how the climate may be influenced, for example, by removing vegetation from the land surface, simplifying the representation of aerosols, changing soil color, removing lakes and wetlands, etc. (see Sect. 3 for further detail), and will facilitate the understanding of the differences between the results obtained here and those in earlier studies of snowball Earth initiation with present-day geography (Yang et al., 2012a, b, c).

Model description
The Community Climate System Model version 3 employed in this work simulates four major components of the climate system, namely atmosphere, ocean, land surface processes and sea ice, and the respective model components are the Community Atmosphere Model version 3 (CAM3), the Parallel Ocean Program (POP, based on version 1.4.3), the Com-munity Land Model version 3 (CLM3) and the Community Sea Ice Model version 5 (CSIM5). These four components are linked through a flux coupler (CPL6), with no flux corrections between the components (Collins et al., 2006a).
Both CAM3 and POP are 3-D general circulation models, which will be run at the recommended (for paleoclimate studies) low resolution of T31 (approximately 3.75 • × 3.75 • ) and "g×3" (100 and 116 grid points in the zonal and meridional direction, respectively), and 26 levels and 25 levels in the vertical, respectively. CLM3, which considers dynamic vegetation (although not relevant here) and explicit river routing, will be run on the same grid as CAM3. CSIM5 includes both thermodynamics and dynamics of the sea ice and it is always run on the same grid as that of the ocean in the coupled mode. The time steps for the atmosphere, ocean, land and sea-ice components are 30, 80 (the default of 120 min was used in initial test runs with the present-day continental configuration; see below), 30 and 60 min, respectively. Fluxes are exchanged every 1 h between atmosphere, land and sea-ice components, and 1 day between the ocean and other components.

Continental configurations
The same continental configurations, one for 720 Ma (Sturtian) constructed by Li et al. (2008) but slightly modified to increase connectivity of the land masses and the other for 570 Ma (approximate Marinoan) constructed by Dalziel (1997), as employed by Liu andPeltier (2010, 2011) are adopted here (Fig. 1). There are three reasons why the 570 Ma continental configuration, rather than the 630 Ma configuration constructed by Li et al. (2008), is employed here. The first is that the 570 Ma continental configuration has been earlier employed by Hyde et al. (2000), who were the first to demonstrate the existence of soft snowball solutions in an AGCM (GENESIS 2). It will be useful to confirm the results obtained there by using an AOGCM, although the current study is just a first step towards that goal since no continental ice sheets will be included. The second is that the 630 Ma continental configuration of Li et al. (2008) is not significantly different from the 720 Ma continental configuration in terms of its influence on climate since continental fragments are concentrated in low latitude in both time slices, while the 570 Ma continental configuration provides a more distinct alternative. The final reason is that the Gaskiers glaciation, though most probably not global in scale, occurred around 580 Ma.
The mean elevation of the topographically "flat" continents relative to sea surface is assumed to be ∼ 400 m, also consistent with the assumption in Liu andPeltier (2010, 2011) which is lower than the present-day value (∼ 740 m) but slightly higher than that (100 m) employed by . It has been commonly accepted that the mean continental elevation had been almost constant through geological time (e.g., Hynes, 2001), but a recent study (Lorentz, Y. Liu et al.: The initiation of Neoproterozoic "snowball" climates in CCSM3 2559 2008) has challenged this view and proposed that the mean continental elevation could have been even higher during the Neoproterozoic. This issue is not expected to be settled anytime soon, but the most important thing relevant to the work here is that the difference in mean climate states obtained for 100 m and 400 m mean elevations is small, ∼ 0.5 • C in terms of globally averaged annual mean surface temperature (results not shown). To help create the river routing map needed for the analyses we have performed, the elevation at the center of the individual continents is raised by 50 m and decreases linearly to 400 m towards the edges/envelope of the continents.
To set up the CCSM3 model for these realistic Neoproterozoic continental configurations, we have generally followed the suggestions provided by the NCAR laboratory (Rosenbloom et al., 2011). The vegetation type on land is set to desert everywhere since surface vegetation has not yet significantly evolved by Neoproterozoic time. Soil color, which determines the dry and saturated soil albedo, is set to a medium value of 4 on a scale of 1 to 8, where 1 is the brightest (Oleson et al., 2004;technical report). The albedo of the soil (color 4) is 0.09-0.18 in the visible band and 0.18-0.36 in the near-infrared band depending on whether the soil is water saturated (corresponding to the low value) or dry (the high value). Soil texture/type which determines the soil thermal and hydrologic properties is set to be that of an average soil ("loam"), consisting of ∼ 15 % clay, 43 % sand and the rest silt. There are no lakes or wetlands assumed to exist on the land which would lower the effective albedo of the surface and therefore impede the transition into the hard snowball regime.
In order for the ocean module, POP, to run properly, the poles of the ocean grid have to reside on land. But because there are no polar continental fragments in either of the continental reconstructions we employ that lie poleward of ∼ 65 • for the 720 Ma continental configuration or in the vicinity of the North Pole for the 570 Ma continental configuration, we have elected to introduce "fictitious" circular continents of radius ∼ 5 • at the pole or poles as sites of the poles of the POP grid (Fig. 1). This strategy is preferred because (1) it involves the least stretch and distortion of the ocean grid, and (2) it is not expected to have a significant influence in determining the hard snowball bifurcation point since, long before entering a hard snowball state, the polar and mid-latitude regions will already be covered by sea ice, and therefore the ocean circulation within the polar region is expected to be weak. This has been partially demonstrated to be true in the analyses of modern snowball Earth initiation by Yang et al. (2012a, b, c).
To find the hard snowball bifurcation point, we need to have a control run for both continental configurations, a control run which has a relatively warm equilibrium climate. For this control run, we set the TSI to be 0.94TSI 0 , pCO 2 to be 2000 ppmv, aerosol optical depth τ to be 0.12, and pCH 4 and pN 2 O to have pre-industrial values 805.6 ppbv 720 Ma 570 Ma Fig. 1. The continental configuration for 720 Ma and 570 Ma, and their zonal land fractions. Grey and blue are land and ocean, respectively. Note that fictitious land has been added at the two poles. and 276.7 ppbv, respectively. Both the atmosphere and the ocean are initiated from static states, with the initial ocean temperature horizontally uniform but vertically decreasing from 24 • C at the surface to 9 • C in the abyss. The salinity of the ocean is assumed to be uniform with a value of 35 psu. However, negligible difference was found when the salinity was increased to 45 psu (results not shown). When the control run reaches equilibrium (after 4000 yr), branch runs are initiated by decreasing pCO 2 with TSI fixed.
Because the ocean bottom is flat and there is no continent in the mid-latitude of the Northern Hemisphere (NH) for either the 720 Ma or the 570 Ma continental configuration, the westerly jet of the ocean in mid-latitudes of the Northern Hemisphere will become extremely intense (like the Antarctic circumpolar current (ACC) of the modern ocean) and eventually the ocean model would fail to converge. This problem cannot be eliminated by decreasing the time steps of the components of CCSM3, nor by increasing the bottom drag (a factor of 100 was tested) of the ocean. Other options available to eliminate this problem include adding fictitious ocean bottom topography (referred to as the RIDGE method hereafter) or some narrow meridionally oriented continental strips in mid-latitudes, and increasing the horizontal tracer diffusivities of the ocean (referred to as the DIF method hereafter). The method we have chosen to employ, and for which the results are presented in what follows, is one based upon the increase of the horizontal tracer diffusivities. Comparison of the control run obtained by this method to that obtained through use of the RIDGE method ( Fig. 3 below) and detailed analyses on the mechanism that causes the difference (Appendix A) indicates that the two methods should give similar results insofar as the snowball Earth bifurcation point is concerned.
We increase both the isopycnal diffusivity (variable "ah" in the model) and the thickness diffusivity (variable "ah_bolus" in the model) by a factor of 5, i.e., from the default 800 m 2 s −1 to 4000 m 2 s −1 . This factor was initially determined by Y. Shin (University of Toronto, personal communication, 2010) to be sufficient to stabilize the model in the context of aqua-planet simulations. This and perhaps even higher values have been demonstrated to be appropriate for the regions where intense currents and thus intense baroclinic instability occurs, for example, around the region where intense westerlies are present (Ferreira et al., 2005). In the newer version of the NCAR model CCSM4, these two diffusivities have been formulated to be vertically varying depending on the stratification following Ferreira et al. (2005) and Danabasoglu and Marshall (2007). For the 3 • resolution employed in the analyses reported herein, both diffusivities can be as large as 4000 m 2 s −1 in the upper ocean, but decrease to 300 m 2 s −1 by a depth of about 2000 m. In our simulations with a flat ocean bottom, very strong westerly ocean currents (zonal speed reaches 70 cm s −1 and the width of the current in the meridional direction is approximately 20 • , not shown) develop in the NH, especially for the 570 Ma continental configuration, warranting the employment of higher tracer diffusivities than the default. A drawback is that this high tracer diffusivity is applied uniformly to the ocean in CCSM3, which may have a significant influence on the ocean circulation. However, as will be described in Appendix A, employing high tracer diffusivities should not have significant influence on the hard snowball Earth bifurcation point. Tracer diffusivities as high as 4000 m 2 s −1 have also been employed by Poulsen et al. (2001) in their study of ocean dynamics during snowball Earth occurrence.

Impact of vegetation and aerosol
Many of the variables in the model, such as vegetation, ozone, CFCs, aerosols, etc., are unconstrained or are known to have been very different from their present-day values during the Neoproterozoic, and must therefore be either simplified in their description or eliminated entirely. To test the degree to which the modification of these variables affects the climate, we choose to use the present-day continental configuration. The advantage of using a present-day continental configuration is that the CCSM3 model by default includes realistic specifications of the variables described above as well as others.
These tests turned out to be important because, as will be described in what follows, the continental vegetation and the nature of the atmospheric aerosol distribution have a significant influence on the climate through modification of the surface albedo and planetary albedo, respectively.
The series of simulations that were performed for the purpose of the above-described investigations is listed in Table 1 in order of the increasing number of variables being modified. The first simulation was performed using the default configuration of CCSM3 which is appropriate for the year 1990; for example, pCO 2 and pCH 4 are 355 ppmv and 1714 ppbv, respectively, significantly higher than those for the pre-industrial period, and TSI = TSI 0 . In the ensuing simulations 2-11, the following changes of variables are progressively introduced: ozone is changed from the distribution appropriate for present-day climatology to the pre-industrial distribution (Rosenbloom et al., 2011), vegetation is completely removed, lakes and wetlands are removed, soil color (which affects the albedo of soil; see above) is changed to a uniform value of 4, soil texture/type is changed to that for an average soil ("loam"; see above for composition. This determines the porocity, hydraulic conductivity, thermal conductivity and heat capacity, etc.; Oleson et al., 2004), pCO 2 is reduced to 300 ppmv, pCH 4 is reduced to 0, pN 2 O is reduced to 0, CFCs are removed, aerosol climatology is replaced by a time-invariable and spatially uniform layer which has an optical depth determined by the variable τ . This representation of the aerosol is the same as that employed in Community Climate Model 3 (CCM3), the immediate predecessor of CAM3. In this representation, the aerosol is assumed to be sulfate aerosol and to exist only in the lowest 3 levels of the model (Kiehl et al., 1996), and the optical properties (specific extinction, single scattering albedo and asymmetry parameter) for this aerosol are calculated as in Kiehl and Briegleb (1993). From Fig. 1 of Kiehl and Briegleb (1993), it can be seen that the sulfate aerosol particles do not absorb any radiative energy from the spectral band between 0.3 and 1.1 µm, but rather scatter this energy away (part of which is scattered back into space; see the asymmetry parameter in this figure). Therefore, it cools the Earth when its concentration increases (through increasing optical depth in the model).
The majority of these simulations were continued for only 100 yr in order to examine their short-term influence on the climate, and the results are shown in Fig. 2a. The time series of globally averaged annual mean surface temperature in this figure demonstrate that changing ozone, removing CFCs, and Table 1. Summary of runs carried out with the present-day continental configuration in order to test the influence of variable changes. Variable "τ " is the optical depth of aerosols. Characters "D", "PI", "R", and "-" stand for "default", "pre-industrial", "removed", and "not applicable", respectively, where the default value (in the brackets) may follow the character "D". "T " is resulting 10 yr mean global mean surface temperature at the end of simulations (either 100 yr or 2000 yr).
Runs Ozone Veg Lake Soil color Soil texture pCO 2 (ppmv) pCH 4 (ppbv)  ) and (b) show the short-term and extended runs, respectively. For part of the long runs 11, 13 and 14, only 10 yr of results were retained for every 100 yr period and these segments have been connected by straight lines in order to create continuous curve shown. Note that the legends describe only the further change of parameters relative to the previous run.
changing soil texture/type all have a slight (< 0.3 • C) cooling effect, while removing vegetation (run 3 in Fig. 2a) has a relatively significant cooling effect (∼ 0.5 • C in 100 yr).
Both removing lakes and wetlands (run 4) and changing soil color (run 5) after the vegetation has already been removed seem to have a negligible effect on the climate. A major impact on the short-term global mean temperature, however, derives from the reduction of pCO 2 , pCH 4 , pN 2 O and the change of the aerosol distribution, which cools the climate by ∼ 0.4 • C, 1.1 • C, 0.6 • C and 1.6 • C in 100 yr, respectively. We must therefore pay greater attention to what might be the most appropriate Neoproterozoic values for pCH 4 , pN 2 O and aerosol in our investigations of the hard snowball bifurcation point for the realistic continental configurations. In order to have a better idea as to what values might be appropriate, we need to consider the longer-term influence on climate.
To this end Fig. 2b shows 2000 yr time series for global mean surface temperature for runs 11, 12, 13 and 14. Clearly, if both pCH 4 and pN 2 O are removed and aerosol is replaced by a uniform layer of sulfate aerosol with its τ value (which is the optical depth of the layer) set to be 0.28 as recommended in the instruction manual for paleoclimate studies provided by the NCAR laboratory (Rosenbloom et al., 2011), the climate will be dramatically colder than that of the present day, with the global mean surface temperature as low as ∼ 4 • C (run 11, Fig. 2b). Since there are no constraints on these three variables for the Neoproterozoic, we will choose the preindustrial values 805.6 ppbv and 276.7 ppbv for pCH 4 and pN 2 O, respectively, to be consistent with Yang et al. (2012a, b, c). These are also the same values as those employed by Voigt and Marotzke (2010) but slightly higher than those (pCH 4 = 650 ppbv and pN 2 O = 270 ppbv) employed by Voigt et al. (2011) in the ECHAM5/MPI-OM simulations.
This will make comparison of the bifurcation points obtained here with those obtained in these previous studies more instructive.
For the choice of τ , the issue is more subtle. A value of 0.28 was suggested by Rosenbloom et al. (2011) for deep paleoclimate studies, but clearly adopting this value will cause the climate to cool by more than 2 • C in global mean surface temperature (compare the result for run 13, which uses a value of 0.24 for τ , with run 12, which employs the presentday aerosol climatology in Fig. 2b, all other conditions remaining the same. Using a higher value for τ such as 0.28 will cool the climate further). This value was preferred in an earlier version (CCM3) of CAM3, whereas in CAM3 or CCSM3 a spatio-temporally varying climatology for aerosol is used as an improvement. Heavens et al. (2012) performed a sensitivity study for the late Permian (251 Ma) with the CCSM4 model and found that the climate obtained with a uniform layer of aerosol with a τ value of 0.28 is ∼ 1.5 • C cooler than that obtained with a prescribed spatio-temporally varying aerosol. The global mean optical depth of their prescribed spatio-temporally varying aerosol, which itself was obtained by a prognostic model for aerosol, was only 0.08. Moreover, the global mean optical depth is only ∼ 0.17 and 0.08 for the present-day and pre-industrial climates, respectively (Mahowald et al., 2011). Therefore, we have good reason to believe the τ value of 0.28 is an overestimate of the optical depth appropriate for the Neoproterozoic climate simulations as well if the CCSM3 model is employed.
Run 12 (Fig. 2b), in which the pre-industrial pCH 4 and pN 2 O and default aerosol climatology are used, may serve as a reference simulation from which an appropriate value for τ can be determined. The global mean surface temperature obtained in this run is ∼ 11.2 • C, about 3.3 • C cooler than the present-day climate. The cooling is mainly due to the combined effects of removing vegetation, reducing pCH 4 , changing soil texture and removing CFCs as discussed above. We found that a close-to-present-day climate can be obtained if the value of τ is reduced to 0.12 (run 14 in Fig. 2b). Since there is no better constraint on this value, we will use both 0.28 and 0.12 in searching for the hard snowball bifurcation point for the Neoproterozoic continental configurations as a sensitivity test. Note that we may not completely rule out the high value because the Neoproterozoic climate is expected to be much drier than today due to its low temperature (see e.g. Fig. 4 of Yang et al., 2012b) and therefore the contribution of dust aerosol might have been more important than is the case today.

Predictions of the impact of continental configuration on the Neoproterozoic hard snowball bifurcation point
The control simulations for both the 720 Ma and 570 Ma continental configurations are shown in Fig. 3. The equilibrium state obtained by both the RIDGE and the DIF methods are close to each other. For the 570 Ma continental configuration, the difference in terms of global mean surface temperature is only ∼ 0.2 • C with the DIF method slightly warmer ( Fig. 3b; see also Table 2). The difference is larger for the 720 Ma continental configuration (the DIF method is ∼ 1.3 • C colder than the RIDGE method). However, the analyses presented in Appendix A indicate that the main influence of the DIF method on ocean circulation is that it readjusts the ocean heat transport so that the temperature of the two hemispheres becomes more symmetric. The analyses in Appendix A also indicate that the critical pCO 2 (below which a hard snowball Earth forms) obtained by the DIF method is approximately 10 % higher than that obtained by the RIDGE method, which is negligible considering the other uncertainties in configuring the model. Therefore, we will focus entirely in what follows on the results obtained by employing the DIF method. Table 2. Summary of the results of runs. T present is the annually averaged global mean surface temperature obtained at equilibrium under pre-industrial solar insolation and greenhouse gas concentrations, i.e., TSI=TSI 0 , pCO 2 = 286.2 ppmv. T control is similar to T present except that it is obtained for TSI = 0.94TSI 0 , pCO 2 = 2000 ppmv, i.e., the control simulations shown in Fig. 3. Critical pCO 2 indicates the value of pCO 2 for which the climate will enter a hard snowball state. This value is only determined to the precision of 10 ppmv. "RIDGE" indicates that a straight ridge of 2 km is added to the ocean floor (see Fig. A1 for their width and locations); "DIF" indicates that ocean bottom is flat but the bottom drag and tracer diffusivities are enhanced (see text).  Branch runs were initiated from the end of the control simulations shown in Fig. 3, and the results are shown in Figs. 4 and 5 for the 720 Ma and 570 Ma continental configurations, respectively. For these branch runs, the TSI is fixed at 0.94TSI 0 . As expected, the Earth enters a hard snowball state more easily than the present-day Earth as determined by Yang et al. (2012a, b) because of the cumulative cool-ing effect due to the removal of vegetation, lakes and wetlands, and changing soil color, etc., as described previously. When the optical depth of the aerosol layer is set to a relatively low value (τ = 0.12), the Earth enters a hard snowball state when pCO 2 is reduced to between 50 and 60 ppmv (Fig. 4a) and between 100 and 110 ppmv (Fig. 5a) for the 720 Ma and 570 Ma continental configurations, respectively. When the high value of 0.28 is employed, the Earth enters a hard snowball at higher values of pCO 2 between 80 and 90 ppmv (Fig. 4b) and between 140 and 150 ppmv (Fig. 5b) for the two continental configurations, an increase of ∼ 50 % and 40 %, respectively. Clearly, the influence of paleogeography on the hard snowball bifurcation point is significant, in fact comparable to that of the assumed optical depth of the aerosol layer. However, the influence of both may be considered small compared to the uncertainties in sea-ice albedo and sea-ice dynamics (e.g., Voigt and Abbot, 2012;Yang et al., 2012a, b, c).
The maximum sea-ice fraction (ratio of sea-ice area to the global ocean area) obtained before entering a hard snowball is ∼ 55 % (Fig. 4) and ∼ 50 % (Fig. 5) for the 720 Ma and 570 Ma continental configurations, respectively, much less than that obtained for the present-day continental configuration (∼ 75 %; Yang et al., 2012a, b) using the same model. This could be due to the increased ocean tracer diffusivities employed for the purpose of the analyses described in this paper. Our results for the maximum sea-ice fraction are similar to those obtained by Voigt et al. (2011) using the ECHAM5/MPI-OM model with the 635 Ma Marinoan continental configuration.
The sea-ice front (defined by the grid-box sea-ice fraction being equal to 50 %) generally reaches 24-27 • N(S) and 27-30 • N(S) for the 720 Ma and 570 Ma continental configurations, respectively (results not shown, but see Fig. 6 for a close analog). However, the local sea-ice front can reach as close as 15 • N for the 720 Ma continental configuration (around 180 • E in Fig. 6a) and 24 • S for the 570 Ma continental configuration (around 50 • W in Fig. 6b)  It may be noted that the critical pCO 2 values obtained for the 720 Ma continental configurations may be underestimated since the simulations with higher pCO 2 values have not reached statistical equilibrium 3000 yr after the branch run was initiated. For example, the run with pCO 2 = 60 ppmv is still drifting very slowly towards colder climate at year 7000 (Fig. 4) and we cannot therefore rule out the possibility that it might eventually enter a hard snowball state. The runs for which pCO 2 = 60 ppmv in Fig. 4a and pCO 2 = 100 ppmv in Fig. 4b have been continued till year 9000, and since the Earth had not entered a hard snowball state by this time the runs were not continued further.  . 6. The 100 yr mean sea-ice thickness (filled contours) and velocities (the arrows) at the end of simulations. The grey area represents the continent, and the hatched regions represent land snow thicker than 1 cm water equivalent. All sea ice with grid-box fraction greater than 1 % is plotted, but the contour lines (red curves) for grid-box fraction of 10 % and 50 % are also superimposed. In these two simulations, TSI = 0.94TSI 0 , pCO 2 = 140 ppmv, τ = 0.12. Figure 7 shows the relationship between global mean surface temperature as well as sea-ice fraction and pCO 2 . Clearly, the global mean surface temperature decreases almost linearly with applied CO 2 forcing (note the logarithmic scale of the x axis of Fig. 7) for both continental configurations. The trend is consistent with that found in Yang et al. (2012b) with the slope reaching approximately −6 K per halving of pCO 2 or ∼ 2.3 K (W m −2 ) −1 if the forcing per halving of pCO 2 is assumed to be −2.66 W m −2 (Collins et al., 2006b). This high climate sensitivity is beyond the upper bound of the climate sensitivity (2.1-4.4 K per doubling of pCO 2 , IPCC, 2007) obtained for the present-day climate among different models, probably due to the action of larger sea-ice and snow albedo feedback here because these characteristics of the model both extend into the tropical region. The relationship between sea-ice fraction and CO 2 forcing is also close to linear (Fig. 7b). Also seen from Fig. 7 is the fact that, subject to the same radiative forcings (solar insolation and greenhouse gas concentrations), the 570 Ma climate is always colder than the 720 Ma climate. The difference in global mean surface temperature is very small (∼ 0.5 • C) for the control run, but increases significantly to ∼ 3.5 • C when pCO 2 is reduced to concentrations equal to or below 286 ppmv. This is perhaps surprising since it was expected that the relatively high albedo of (bare) land would have a stronger cooling effect when the continents are located at low latitude rather than at high latitude (e.g., Kirschvink, 1992). Therefore, there are clearly other internal processes that are acting more strongly than the radiative forcings due to land albedo, and it is instructive to identify these processes.

Influence of emissivity, albedo and heat transport
For the control run, we find that the 570 Ma climate is only seasonally colder than the 720 Ma climate, especially during the austral winter (see Fig. 8a and b). Clearly the southern part of the Southern Hemisphere (SH) of the 570 Ma continental configuration has a much stronger seasonal cycle. This very strong seasonal cycle can be explained by the lower heat capacity of land, but it cannot explain the lower annual mean temperature in this region. To identify the processes that contribute to the lower annual mean surface temperature of the 570 Ma climate, we utilize the same 1-D EBM as in Heine- ( where Q (φ), α (φ), H (φ) and (φ) are respectively the incoming solar insolation at the top of atmosphere, the planetary albedo, divergence of total (atmosphere + ocean) meridional heat transport and atmospheric emissivity. The latter three of these quantities can be diagnosed, respectively, as the ratio of upward to downward shortwave radiative flux at the top of the atmosphere, the net of shortwave plus longwave radiative fluxes (upward is positive) at the top of the atmosphere, and the ratio of upward longwave radiative flux at the top of the atmosphere to that at the bottom of the atmosphere. All of these quantities (including Q (φ)) are calculated from the zonal mean radiative fluxes averaged over the last 100 yr of the simulations being analyzed. σ = 5.67 × 10 −8 W m −2 K −4 is the Stefan-Boltzmann constant. The zonal mean surface temperature T s calculated from Eq. (1) is almost identical to that obtained directly from the model in an equilibrium state (not shown). Therefore, this tool is very useful in analyzing the contribution of each parameter (α (φ), H (φ) or (φ)) to the difference in surface temperature between two simulations at each latitude.
The results from the 1-D EBM analyses are shown in Fig. 9. Please note that the temperature difference in Fig. 9 is weighted by the ratio of the surface area of a given latitudinal band (widths of 1 • ) to the total Earth surface area, so that the difference in global mean surface temperature will just be the sum of all the values at each latitude. In this way, we can directly see how much the temperature difference at each latitude contributes to the global mean temperature difference. From Fig. 9a (black curve), it is clear that the cooler global mean surface temperature of the 570 Ma climate relative to that of the 720 Ma climate is primarily due to the much colder temperature in the southern part of the SH. This colder temperature is mainly due to the higher emissivity (low water vapor effect), and higher albedo of the 570 Ma continental configuration in this region (Fig. 9a). The lower total atmosphere plus ocean heat transport for the 570 Ma continental configuration only contributes slightly in the south polar region. The reason as to why the emissivity, albedo and heat transport are different between the two continental configurations will be described in greater detail below for the branch runs.
For the branch runs (taking pCO 2 = 140 ppmv as an example because this value has been tested for both continental configurations), however, the surface temperature for the 570 Ma model is lower than that for the 720 Ma model in all seasons ( Fig. 8c and d). Moreover, the 570 Ma climate is colder than the 720 Ma climate at almost all latitudes except between 40 and 50 • N, but most significantly for latitudes < 35 • N (see black curve in Fig. 9b, and also see Fig. 10a for absolute values). For the southern part of the SH (> 45 • S), all three processes discussed above contribute to the colder surface temperature of the 570 Ma cli-mate, among which the emissivity contributes the most and planetary albedo contributes a similar effect to that of the heat transport (Fig. 9b). The 570 Ma atmosphere has higher emissivity due to its lower water content (Fig. 10d) and smaller longwave cloud forcing (Fig. 10f). It also has higher planetary albedo (Fig. 10c), which may be due to its higher surface albedo (Fig. 10b).
The difference in surface albedo between the 570 Ma and 720 Ma models is most significant during austral summer ( Fig. 11a and e). This is a consequence of the fact that, during summer, both grid-box sea-ice fraction and the snow depth on sea ice decrease over most of the SH. For example, at 60 • S where the contribution of albedo to temperature peaks locally, the grid-box sea-ice fraction reduces from 100 % in winter to 94 % ( Fig. 11b and f) and the snow depth is reduced from 13.6 cm to 7.5 cm ( Fig. 11c and g), with the latter translating into a reduction of snow coverage from ∼ 87 % to 79 % (Briegleb et al., 2004). These two effects therefore contribute more than 10 % of change in seasonal snow cover over that region. In contrast, the snow depth on land from 90 • S to 40 • S is very thick and does not have any significant seasonal variation ( Fig. 11d and h). Since it is mostly covered by land around 60 • S in the 570 Ma model (see the zonal fraction of land in Fig. 1), the surface albedo changes very little with season as shown in Fig. 11. This difference in effective snow albedo on land and on sea ice is smaller when pCO 2 is higher, e.g., 2000 ppmv, because the snow depth on land will have as significant seasonal variations (from over 80 cm during austral winter to less than 1 cm during summer, not shown) as the sea ice. This is most probably the reason why albedo contributes less to the surface temperature difference over the southern SH between the 570 Ma and 720 Ma models in the control runs (Fig. 9a).
For the southern SH (> 45 • S), the total (atmosphere + ocean) heat transport in the 570 Ma model is smaller than that in the 720 Ma model primarily due to the smaller ocean heat transport (not shown), which follows from the fact that the former model has reduced ocean area as compared to the latter.
For the tropical and subtropical regions (35 • S-35 • N), the reasons are more complex. Around the sea-ice edges (∼ 30 • S and ∼ 30 • N), both albedo and emissivity favor a cooler 570 Ma climate, especially the former. The higher planetary albedo in the 570 Ma model is most likely due to the more equatorward sea-ice edges (see Fig. 11b and f), and land snow cover in the SH (Fig. 11d and h) and the more lowlevel cloud cover in the NH. However the total heat transport (green curve in Fig. 9b) largely cancels the effect of albedo, and this is primarily due to the negative feedback of ocean circulation consistent with what has been found in previous studies (Poulsen and Jacob, 2004;Voigt and Abbot, 2012;Yang et al., 2012a, b). The reason that the 570 Ma model has more low-level cloud cover around 30 • N is very likely due to its larger ocean area, since water vapor is more available (note the significantly reduced low-level cloud cover over land in Fig. 12). This low-level cloud cover contrast between the 570 Ma and 720 Ma continental configurations around 30 • N (Fig. 12c) is also consistent with the shortwave cloud forcing in this region as shown in Fig. 10e. Note that because the shortwave cloud forcing is dependent on the surface albedo, the negative peak of this forcing is always more equatorward of the sea-ice edges as in Fig. 10e. This negative shortwave cloud forcing oceanward of the sea-ice edges, we think, provides a positive feedback for sea-ice expansion, and this feedback is stronger when there is more ocean around the sea-ice edge. This feedback is as strong as the negative feedback due to atmosphere plus ocean heat transport (see Fig. 9b).
Although cloud parameterizations are very different among different GCMs (e.g. Wyant et al., 2006;Zhao, 2013), this feature of higher low-cloud cover over ocean is physically sound and consistent with observed low-cloud cover for the modern climate (e.g., Warren et al., 2013). Different GCMs should generally reproduce this feature of the ob-served low-cloud cover but may differ dramatically in the magnitude of the radiative forcing. Therefore, we expect that the sign of the low-cloud forcing near the sea-ice edges obtained here to be reproducible in other GCMs, but the magnitude of the forcing may be either smaller or larger than obtained in the present study.
Over the tropical region, the reasons are more mixed. Since it has a relatively small contribution to the global mean surface temperature difference, we do not intend to analyze the reasons in greater detail here. However, it is noticeable that the surface albedo of the 720 Ma model is only slightly higher than that of the 570 Ma model because the soil is close to being water saturated, for which the soil albedo is only about half that when the soil is dry (see above). By comparing Fig. 10c to b, it can be inferred that the planetary albedo in this region is dominated by the cloud albedo. Therefore, the albedo effect of the more equatorial distribution of land in the 720 Ma model does not have any significant influence on the climate for the model settings presented and employed herein.
In summary, the 570 Ma climate is colder than the 720 Ma climate primarily because of the higher emissivity and higher planetary albedo of the former. The stronger total heat transport (mostly due to the ocean) towards the ice edges in the 570 Ma model is very important in canceling the sea-ice albedo feedback, but clearly contributes to the slight cooling over the tropical region (Fig. 9b). The difference in planetary albedo of the southern part of the SH is due to significantly larger seasonal variations of sea ice (and the snow depth/cover on it) than that of snow on land, while the difference around 30 • N is partially due to the higher low-cloud cover over oceans than over land. The different emissivities should be due to the different water content in the atmosphere for the two continental configurations; however, it is not absolutely clear whether this should be considered as a pure feedback mechanism or a partial triggering mechanism especially in the southern part of the SH.

Enhanced equatorward sea-ice transportation
The annual mean sea-ice edges in the 570 Ma model are closer to the Equator (by ∼ 5 • ) than those in the 720 Ma model within regions far from the continents (Fig. 13a; see also Fig. 6). To identify the reasons for this, all the quantities shown in Fig. 13 are zonally averaged over only the longitude belt 90 • W-30 • E, i.e., over the portion of the ocean outside of the 720 Ma continents. Since the SH of the 570 Ma continental configuration does not have ocean within this longitudinal range, some of the quantities are zero there and we will focus on the NH only. Sea ice is thicker for all latitudes in the 570 Ma model than that in the 720 Ma model (Fig. 13a), especially in the mid-latitudes. Clearly, the more rapid thermal growth of sea ice of the 570 Ma model (Fig. 13c)  due to the colder shallow ocean temperature of the 570 Ma model. For example, the ocean temperature at 200 m depth is approximately −1.8 • C and −1.2 • C for the 570 Ma and 720 Ma models, respectively. However, for the region between the ice fronts (∼ 21 • N and ∼ 26 • N for the 570 Ma model and 720 Ma model, respectively; see Fig. 13a) and 30 • N, the thicker sea ice of the 570 Ma model is obviously due to larger ice transport from higher latitudes (Fig. 13d). This larger transport is related to both thicker sea ice and larger meridional ice velocity between 25 and 50 • N (Fig. 13f), which leads to larger divergence both north of 30 • N (positive) and south of 30 • N (negative). We may consider that the portion of enhanced ice transport towards low latitude due to thicker sea ice in the high latitude is the result of the colder climate, while the portion due to large meridional ice velocity (and divergence) is not. Instead, it is clearly connected to the continental configuration.
First, it is clear on the basis of Fig. 13g and h that the wind speed in both the 570 Ma and 720 Ma climates are similar. The wind stresses are also similar with those of the 570 Ma model being slightly larger (not shown). Moreover, the direction of the meridional wind stress is opposite to the ice flow, meaning that the Hadley circulation is not playing a role. So there is no atmospheric cause of the enhanced ice velocities of the 570 Ma model. Since the wind stresses are similar, we may safely conclude on the basis of Fig. 6 that the much smaller zonal velocity of sea ice in the 720 Ma model (Fig. 13e) is due to the presence of continents in the subtropics. However, the smaller meridional velocity of sea ice cannot be directly related to the continents since the values shown in Fig. 13 are for the bulk ocean area only; rather, it is due to the larger southward Coriolis force associated with the larger zonal velocity of sea ice. This is shown to be the case in Fig. 13b. The difference in Coriolis stresses (in the meridional direction) between the 720 Ma and 570 Ma models is at least one order of magnitude larger than the difference in wind stresses and ocean stresses (not shown).
Therefore, in the absence of continents and shallow ocean bottom topography, exceptionally strong zonal sea-ice (and ocean) currents will develop under the action of westerly wind forcing. This strong zonal current promotes sea-ice transport (Ekman flow) towards the Equator due to the larger Coriolis force. This means that the influence of the continental configuration on sea-ice flow is another factor that may drive the climate towards a colder state (if there is no continental mass in the latitudinal bands of the westerlies). To quantify the strength of the influence of this cooling mechanism on climate, further simulations will need to be performed in which the influence of all other processes should be excluded. Note that this equatorward Ekman flow of sea ice becomes significant only when the climate is cold enough so that annual mean sea-ice front enters the westerlies zone (< 45 • S or 45 • N, Fig. 6).

Comparison with the results of previous publications
Our finding that high-latitude continents are more conducive to global glaciation contradicts the conclusion by Lewis et al. (2003) and Voigt et al. (2011). Voigt et al. (2011) came January July Fig. 11. Monthly mean surface albedo, grid-box sea-ice fraction, snow depth on sea ice and snow depth on land during January (left column) and July (right column) for both the 720 Ma and 570 Ma continental configurations. Note that the snow depth on sea ice is shown only up to 20 cm in order to make its seasonal variation around the 60 • S latitude more evident. The results are averaged over the last 100 yr of the simulations; TSI = 0.94TSI 0 , τ = 0.12, pCO 2 = 140 ppmv for both of the simulations.
to their conclusion based on the simulation of results for a single Neoproterozoic continental configuration and compared it with the results for present-day geography (Voigt and Marotzke, 2010); however they did not quantify the cooling effect of removing vegetation. Moreover, the case of the high-latitude continental configuration (for the present day) which they considered obviously cannot develop strong zonal ocean/sea-ice currents under the westerlies in either hemisphere, thereby missing a probably important cooling mechanism as discussed above. Furthermore, they have specified the albedo of land surface to be 0.272 uniformly, which is equivalent to the mean albedo for dry soil (0.18 for the visible and 0.36 for the infrared) employed here. For water-saturated soil, however, our albedo is reduced by half. Given the non-negligible precipitation on tropical land (not shown), the effective land surface albedo within 10 • of the Equator in our model is only slightly above 0.18, probably much lower than that em-ployed by Voigt et al. (2011). Moreover, the sea-ice albedo in the ECHAM5/MPI-OM model employed by Voigt and Marotzke (2010) and Voigt et al. (2011) is much higher than that in the CCSM3 model (Yang et al., 2012b), which may reduce the albedo contrast between snow-covered sea ice and snow on land as observed in Fig. 10b. However, our finding of the cooling effect of land-snow albedo for the high-latitude continental configuration is consistent with that by Poulsen et al. (2002), who employed the AOGCM FOAM1.4, even though the sea-ice albedo in FOAM (∼ 0.6) is higher than that (∼ 0.5) in CCSM3.
The model employed by Lewis et al. (2003) is

Conclusions
We have studied in detail the hard snowball bifurcation point for two realistic Neoproterozoic continental configurations using CCSM3, which is a further improvement on this subject relative to the work of Yang et al. (2012a, b), who employed the present-day continental configuration. The bifurcation point is dependent equally strongly on both the optical depth (τ ) of aerosol and the continental configuration. When a larger value for τ , 0.28, is adopted, the critical values of pCO 2 , below which a hard snowball will be realized, are be-tween 80 and 90 ppmv and between 140 and 150 ppmv for the 720 Ma and 570 Ma continental configurations, respectively. However, with a lower value, 0.12, which approximately reproduces the present-day climate (in terms of global mean surface temperature), the critical values of pCO 2 become 50-60 ppmv and 100-110 ppmv for the two continental configurations, respectively. For either case, they are higher than the values obtained by Yang et al. (2012a, b) probably due to the absence of vegetation during the Neoproterozoic and the consequently higher surface albedo of continental surfaces. These values are nevertheless much lower than those obtained using the ECHAM5/MPI-OM model in its standard configuration (∼ 500 ppmv; Voigt et al., 2011), and still somewhat lower than the values obtained when the seaice albedo in that model is reduced from 0.75 to a more appropriate value of 0.45 (∼ 204 ppmv;Voigt and Abbot, 2012) for the reasons that we have previously explained in Yang et al. (2012a, b, c). The remaining difference may be due to the significant difference in sea-ice dynamics between the two models (Voigt and Abbot, 2012), but the present study suggests that this could also be due in part to the different parameterization of aerosols. These new results are similar to those obtained previously using CCSM4 but with the present-day continents (70-100 ppmv) (Yang et al., 2012c). The high-latitude continental configuration (570 Ma) is found to be more conducive to global glaciation than the low-latitude continent distribution (720 Ma) in the CCSM3 model. Three factors contribute to the cooler climate of the 570 Ma configuration: in the mid-to high latitudes, a more significant seasonal variation of sea ice (and of the snow on it) than snow on land favors a warmer 720 Ma climate; around the ice edges (∼ 30 • ) in the near snowball state, the higher low-cloud coverage and thus more negative cloud forcing when there is more ocean within this region (for example, the NH of 570 Ma continental configuration) favors a colder 570 Ma climate; and finally the enhanced sea-ice transport towards the Equator when there is no continent blocking the eastwardly directed ocean/sea-ice currents under the westerlies (again, the NH of 570 Ma continental configuration) also favors a colder 570 Ma climate. An additional fourth factor that may not be directly related to the continental configurations is the greenhouse effect of water vapor (affecting the emissivity of the atmosphere). However this factor could be considered to simply constitute a pure positive feedback.
Among the three factors, the first one is probably model dependent and may not be constrained by present-day observations as it becomes significant only when the climate is very cold. The second factor should be robust across different GCMs since it is consistent with the observed present-day contrast in low-cloud coverage over land and over ocean, and we cannot identify any obvious reason why this characteristic should change in different continental configurations or different climate states as long as open ocean continues to exist. However, the magnitude of the forcing due to this low-cloud coverage may vary dramatically among different GCMs. We consider the third factor to be robust across GCMs as well since it is physically straightforward to deduce. However, how much this mechanism contributes to the cooler climate of a high-latitude continental configuration compared to that of a more equatorial continental configuration has not been quantified in the present work.

Stabilizing the ocean and the influence of increased horizontal thermal diffusivities
To stabilize the ocean currents which are far too fast with a flat-ocean-bottom and continent-free (in the Northern Hemisphere here) world, one option would be to add narrow meridionally oriented continents in mid-latitudes. However, this is not a preferred solution because this would be ex-pected to influence the climate significantly insofar as seaice transport is concerned and would make the control model of little use for the later analyses we will be performing in which explicit continental ice sheets are introduced.
We are then left with the possibilities of adding fictitious mid-ocean ridges or increasing tracer diffusivities. We first tested the response of the model to the addition of such fictitious ridges (referred to as the RIDGE method hereafter), examples of which are shown in Fig. A1. Both straight ridges of different height (Fig. A1a and c) and more complex presentday mid-ocean ridges (Fig. A1b) but with the ridges of the current Southern Hemisphere flipped to the north have been investigated. It was found that the addition of such ridges does postpone the onset of model instability, and the effect is stronger with higher ridges (results not shown). When the height of the straight ridge is 2 km, we were able to obtain a control run that reaches equilibrium. However, when branching off from this control run to obtain a colder climate, the model again failed to converge. The same occurred in the version of the model with the more complex set of ridges. We did not further increase the height of the straight ridge because, although the mid-ocean ridges of the present-day ocean can be higher than 2 km, the distance between the top of the ridges and the ocean surface is usually approximately 3 km or more (Fig. A1d). Even considering the sea level drop of 800-1000 m due to land ice sheet formation during snowball Earth occurrence (Liu and Peltier, 2013a, b), the distance is still around 2 km. Therefore, the method upon which we have finally settled is to increase the horizontal tracer diffusivities of the ocean (referred to as DIF method below). More specifically, we increase both the isopycnal diffusivity (variable "ah" in the model) and the thickness diffusivity (variable "ah_bolus" in the model) by a factor of 5.
Since the default values for both diffusivities in CCSM3 at g×3 resolution are 800 m 2 s −1 , the enhanced diffusivities in the DIF method here are then 4000 m 2 s −1 . This and perhaps even higher values have been demonstrated to be appropriate for the regions where intense currents and thus intense baroclinic instability occur, for example, around the region where intense westerlies are present under modern climate conditions (Ferreira et al., 2005). In the more recent version of the NCAR model CCSM4, these two diffusivities have been formulated to be vertically varying depending on the stratification following Ferreira et al. (2005) and Danabasoglu and Marshall (2007). For the 3 • resolution employed in the analyses reported herein, both diffusivities can be as large as 4000 m 2 s −1 in the upper ocean, but decrease to 300 m 2 s −1 by a depth of about 2000 m. In our simulations with a flat ocean bottom, very strong westerly ocean currents (zonal speed reaches 70 cm s −1 and the width of the current in the meridional direction is approximately 20 • , not shown) develop in the NH, especially for the 570 Ma continental configuration, warranting the implementation of higher tracer diffusivities than the default. A drawback is that this high tracer diffusivity is applied uniformly to the ocean in CCSM3 which may have a significant influence on the ocean circulation. However, as will be described in what follows, employing high tracer diffusivities should not have significant influence on the hard snowball Earth bifurcation point. Tracer diffusivities as high as 4000 m 2 s −1 have also been employed by Poulsen et al. (2001) in their study of ocean dynamics during snowball Earth occurrence.
In Fig. A2, the sea-ice thickness and velocity at equilibrium are shown for both methods for both continental configurations. Clearly the most significant difference between the two methods is that the sea-ice extent in the SH is much larger for the DIF method, which causes the surface temperature of the SH to be lower ( Fig. A3a and d). However, the surface temperature of the NH is higher for the DIF method especially for the 570 Ma continental configuration (Fig. A3d). The reason why the DIF method acts asymmetrically on the two hemispheres is connected to its influence on the ocean circulation and associated ocean heat transport. Ocean heat transport shows large changes ( Fig. A3c and f) when switching from the RIDGE method to the DIF method, and the sign is consistent with that of the surface temperature changes, decreasing in the SH and increasing in the NH. The atmospheric heat transport shows the opposite changes to counter the change in ocean heat transport ( The grey area is the continent, and the hatched regions represent land snow thicker than 1 cm. All sea ice with grid-box fraction greater than 1 % is plotted, but the contour lines (red curves) for grid-box fraction of 10 % and 50 % are also superimposed. The upper and lower panels were obtained by RIDGE method and DIF method, respectively. Note that the maximum value for the 570 Ma continental configuration (c and d) is much larger than that for the 720 Ma continental configuration (a and b) due to eddy-induced circulation (not shown). TSI = 0.94TSI 0 , τ = 0.12, pCO 2 = 2000 ppmv for all four simulations. e). However, clearly the change of ocean heat transport has a much bigger influence in the south polar regions. The change of ocean heat transport is naturally linked to the change of the meridional overturning circulation (MOC). Figure A4a and c show the MOC in the RIDGE method for the 720 Ma and 570 Ma continental configuration, respectively. Inspection will show that the wind-driven subtropical cells (STCs) of the upper ocean are similar to the present-day ocean. A strong channel flow develops at the middle latitudes of the NH for both these configurations of the model. Meanwhile, there exists a strong, counterclockwise deep-water cell from the mid-latitudes of the NH all the way to the high latitudes of the SH. Such a feature is similar to the present-day Earth and to that obtained in the Double Drake experiment of Ferreira et al. (2010) but with the NH and SH switched. It is clearly due to the so-called "Drake passage effect" (Toggweiler and Samuels, 1995;Toggweiler and Bjornsson, 2000), which is such that a zonal channel flow in one hemisphere enhances the deep-water formation in the other hemisphere. The effect is such that more heat is transported by the ocean to the hemisphere where deep water forms and that there is a significant net heat transport across the Equator from one hemisphere to the other hemisphere, as shown in Fig. A3c and f.
Because the 570 Ma continental configuration does not have any continental fragment further north than 25 • N, the channel flow is much stronger than that of the 720 Ma continental configuration (not shown) and the deep-water formation in the SH is also much stronger. Indeed, the deep-water cell for the 570 Ma continental configuration (Fig. A4c) has a strength of ∼ 60 Sverdrup (Sv) compared to ∼ 35 Sv for the 720 Ma continental configuration (Fig. A4a). This is also the reason why switching from the RIDGE method for the control of the ocean circulation to the DIF method has a larger impact for the 570 Ma continental configuration as will be described further in what follows.
Comparison of Fig. A4c and d clearly indicates that the deep-water formation process is greatly weakened due to the increased tracer diffusivities. This is expected since it is known that increased isopycnal and eddy-induced diffusivities act to weaken the effect of the channel flow on deepwater formation (e.g., Gnanadesikan et al., 2003Gnanadesikan et al., , 2004. The weakening of deep-water formation in the SH causes less heat to be transported to the SH (Fig. A3f) and thus the atmosphere there to be colder (Fig. A3d), while the effect on the NH is opposite. The compensation by the inter-hemispheric atmospheric heat transport is reduced (Fig. A3e) since the ocean heat transport is now more symmetric. The changes for the 720 Ma continental configuration are much less significant since it had weaker deep-water formation in the RIDGE method. Therefore, the primary influence of enhanced tracer diffusivities on the climate system is to weaken the deepwater formation process, and redirect ocean heat transport. The consequence is that the climate in the two hemispheres is more symmetric (red curves in Fig. A3a and d).
For the 570 Ma continental configuration, the asymmetric change in the two hemispheres almost completely compensates so that the global mean surface temperature stays almost the same as that predicted by the RIDGE method ( Fig. 3b in the main text). For the 720 Ma continental configuration, the cooling effect in the SH is similar but there is no obvious warming in the NH (Fig. A3a) so that the global mean surface temperature is ∼ 1.3 • C lower than that predicted by the RIDGE method (Fig. 3a in the main text).
A significant question concerns the issue as to whether use of the DIF method will have a significant influence on the hard snowball bifurcation point. We expect not. As the Earth becomes increasingly close to entering a hard snowball state, the climate in the two hemispheres tends to be more and more symmetric (see Fig. 1b for the atmosphere and Fig. 2 for the ocean of Yang et al., 2012b) and the deep ocean circulation becomes less and less important, so that the influence of the DIF method is therefore expected to be minimal. The fact that switching from the RIDGE method to the DIF method does not have any significant effect on either the NH surface temperature (Fig. A3a) or sea-ice front position (compare Fig. A2a and b) for the 720 Ma continental configuration may be taken as an analog for the influence of the DIF method on an already relatively symmetric climate. However, because the upper ocean circulation (Fig. A4) and ocean heat transport ( Fig. A3c and f) are slightly enhanced in the tropical region when tracer diffusivities are increased, the DIF method might have a slight effect of preventing the advance of the sea-ice fronts towards the Equator once they enter the tropical region.
As an effort to further quantify the influence of enhanced tracer diffusivities on the identification of bifurcation points, two simulations employing slightly lower values, 3200 m 2 s −1 and 2400 m 2 s −1 , were carried out for the 570 Ma continental configuration with pCO 2 = 140 ppmv. Compared to the simulation with tracer diffusivities of 4000 m 2 s −1 (blue curve in Fig. 5b), the global mean surface temperature is lower by 0.14 and 0.23 • C, respectively. Therefore, enhanced tracer diffusivities do produce warmer climate due to enhanced ocean heat transport, thus leading to probably an underestimate of the critical pCO 2 for hard snowball Earth formation, but the magnitude of this effect is small compared to other uncertainties of the model such as sea-ice albedo, sea-ice dynamics, cloud parameterization, etc.