Large equatorial seasonal cycle during Marinoan snowball Earth

See allHide authors and affiliations

Science Advances  03 Jun 2020:
Vol. 6, no. 23, eaay2471
DOI: 10.1126/sciadv.aay2471


In the equatorial regions on Earth today, the seasonal cycle of the monthly mean surface air temperature is <10°C. However, deep (>1 m) sand wedges were found near the paleoequator in the Marinoan glaciogenic deposits at ~635 million years ago, indicating a large seasonal cycle (probably >30°C). Through numerical simulations, we show that the equatorial seasonal cycle could reach >30°C at various continental locations if the oceans are completely frozen over, as would have been the case for a snowball Earth, or could reach ~20°C if the oceans are not completely frozen over, as would have been the case for a waterbelt Earth. These values are obtained at the maximum eccentricity of the Earth orbit, i.e., 0.0679, and will be approximately 10°C smaller if the present-day eccentricity is used. For these seasonal cycles, theoretical calculations show that the deep sand wedges form readily in a snowball Earth while hardly form in a waterbelt Earth.


Deep (>1 m) sand wedges have been recognized in the Marinoan [~635 million years (Ma)] glacial strata on various continents, a number of which probably formed at very low paleolatitudes (≤10°) (1). Sand wedges form when cold brittle ground is cooled rapidly so that cracking occurs due to thermal contraction (2). The rapid cooling events need to last for two or more days for their effect to penetrate deep into the ground. Continuous (slow) cooling or coldness is necessary to let the cracks widen or remain open (3). The open cracks are then filled with windblown sands under dry conditions (4). Therefore, both cold annual mean temperature and large seasonal cycle are required for deep sand wedges to form. These initial sand wedges tend to grow in subsequent years because they have lower tensile strength than the surrounding environment; new cracks will likely form within these existing wedges.

Deep sand wedges only form in high latitudes on present-day Earth. Although temperatures can be cold at high elevations of low-latitude mountains, sand wedges have not been observed to form there because of weak seasonal cycle of temperature (5). On the basis of the geographical distribution of present-day sand wedges, the amplitude of seasonal cycle (ASC; defined as the difference in monthly mean surface air temperature between the warmest and the coldest months) required to form a sand wedge has been estimated to be >30°C (5, 6). The finding of sand wedges near the paleoequator therefore means that the equatorial region during the Marinoan was dry and cold and had a large seasonal cycle.

Widespread glacial deposits near sea level at low-latitude regions occurred multiple times in the Late Neoproterozoic [0.75 to 0.54 billion years (Ga)], leading to the proposal that all or most of the continents were covered by ice sheets during these glacial events (7). In the sensu stricto snowball Earth hypothesis (8), sea ice was global and thick (~1000 m), and the climate was extremely cold even at the equator. An alternative hypothesis has also been proposed, i.e., the waterbelt Earth (often also called a slushball or soft snowball Earth), in which the sea ice was thin and did not cover the tropical oceans (911). Both hypotheses face their own challenges [e.g., (1214)], but they both give cold temperature within the equatorial region; the annual mean surface temperatures averaged over equatorial land regions are approximately −40° and 0°C for the snowball and waterbelt Earth, respectively. They may thus potentially explain the existence of sand wedges there.

A large ASC near the equator was, however, not found for either the snowball or waterbelt Earth in previous climate modeling studies (11, 15); these studies generally found small ASC (~10°C or smaller), most likely inadequate for producing sand wedges (3, 16). Maloof et al. (17) tried to bypass this problem by suggesting that the deep sand wedges could have formed on a brittle ground under diurnal cycles. Their theoretical calculations showed that this was possible, but the mechanism cannot be confirmed by observations on Earth and has, thus, been questioned (5). Ewing et al. (18) argued on the basis of observational evidence that some of the equatorial sand wedges formed in the active layer under a relatively warm climate and weak seasonal cycles. However, the maximum ASC in the equatorial region they obtained with a one-dimensional (1D) vertical heat transfer model was again very small (~8°C); it is similar to the ASC of present-day equatorial regions where no sand wedges are known to form (5).

The previous climate modeling studies focused on the initiation and deglaciation of a snowball Earth (15, 1922) or on the general climate during a snowball Earth under the present-day Earth orbit (11, 23, 24). As we will demonstrate below, their model configurations and experimental designs were not the most suitable for investigating the equatorial ASC during snowball Earth events. In this study, we examine the equatorial ASC for either a snowball or a waterbelt Earth in a general circulation model (GCM; see Materials and Methods). Our approach differs from previous modeling studies in that we explore different orbital parameters than what exists today and the effects of surface topography. We show that the ASC is spatially highly variable even at the same latitude and can have large values at many equatorial locations. Using the same method as that of Maloof et al. (17), we also demonstrate that these ASCs are large enough to form deep sand wedges.


Equatorial seasonal cycle

Our model shows that an equatorial ASC of much greater than 20° can be achieved for a snowball Earth when the orbital parameters of the Earth are varied even just within the range of those of the Pleistocene (25). The ASC of surface air temperature is largely constrained by the seasonal cycle of local solar insolation and the heat capacity of the surface. For the equatorial region (defined herein as within ±11°, or the width of three grid cells on either side of the equator in our model, where each grid cell is approximately 3.75°×°3.75 in size), the seasonal variations in solar insolation are insensitive to the variations of obliquity within the range of 22.5° to 24.5° (fig. S1, A and D) but are much more sensitive to variations in eccentricity (0.016 to 0.0679) and precession (fig. S1). Most studies used present-day orbital configurations [e.g., (20, 21, 23)], which have a small eccentricity (0.0167). Such a small value gives a maximum equatorial ASC smaller than 23°C during a snowball Earth (Fig. 1A). The maximum ASC near the equator in the favorable hemisphere (depending on the precession) increases steadily with eccentricity, approaching 30°C when the eccentricity is 0.057 (Fig. 1, A to C). Likewise, the equatorial ASC in a waterbelt Earth [see its climate in (20)] increases in general when the eccentricity is increased (Fig. 1, D to F). However, the maximum ASC over the region is still much less than 20°C even when the eccentricity is 0.057 (Fig. 1F). Therefore, the equatorial seasonal cycle during snowball Earth events could indeed be much stronger than that of the present day, encouraging the growth of sand wedges.

Fig. 1 Amplitude of the seasonal cycle of surface air temperatures in the tropical regions for different eccentricities.

(A to C) Snowball Earth where the entire ocean is covered by thick ice (white color), simulated by CAM3. (D to F) Waterbelt Earth where the tropical ocean is ice free (blue color), simulated by CCSM3. The ice sheet configuration in Fig. 2B is used in all simulations here. The eccentricity (as denoted by ecc in A to C) of the Earth increases from 0.017 at the top to 0.057 at the bottom in both columns. Obliquity = 24.5°, precession = 270° for all simulations. pCO2 = 0.1 bar and 1.4 × 10−4 bar in the left and right panels, respectively. The white patches in the left panels represent areas covered by both land ice (see Fig. 2B) and thick sea ice, while only the former exists in the right panels.

It is found that the snowball Earth climate takes at least a few hundred years to reach equilibrium (fig. S2); the pattern of snow depth takes a long time to stabilize because of the weak hydrological cycle. Running the model for long enough is important, at least for the atmospheric GCM used here, i.e., Community Atmosphere Model version 3 (CAM3); the equatorial ASC can increase by as much as 12°C when the model duration is increased from 50 to 300 years, as the snow within the tropical region is gradually sublimed away. This problem does not exist for the waterbelt Earth because its climate has to be simulated by a fully coupled atmosphere-ocean GCM. These models are always run for at least a few hundred years for the surface climate to reach equilibrium (fig. S2, B and D).

For sand wedges to form, there needs to be exposed land surfaces, i.e., not covered by ice sheets or thick snow. Some land surface may not be covered by ice during a snowball Earth. For example, entire tropical continents may not be ice covered (Fig. 2B) in the very early stage of a snowball Earth (26), or in the late sage, depending on the greenhouse gas concentrations and orbital configurations (27); the coastal regions may be ice free even if entire continents are covered by ice sheets (Fig. 2A), like the Dry Valleys of the present-day Antarctica, where ice flow is restricted by nearby topography and by katabatic winds [e.g., (28)]. In a waterbelt Earth, which is much warmer than a snowball Earth, the tropical continents should be exposed more easily.

Fig. 2 Surface topography of continental ice sheets (blue-yellow-pink color) and the exposed periglacial region (gray).

The ice sheets have developed in the tropical region in (A) while not yet in (B). The width of the exposed coastal region is one grid cell (~3.75°) in the climate model, but it appears larger here because of imprecise contour plotting of the ice sheets. The precise area of the exposed region in the climate model can be found in Figs. 3 and 5. In (C), a straight mountain range is added compared to (B). Its peak height is 2 or 4 km in different experiments (tables S1 and S2), and it is not covered by ice.

Here, we prescribe two types of continental ice sheet configurations as shown in Fig. 2; all the coastal regions are ice free in these configurations. As will be shown later, large fraction of the coastal regions will accumulate thick snow during the simulations, meaning that they should have been covered by ice sheet. However, some part of the coastal regions will remain snow free, and we examine the ASC over these regions. We assume that it is always possible that the topography around some of these snow-free coastal regions is in such a way that ice flow from inland or neighboring region is diverted away from them (one example is shown in fig. S3, in which a coupled simulation of an ice sheet model and CAM3 has been carried out; see Materials and Methods for more details). The two ice sheet configurations are tested for both the snowball and waterbelt Earth. The results regarding the ASC over ice- and snow-free regions are described in the following two sessions.

Seasonal cycle on exposed land surfaces in a snowball Earth

For the snowball Earth, a large fraction (>30%) of the initially exposed land surfaces remains exposed to air (i.e., snow is thin; Fig. 3) after an integration of at least 600 years (fig. S2). The land surface is considered exposed if the annual mean snow thickness is <5 cm, which should be thin enough so as not to insulate the ground from seasonal variations of surface temperature [e.g. (16)]. At some locations, snow remains thin throughout the simulation. At other locations, however, the snow depth may oscillate markedly during the first few hundred years of simulation but equilibrates to small values (e.g., <2 cm) at the end (fig. S4 and movie S1). This again indicates the importance of continuing the model runs for a sufficiently long duration.

Fig. 3 Amplitudes of seasonal cycles in terms of the monthly mean surface air temperature in a snowball Earth scenario.

Precession is 90° for (A) and (C), and 270° for (B) and (D). Only values over exposed land areas where the annual mean snow thickness is less than 5 cm are shown. The values are averaged over the past 100 years of the respective simulations. pCO2 = 0.1 bar, eccentricity = 0.0679, obliquity = 24.5°, precession = 270°. Each square is one grid cell in the model, whose size is approximately 3.75° × 3.75°. The numbers 1 to 6 in (D) indicate the locations for which the time series of snow depths are shown in fig. S2. The variations in ASC of soil temperature with depth and the seasonal temperature variation at 1-m depth for locations 1 to 3 are shown in (E) and (F), respectively. Note that the colorbar here is different from that in Fig. 1.

The ASC is greater than 24°C at many of these exposed locations and even greater than 30°C at a few places (yellow and orange cells in Fig. 3, see also Fig. 4) when the eccentricity is 0.0679. For the Northern Hemisphere, the strongest (weakest) seasonal cycle is obtained when the precession is at approximately 270° (90°), at which the summer occurs at the perihelion (aphelion). The opposite is true for the Southern Hemisphere (Figs. 3 and 4). The ASC has interannual variability, with a standard deviation generally between 1° and 2°C. Note that each square in Fig. 3 represents an area of approximately 400 km by 400 km, large enough to accommodate all of the Whyalla Sandstone in the Gawler Craton of Australia, where the most spectacular sand wedges were found (5).

Fig. 4 The sensitivity of the maximum amplitudes of seasonal cycle at each latitude to different parameters.

(A) Sensitivity to the atmospheric CO2 concentration, with eccentricity = 0.0679, obliquity = 24.5°, and precession = 270°. (B) Sensitivity to obliquity, pCO2 = 0.1 bar, eccentricity = 0.0679, and precession = 270°. (C) Sensitivity to precession, pCO2 = 0.1 bar, eccentricity = 0.0679, and obliquity = 24.5°. (D) Sensitivity to eccentricity, pCO2 = 0.1 bar, obliquity = 24.5°, and precession = 180°. The data points (filled circles) are based on snowball Earth simulations with continental ice sheet distributions as shown in Fig. 2B. The filled squares in (B) are obtained for the waterbelt Earth scenario (see Fig. 5D), and the orbital parameters are the same as those in the snowball Earth simulations in the panel, but pCO2 = 1.4 × 10−4 bar for the M0 and M1 cases and 2.8 × 10−4 bar for the M2 case. The columns denoted by M0, M1, and M2 are the same except that in M1 and M2, a north-south–oriented mountain of 2 and 4 km above sea level (see Fig. 2C) is added.

The ASC of soil temperature is smaller than that of the surface air temperature. It decreases rapidly in the first few centimeters and then relatively slowly with depth in a near-exponential manner. For a few points chosen in Fig. 3D, the ASCs are all greater than 20°C at depths of near 3 m (Fig. 3E) and greater than 28°C for one of them. The minimum monthly temperature is between −35° and −40°C, and the maximum is between −18° and −9°C at a depth of 1 m (Fig. 3F). This large ASC at depth makes the formation of deep sand wedges possible.

The ASC is weakly dependent on the atmospheric CO2 concentration, pCO2; as pCO2 increases, the maximum ASC within the equatorial region decreases slightly but remains large (Fig. 4A). The annual mean surface temperature is below −20°C in the equatorial regions even for a high pCO2 of 0.1 bar (fig. S5), meaning that the whole region is permafrost. However, the maximum hourly temperatures can be greater than 0°C at various locations over the year (fig. S5), which may explain the occurrence of fluvial deposits that are sometimes associated with the sand wedges (18). The maximum ASC is not sensitive to obliquity for the range considered here (22.5° to 24.5°; Fig. 4B) but is sensitive to eccentricity (0.017 to 0.0679; Fig. 4D), as expected from the seasonal cycle of local solar insolation (fig. S1). Precession determines the hemisphere in which the maximum ASC occurs (Fig. 4C). Overall, the maximum ASC can be close to or larger than 30°C over a wide range of orbital parameters and pCO2, as long as the eccentricity is greater than ~0.045.

Seasonal cycle on exposed land surfaces in a waterbelt Earth

In a waterbelt Earth, the snow thickness is small in most equatorial regions due to the warmer climate [e.g., (29)], but the ASC is in the range of 6° to 20°C even when the eccentricity is at its maximum (0.0679; Fig. 5). The maximum ASCs in the six equatorial latitudinal bands (squares in Fig. 4B) are more than 10°C smaller than those for a snowball Earth. Adding a hypothetical mountain to the middle of the continent (Fig. 2C) increases the maximum ASC slightly, but only to ~23°C (Fig. 4B; see also fig. S7). The annual mean surface air temperatures in the equatorial regions are high, in the range of −8° to 8°C, so the ground surface is seasonally frozen in most parts of the equatorial regions in a waterbelt Earth.

Fig. 5 Amplitudes of the seasonal cycle in terms of monthly mean surface air temperatures in a waterbelt Earth scenario.

The deep blue region is ocean, and the hatched area is where the grid-cell sea ice fraction is greater than 50%. The results are averaged for the past 10 years of the simulations of the fully coupled model, CCSM3. Eccentricity = 0.0679, obliquity = 24.5°, and the precession is indicated at the lower left corner of each panel. pCO2 = 1.0 × 10−3 bar in (A), 7.0 × 10−4 bar in (B), and 1.4 × 10−4 bar in (C) and (D). The climates will enter a snowball state when pCO2 is lowered further. The variations in ASC of soil temperatures with depth and the seasonal temperature variations at 1-m depth for points 7 to 9 are shown in (E) and (F), respectively.

The ASC of soil temperature in a waterbelt Earth is mostly less than 15°C at depths of 1 m and deeper (Fig. 5E). The minimum monthly temperature at 1-m depth could reach −16.5°C for point 8 in Fig. 5D but barely reached −13°C for other points. The maximum monthly temperature at 1-m depth is between −3° and 0°C. Therefore, these sites have a permafrost subsurface and an active surface layer.

The seasonal cycle of surface radiative forcing is similar in magnitude between a snowball and a waterbelt Earth (fig. S8, D and H). The relatively small ASC in a waterbelt Earth should be due to the following reasons: (i) the heat capacity of the wet soil layer in a waterbelt Earth is larger than that of the frozen soil layer in a snowball Earth (compare fig. S8, C and G), since the specific heat of water is almost twice that of ice and four times that of clay (see Materials and Methods). Note that the phase change between soil water and ice in a waterbelt Earth should also increase the effective heat capacity of the soil layer but is not counted in fig. S8G. (ii) The latent heat flux due to evaporation of water over land surface in a waterbelt Earth is nearly three times stronger than that in a snowball Earth (fig. S8, D and H); its seasonal cycle is in antiphase with the seasonal cycle of radiative forcing and cancels approximately one-third of the latter (fig. S8H). (iii) The seasonal cycle of sensible heat flux is similar to that of the latent heat flux in a waterbelt Earth and is also significantly larger than that in a snowball Earth (fig. S8, D and H). This is likely due to the more unstable atmospheric temperature profile in a waterbelt Earth than in a snowball Earth. Because of these reasons, the seasonal variation of surface radiative flux in a snowball Earth has to be balanced by a large seasonal variation of the outgoing longwave radiation (fig. S8D), namely, the surface temperature.

Crack formation

A relatively simple model as used by Maloof et al. (17) can be used to test whether deep cracks are able to form under the forcing of ASC calculated above. In this model, the occurrence of a crack is determined by a crack edge intensity factor (K). If it is larger than the fracture toughness of the material, then cracking occurs. The fracture toughness is negatively correlated with temperature and grain size of soil particle and changes nonmonotonically with the ice content in the soil. Therefore, a wide range between 0.1 and 0.8 was specified in (17). The rheology of the soil is also uncertain. Three deformation mechanisms were proposed in (17), i.e., dislocation climb-limited steady-state creep (clime limited), grain boundary sliding–accommodated basal slip–limited steady-state creep (GBS limited), and drag-limited transient creep. Only the first two are adopted here because the third one was proposed to be more appropriate for deformation under diurnal temperature variations. The intensity factor is calculated using the Eq. 12 in (17) and all the parameter values there, except that the soil temperatures herein come from climate model simulations rather than being specified as an idealized sin function.

The intensity factors for points 3 and 6 in a snowball Earth (Fig. 3D) and points 8 and 9 in a waterbelt Earth (Fig. 5D) are shown in Fig. 6. They achieve maximum values when the temperatures are approaching the minimum (red boxes in Fig. 6, A to D). Their values are quite sensitive to the creeping mechanisms assumed, and the GBS-limited creeping gives much lower intensity factor than the climb-limited creeping. For point 3, which has an ASC of 29.8°C (Fig. 3E), cracking can occur at depth >2.5 m for either of the creeping mechanisms (Fig. 6E). The intensity factors at 2.5-m depth for points 1 and 2 are only slightly smaller than that for point 3. For point 6, which has a surface ASC of 17.7°C, deep cracking is much less likely if the GBS-limited creeping dominates in the soil layer (red curves in Fig. 6F); the relatively deep snow during winter (as can be inferred from the very small diurnal cycles of temperature in Fig. 6B) prevents the soil temperature to drop rapidly.

Fig. 6 Soil cracking calculated for a few typical points.

(A to D) are times series of 2-hourly soil temperature at 20-cm depth and (E to H) are the corresponding depth profiles of maximum intensity factor. The locations of points 3 and 6 can be found in Fig. 3D and those of points 8 and 9 in Fig. 5D. The time periods at which the intensity factors are maximum are indicated by the red boxes on the left. The two dashed lines in the right panels indicate the range of fracture toughness of frozen soils. The solid curves labeled ice in the right panels are the intensity factors calculated assuming the creeping mechanisms of polycrystalline ice, while the dashed curves labeled soil1 are for frozen soils that have viscosity 10 times that for pure ice. The viscosity of ice-cemented sand could have viscosity 100 times that of pure ice (17).

For points 8 (ASC = 18.3°C) and 9 (ASC = 13.0°C) in a waterbelt Earth, the temperature below 20-cm depth is capped at 0°C, due to the buffering effect of the melting of soil ice (Fig. 6, C and D). Cracking at depth >1 m can hardly occur if the GBS-limited creeping dominates (red curves in Fig. 6, G and H). During an especially cold event, cracking may be able to occur at 1-m depth if the fracture toughness of the soil is exceptionally weak (Fig. 6, D and H). From the intensity factor calculated for these and all other points in Figs. 3D and 5D, deep cracking can occur easily in either a snowball or a waterbelt Earth if the climb-limited creep dominates (blue curves in Fig. 6, E to H). Given the fact that deep sand wedges were not found so commonly in the Neoproterozoic strata, the GBS-limited creeping was probably the dominating mechanism of soil deformation. Then, deep sand wedges could form near the equator in a snowball Earth as long as the ASC is greater than ~24°C, but could hardly form in a waterbelt Earth.


Uncertainties in the conditions for sand wedge formation

On the basis of information from a few sites in Antarctica and in the high Arctic, Williams et al. (5) inferred that an ASC of >30°C or even ~40°C is required to form deep sand wedges in a permafrost. Such inference, however, may not be extrapolated to the situation in a snowball Earth. On the present-day Earth, the permafrost that has a cold annual mean surface temperature (e.g., −20°C; this would be close to the equatorial temperature in a snowball Earth; see Fig. 6, A and B) is located in the high latitudes, where the seasonal cycle is always strong. Therefore, it cannot be known from modern observations whether sand wedges can form or not at weaker seasonal cycles in a cold permafrost. Our relatively simple theoretical calculations show that sand wedges may be able to form at ASCs smaller than 30°C and greater than 24°C. Such ASCs are entirely possible within the equatorial region during a snowball Earth, especially when the interannual variability of ASCs is considered.

Relict sand wedges have long been suggested to form on seasonally frozen ground that has much higher mean annual temperatures (4). For example, sand wedges were found in France (30) and Ordos in China (31), and some of them are very deep (>3 m) in the former location. They were formed during the Pleistocene glaciations, with estimated annual mean surface temperatures near or only a few degrees below 0°C. Not until recently have direct observations of contemporary formation of thermal contraction cracks (at 1-m depth) been made in the Northwest Territories, Canada, where the annual mean surface temperature is −4.1°C (16). These cases have greatly relaxed the climatic conditions for sand wedge formation.

Although there are no reliable observational constraints on the ASC, simulations of the climate during the coldest period of the Pleistocene [e.g., the last glacial maximum (LGM), ~21,000 years ago] by many climate models show a mild ASC (17.0° ± 3.5°C at the center of France; 1σ uncertainty) over all of France (fig. S9). For Ordos in China, the modeled ASC is indeed large (32.5° ± 3.5°C). The climate models are not yet perfect, but the consistency among the models indicates that the formation of sand wedges may require substantially smaller ASC values depending on the sedimentological and climate settings. Therefore, this again indicates that the conditions required for sand wedge formation as put forward by Williams et al. (5) may have been exaggerated. Dry environments are available in a waterbelt Earth, e.g., points 7 to 9 in Fig. 5 (fig. S10), but the equatorial ASC is only marginally large enough even if the threshold ASC is relaxed to that in France during the LGM.

Implications for the snowball Earth hypotheses

The deep sand wedges in southern Australia have a paleolatitude of 8.8° ± 3.2° (32). Ewing et al. (18) obtained a weak ASC (~8°C) near 10° latitude with a simple model of ground temperature and, thus, questioned the validity of the snowball Earth theory. Our GCM simulations show that one important reason that they obtained this small ASC is because the climate was assumed to be zonally uniform in their model. Moreover, they determined the values of some key parameters (e.g., soil heat capacity) in their formulae by fitting the results to the diurnal and seasonal temperature variations over the present-day land, which will inevitably give weak seasonal cycle. These prevented them from obtaining a larger ASC within the equatorial region.

Our results show that a large ASC (>30°C) is fully possible at equatorial latitudes in a snowball Earth, even if Earth had a normal obliquity (i.e., between 22.5° and 24.5°). Therefore, the existence of equatorial sand wedges poses no danger to the snowball Earth scenario. The formation of equatorial sand wedges is much less likely in a waterbelt Earth according to our theoretical calculations (Fig. 6, G and H), but the possibility may not be entirely ruled out before a complete explanation for the relict sand wedges in France is available. The lower the paleolatitudes of sand wedges, the less likely that Marinoan glaciation was a waterbelt Earth (Fig. 4B). Therefore, further examination of the paleolatitudes of the reported Marinoan frost-shattered carbonates or sand wedges is warranted (27, 33).


Continental distribution

The continental distribution used in this study is that reconstructed by Li et al. (34) for 720 Ma ago. Because the reconstructed continental distribution for 650 Ma is not so different from that for 720 Ma, the latter is taken to be appropriate for both the Sturtian and Marinoan glaciations. There are still large differences in the details between the continental distributions reconstructed by different groups [e.g., (35, 36)]. Therefore, we have focused on discussing the general features of seasonal cycles in the equatorial region during a snowball Earth, rather than the seasonal cycles on a specific continent.

Continental ice sheets on a snowball Earth

The most comprehensive way of examining the exposure of land surface would be to simulate the ice sheet dynamics and climate at the same time by using a coupled GCM and ice sheet model. However, this will be extremely time-consuming, not only because it will take a very long time [a few hundred thousand of years; e.g., (26, 37)] for the model to reach equilibrium but also a large amount of work has to be done manually to design surface topography in some local region. The latter part is needed to prevent the region from becoming ice covered due to ice flow, even though it is not able to grow ice locally because of either high temperature or too little precipitation. This particular region cannot be identified precisely in advance, since we cannot predict which (small) region will have the largest ASC. Therefore, we took two different configurations of ice sheets that were simulated by an ice sheet model coupled to an energy balance model (EBM) and prescribe them as boundary conditions in the GCM. The prescribed ice sheets may affect the detailed distribution of ASCs over exposed land surface but should not affect the overall characteristics of how they change with eccentricity and other forcings.

The continental ice sheets specified in the GCMs are obtained by a 3D thermomechanical dynamic ice sheet model, University of Toronto (UofT) GSM, coupled with a 2D global EBM (10). The shallow ice dynamical equation is solved using the finite difference method in the UofT GSM, with a horizontal resolution of 50 km and 35 layers in the vertical direction. The EBM is solved using the spectral method with a resolution of triangular truncation of the spherical harmonic functions at degree 11. The time step for the ice sheet model is generally 25 years and is coupled to the EBM once every 500 years. The model is run for 1.2 million years. The solar constant is set to be at 94% of the present-day value, and pCO2 is set to 7.0 × 10−5 bar.

Starting from an ice-free state, the coupled model is integrated forward for more than 500 thousand years. The ice sheet thickness at two different stages of the ice sheet evolution is selected from the time series: One is at the equilibrium state (Fig. 2A; referred to as FULL hereafter), and the other is at an intermediate stage where the continental ice has expanded to near 20° latitude (Fig. 2B; referred to as PARTIAL hereafter). Slight modifications are made to the ice sheet thickness distribution so that the coastal regions (a width of one grid cell) have no ice. The area of the coastal region in Fig. 2A is 8.9% of the total area of the continents. Such a scenario would have a negligible effect on the global climate because of the small area but will facilitate identification of all possible regions for potential sand wedge formation in a single model run.

In one case (fig. S3), the ice sheet model is coupled to a GCM CAM3 (see below) to obtain the ice sheets in a snowball Earth. In this simulation, a few narrow (width, 50 km) mountain ranges are designed near the coast (fig. S3A) to prevent ice flow from covering certain coastal regions (fig. S3B), as would happen in the real world. The CAM3 is first run for 150 years with the land surface of mid to high latitudes (>27°) specified as glacier. The climate fields averaged over the past 10 years are then used to drive the ice sheet model, UofT GSM, for 1 million years. The ice sheet distribution obtained by the ice sheet model is input into the climate model to obtain the new climate fields. The whole cycle is repeated for three times so that the coupled system reaches quasi-equilibrium. The ice sheet distribution and ASC are, thus, obtained simultaneously in this simulation. However, the coupling involves large amount of human intervention, because we do not know where to put the mountain ranges beforehand and, thus, have to do it by way of trial and error. Moreover, the locations will be different for different orbital configurations and CO2 concentrations. Therefore, it is not practical to do such simulation for all the experiments tested in this study. The results of this simulation do show that it is possible for coastal regions to remain non–ice covered and have large equatorial ASC. Therefore, it should give us enough confidence that the results from the simplified approach, i.e., the ice sheets being prescribed and all the coastal regions being exposed, are trustworthy.

Waterbelt Earth simulations using CCSM3

The National Center for Atmospheric Research (NCAR) Community Climate System Model version 3 (CCSM3) was used here to estimate the seasonal cycle of surface air temperatures in a waterbelt Earth. The model simulates the dynamics and thermodynamics of four components of the climate system, namely, atmosphere, ocean, land surface processes, and sea ice. The respective model components are the CAM3, the Parallel Ocean Program (POP; based on version 1.4.3), the Community Land Model version 3 (CLM3), and the Community Sea Ice Model version 5 (CSIM5). These four components were linked through a flux coupler (CPL6), with no flux corrections between the components (38).

Both CAM3 and POP are 3D GCMs and are run at resolutions of T31 (approximately 3.75° × 3.75°) and “gx3” (100 and 116 grid points in the zonal and meridional directions, respectively) in the horizontal and 26 levels and 25 levels in the vertical direction, respectively. The top of the atmosphere is 3 hPa, and the depth of the ocean is 4 km. CLM3 is run on the same grid as CAM3. The surface albedo of land, if not ice covered, is 0.09 to 0.18 in the visible band and is 0.18 to 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). If the land is covered by land ice, then the albedo is 0.8 in the visible band and 0.55 in the near-infrared band. The albedo of new snow is higher, 0.95 and 0.65 in the visible and near-infrared bands, respectively, but can decrease to minimum values of 0.76 and 0.325 as the snow ages (39). The soil is assumed to be an average soil (i.e., loam), which contains 15% clay, 43% sand, and 43% silt. CSIM5 includes both the thermodynamics and dynamics of the sea ice and is run on the same horizontal grid as that of the ocean. The time steps for the atmosphere, ocean, land, and sea ice components are 30, 120, 30, and 60 min, respectively. Fluxes are exchanged every 1 hour between the atmosphere, land, and sea ice components and every day between the ocean and other components. These setups are the same as those in Liu et al. (20, 29), where the initiation of a snowball Earth is studied.

The solar constant is fixed at 94% of the present-day value, 1285 W m−2, and the atmospheric concentrations of CH4, N2O, and chlorofluorocarbons (CFCs) are fixed at 0.8056, 0.2767, and 0 parts per million by volume (ppmv), respectively. The orbital parameters, e.g., the eccentricity, obliquity, and precession, can be easily modified in CAM3 or CCSM3. The precession is defined as the longitude of the vernal equinox measured counterclockwise from the perihelion when the orbit is viewed from above. The initial states of the ocean and sea ice are from the simulations in Liu et al. (20) for pCO2 = 7.0 × 10−4 bar. The atmosphere is initialized from a static state with vertically and zonally uniform temperatures, which vary from 15°C at the equator to −25°C at the poles. The coldest climate for a waterbelt Earth is obtained for each of the two ice sheet distributions in Fig. 2, when pCO2 levels are approximately 7.0 × 10−4 and 1.4 × 10−4 bar, respectively. For some orbital configurations, Earth enters a snowball Earth stage more easily than for others, and higher pCO2 (up to 1.4 × 10−3 bar) is used for these configurations (table S1).

The model is integrated for 600 years for the surface climate to reach quasi-equilibrium (fig. S2). Unlike the simulations for a snowball Earth (see below), there is almost no permanent snow accumulation in the equatorial region for a waterbelt Earth, so it is not necessary to integrate the model for very long. The final 10 years are averaged to calculate the seasonal cycle of temperature in Figs. 4 and 5 in the main text. The amplitude of the seasonal cycle will be even smaller if more years of data are averaged. A total of 17 simulations are carried out for different orbital, ice sheet, and mountain configurations (table S1).

Snowball Earth simulations using CAM3

Because there is no open ocean, and the sea ice is very thick and is almost static for a snowball earth, it is sufficient to model such a climate with only the atmospheric (CAM3) and land components (CLM3) of the CCSM3 model. The resolutions for the two components are the same as those used in the CCSM3 simulations above. The solar constant is also fixed to be 94% of the present-day value. Because the sea ice is thick (~1000 m) and moves very slowly, it can be considered static and treated as flat land glaciers in the model, while in the waterbelt Earth simulations above, the sea-ice dynamics is explicitly modeled. In addition to the two continental ice sheet configurations considered in the CCSM3 simulations, a third configuration (referred to as FLAT hereafter) is considered here, in which the entire globe is covered by glaciers, but the surface altitude is uniformly 0. This configuration is used to demonstrate why large equatorial seasonal cycles were not obtained in previous studies (21, 23, 40).

The atmosphere is initialized from a static state with vertically and zonally uniform temperatures, which vary from 15°C at the equator to −25°C at the poles. Because the snow depth on exposed land surfaces accumulates or melts very slowly in a snowball Earth, the CAM3 simulations are carried out for at least 600 years to make sure that the snow thickness over exposed land surfaces reaches an equilibrium state. For very low pCO2, e.g., 1.0 × 10−4 bar, the simulation may be carried out for as long as 1500 years (fig. S2). The final 100 years of data are averaged to calculate the seasonal cycle of temperature as well as other climate variables, except that in Fig. 1 (E to H), where 10 years of data are used.

A snowball Earth can exist for a large range of pCO2, from ~1.0 × 10−4 to 0.1 bar in the NCAR models (19, 29). Four values of pCO2, 1.0 × 10−4, 1.0 × 10−3, 1.0 × 10−2, and 0.1 bar, are tested here. Because seasonal cycles of temperature are not highly sensitive to pCO2 (Fig. 4A), it is set to 0.1 bar in most simulations, since the climate as well as snow depth reach equilibrium sooner with higher pCO2 (fig. S2). The maximum ASC that can be obtained within the equatorial region decreases with increasing pCO2 (Fig. 4A); it changes from 37°C when pCO2 is 1.0 × 10−4 bar to 32°C when pCO2 is 0.1 bar. Therefore, the ASCs presented in the paper are probably lower-bound values. A total of 60 simulations are carried out for the different pCO2 values and orbital and ice sheet configurations (tables S2).

The heat capacity of the soil layer

The soil layer in the land module, CLM3, has a thickness of 3.34 m and is divided into 10 levels. The thickness of each level increases exponentially from top to bottom (39). The heat capacity of the soil layer is the sum of three partshc=hcsoil+hcwater+hcicewherehcsoil=hcsand*percsand+hcclay*percclaypercsand+percclayis the heat capacity of the soil itself [see Eq. 6.68 of Oleson et al. (39)]. percsand and percclay are the percentages of sand (43%) and clay in the soil (15%). hcsand, hcclay, hcwater, and hcice are the heat capacities of sand, clay, soil water, and ice, respectively. The specific heat capacities of sand, clay, soil water, and ice are 0.788, 0.883, 4.188, and 2.117 kJ kg−1 K−1 in the model, respectively. Given these values, the heat capacity of all three constituents can be calculated for each layer, since the amount of soil water (unit: kg m−2) and soil ice (unit: kg m−2) are output by the model. The equivalent heat capacity due to the phase change of soil ice is difficult to estimate. The monthly mean heat flux due to the phase change between ice and water (calculated with soil water and soil ice output every 30 min) is only up to ~3 W m−2 in the summer months, which is small compared with other energy fluxes (fig. S8, D and H). It is therefore omitted here.

Simulation of soil cracking

Soil cracking is simulated using the same method as in Maloof et al. (17), in which the soil was assumed to be viscoelastic. The dislocation climb-limited steady-state creep (climb limited) and grain boundary sliding–accommodated basal slip–limited steady-state creep (GBS limited) are the two major creeping mechanisms for the polycrystalline ice. The effective viscosity associated with both mechanisms is temperature and stress dependent. The GBS-limited creep is also dependent on the grain size, which is assumed to be 0.2 mm herein. All rheological parameters come from Table 1 of (17). They also proposed a third mechanism at diurnal time scale, i.e., the drag-limited transient creep. Since we are concerned with cracking under seasonal cycles, it is ignored in this study. Moreover, they considered two viscosities of the ice-cemented soil, which are 10 and 100 times larger than that of the polycrystalline ice, respectively, since mixed soil and ice is much stiffer than pure ice. We show results only for the viscosity of pure ice and 10 times of this value as a conservative estimate. The Young’s modulus and Poisson’s ratio of the soil are assumed to be 1010 Pa and 0.3, respectively. The coefficient of thermal expansion is set to be 2.3 × 10−5 K−1. These are also the same as those in (17).

The stress in the soil due to temperature variations is calculated using Eq. 10 of (17) and integrated numerically with the fourth-order Runge-Kutta method with a time step of 15 min. Then, the crack edge intensity factor (K) is calculated using Eq. 12 of (17). The only difference between our calculation and that in (17) is that we use the temperature output from either the CAM3 or CCSM3 model. The instant temperature is output at every time step (30 min) and each soil layer for 10 years. It is found that the peak-to-peak amplitude of the diurnal cycle in a snowball Earth is ~20°C and is no larger than that in a waterbelt Earth. This amplitude is half that assumed in (17). This smaller amplitude will reduce the intensity factor by a factor of approximately 2 and make the cracking by diurnal cycles more difficult. But, it does not overturn the major conclusion of (17) because of the large uncertainty in the viscosity and the fracture toughness of ice-cemented soil.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank P. Hoffman for the helpful discussions when preparing the manuscript and three anonymous reviewers whose comments helped improve the manuscript tremendously. The computations were supported by the High-performance Computing Platform of Peking University. Funding: Y.L. is supported by the National Natural Science Foundation of China under grant 41875090. Y.L., J.Y., and Y.H. are supported by the Chinese-Israel grant 41761144072. Y.H. is also supported by the National Natural Science Foundation of China under grant 41888101. Author contributions: Y.L. proposed the project. Y.L. and J.Y. performed modeling analyses. H.B., B.S., and Y.H. participated in designing the project and developing the manuscript. All authors contributed to the discussions and manuscript revision. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the corresponding author Y. L. according to the material transfer agreement.

Stay Connected to Science Advances

Navigate This Article