The sea-level fingerprint of a Snowball Earth deglaciation

Earth and Planetary Science Letters 399 (2014) 74–85

Contents lists available at ScienceDirect

Earth and Planetary Science Letters www.elsevier.com/locate/epsl

The sea-level fingerprint of a Snowball Earth deglaciation Jessica R. Creveling a,∗ , Jerry X. Mitrovica b a b

Geological and Planetary Sciences, California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA Earth and Planetary Sciences, Harvard University, 20 Oxford Street, Cambridge, MA 02138, USA

a r t i c l e

i n f o

Article history: Received 28 December 2013 Received in revised form 18 April 2014 Accepted 21 April 2014 Available online 28 May 2014 Editor: Y. Ricard Keywords: Snowball sea level cap carbonate fingerprint deglaciation

a b s t r a c t Cap dolostones are thought to represent deposition from seas transgressing over formerly glaciated continental margins during Marinoan Snowball deglaciation. Nevertheless, facies associations within some cap dolostones indicate that an episode of regional regression punctuated these transgressive sequence tracts. To date, inferences of sea-level change during and after the Marinoan Snowball deglaciation have been interpreted using simple, qualitative arguments. In the present study, we explore the full spatio-temporal variability of sea-level change during Snowball deglaciation and its aftermath using a gravitationally self-consistent theory that accounts for the deformational, gravitational and rotational perturbations to sea level on a viscoelastic Earth model. The theory is applied to model Marinoan Snowball deglaciation on a generalized Ediacaran paleogeography with a synthetic continental ice-sheet distribution. We find that sea-level change following a synchronous, rapid (2 kyr) collapse of Snowball ice cover will exhibit significant geographic variability, including site-specific histories that are characterized by syn-deglacial sea-level fall or stillstand. Moreover, some sites that experience syn-deglacial transgression will continue to experience transgression in the post-deglacial phase. Taken together, these results suggest that sea-level change recorded by strata capping Snowball glaciogenic units may reflect a more complicated trajectory than previously thought, including deposition that was not limited to the deglaciation phase. These simulations, as well as others that explore the response to asynchronous melting and deglaciation phases of longer duration, demonstrate that an episode of regional regression interrupting a cap dolostone transgressive sequence tract may reflect one of several processes (or their combination): (1) near field adjustment associated with rapid local melting during an otherwise global hiatus in deglaciation; (2) post-glacial uplift of sites during a period of slowing deglaciation, and (3) a transition, at some sites, from a sea-level fall dominated by post-glacial uplift to a phase of sea-level rise due to eustasy and peripheral bulge subsidence throughout an extended (order 50 kyr or greater) Snowball deglaciation. © 2014 Elsevier B.V. All rights reserved.

1. Introduction Glaciogenic diamictites interpose low-latitude marine strata on nearly every Cryogenian–Ediacaran paleocontinent (Evans, 2000; Hoffman and Li, 2009; Hoffman, 2010; Li et al., 2013). ‘Cap carbonates’ sharply overlie most late Neoproterozoic Marinoan Snowball glacial deposits (Hoffman et al., 1998). Traditionally, a cap carbonate encompasses a basal dolostone unit and an overlying limestone unit (Hoffman et al., 1998); however, cap carbonates may also include siliciclastics within stratigraphic successions recording postglacial sea-level change (e.g., the Brachina Formation, Australia, to name just one example; Rose and Maloof, 2010). Stratigraphic relationships indicate that cap dolostones represent the transgressive systems tract of a broader cap depositional sequence (Hoffman and

*

Corresponding author. Tel.: +1 (626) 395 1783. E-mail address: [email protected] (J.R. Creveling).

http://dx.doi.org/10.1016/j.epsl.2014.04.029 0012-821X/© 2014 Elsevier B.V. All rights reserved.

Schrag, 2002). Further, the consistent vertical succession of sedimentary structures and composite δ13 Ccarb chemostratigraphic data across reconstructed slope-to-platform paleoenvironments suggest a time-transgressive (diachronous) model for cap dolostone deposition (Hoffman et al., 2007; Rose and Maloof, 2010). Under the assumption that sea-level transgression requires contemporaneous deglacial melting, cap dolostones have been interpreted to record deposition over a time-scale confined to the Snowball deglaciation (Hoffman et al., 2007). Theoretical (e.g., Hoffman et al., 1998) and climate model-based (Hyde et al., 2000) predictions of Snowball deglaciation posit a melt timescale of 2–10 kyr; in contrast, paleomagnetic polarity reversals preserved within some cap dolostones imply deposition that lasted >100 kyr (Trindade et al., 2003; Kilner et al., 2005; Hoffman et al., 2007). By comparison to the bathymetric profile of modern carbonate platforms, Hoffman (2010) estimated an ∼1–1.5 km Marinoan glacioeustatic rise on the Namibian margin.

J.R. Creveling, J.X. Mitrovica / Earth and Planetary Science Letters 399 (2014) 74–85

Despite broad consensus that cap dolostones record post-glacial transgression, episodes of an early sea-level fall interrupting this transgression have been inferred from the stratigraphic transition from below wave-base limestone turbidites to low-angle, swaley and wave-ripple cross-stratified peloidal dolograinstone in cap dolostones on the Congo and Kalahari cratons (and, perhaps, Australia based on a reinterpretation of sedimentological descriptions of Kennedy (1996); Hoffman et al., 2007; Hoffman and Macdonald, 2010). In contrast (or, possibly, in addition) to sea-level fall in the cap dolostone of Australia, Rose et al. (2013) concluded that regression punctuated the transgressive syn-deglacial siliciclastic lithofacies of the underlying Elatina Formation, Flinders Range, Australia. These studies attributed regression to the loss of gravitational attraction of the sea surface to local, waning ice sheets, although Rose et al. (2013) acknowledged the possibility of, and discussed the temporal implications for, isostatic rebound in contributing to a local sea-level fall. Moreover, Hoffman and Macdonald (2010) argued that a rapid and early sea-level fall would have decreased lithostatic pressure, contributing to pore-fluid over-pressurization and the subsequent formation of bed-parallel sheet-cracks filled with isopachous cements. Sheet-crack cements have been identified in the basal-meters of cap dolostones on multiple cratons (Hoffman, 2011), indicating, based on the model of Hoffman and Macdonald (2010), that regionalized early melt (and localized sealevel fall) preceded the eustatic transgression at each of the localities where this sedimentary structure appears. The logical corollary of these assumptions is that Marinoan ice-sheets vanished asynchronously, not in unison (Hoffman and Macdonald, 2010). However, many well-studied cap dolostone successions do not exhibit sedimentary evidence for sea-level fall at the base of cap dolostones, including some that host sheet-crack cements. A number of questions arise from the above stratigraphic studies that have relevance for the interpretation of the geological record of Marinoan Snowball deglaciation. What is the plausible range of geographic variability in regional sea-level change driven by the deglaciation? Is this geographic variability a strong function of the duration of the deglaciation? Can a local geological inference of the magnitude of transgression provide a robust estimate of the globally averaged (eustatic) sea-level rise associated with the deglaciation? Finally, what specific circumstances could lead to a pronounced regional regression prior to, and perhaps also during, glacioeustatic transgression? In this study, we explore the spatio-temporal variability of sea-level change driven by Marinoan Snowball deglaciation using a gravitationally self-consistent theory and numerical algorithm that accounts for the deformational, gravitational and rotational perturbations to sea level on a viscoelastic Earth model and time-dependent shoreline migration (Mitrovica and Milne, 2003; Kendall et al., 2005). Our numerical model of the Marinoan Snowball deglaciation is configured with an Ediacaran paleogeography and a synthetic continental ice-sheet distribution. Using the model, we explore the sensitivity of the predictions to variations in both the relative synchronicity of regional ice melting and the duration of the global deglaciation phase. In the discussion below, we frame our predictions in terms of the physics of sea-level change at ‘nearfield’ versus ‘far-field’ localities. Sites in the near field of a specific region of ancient ice cover are located within ∼1000–2000 km of the margin of the ice sheet, and far-field sites are located beyond this zone. (Within a pan-continental Snowball glaciation, sites in the far-field of all ice complexes would be located in oceans; nevertheless, this designation is useful when examining both broadscale sea-level trends and melt-water volume balance.) Furthermore, in discussing time-evolving sea-level trends, we apply the terms ‘syn-deglacial’ when referencing time during a melting event and ‘post-deglacial’ to signify the time after the complete melt of global ice sheets; together, these define the ‘post-glacial’ interval.

75

Liu and Peltier (2013) describe a preliminary study of sea-level change associated with Snowball glaciation. Their analysis focused on the net sea-level fall across a multi-million year glaciation, i.e., after isostatic equilibrium is achieved in the glaciated state, and, as a consequence, the spatial variability they predict is muted (see their Fig. 12). In addition, they adopted a model of paleogeography at 570 Ma for the Marinoan (∼635 Ma) Snowball event and assumed that the position of shorelines was fixed in time (i.e., coastlines were characterized by vertical cliffs). Notwithstanding these approximations, their analysis provides a (slow) glaciation-phase complement to the detailed spatio-temporal predictions of synand post-deglacial sea-level change described herein. The timedependent crustal deformations, perturbations to the gravity field and rotational state, and eustatic sea-level variations that drive syn- and post-deglacial sea-level changes dictate the sequence stratigraphic architecture of the abundant and widely distributed cap successions and thus our numerical predictions provide an important framework for interpreting the geological record. 2. Modeling the Marinoan Snowball Earth deglaciation The configuration of Ediacaran paleocontinents remains uncertain (compare, for instance, Torsvik, 2003; Meert and Torsvik, 2004; Hoffman and Li, 2009; and Li et al., 2013). Our simulations of sea-level change in response to Snowball deglaciation were guided by the 635 Ma paleogeographic reconstruction of Li et al. (2013), but for model simplicity we assume 12 paleocontinents that, in some cases, represent an amalgamation of multiple paleocratons (Fig. 1A). Each shaded region in Fig. 1A represents the areal extent of a continent. (Adopting a paleogeography in which the continents were fully assembled into a supercontinent shows relatively small differences in the calculations described below.) Global topography prior to the Snowball glaciation is prescribed as follows: the elevation of continental interiors is set to 850 m and decreases linearly to 0 m within 350 km of the shoreline; offshore, the bathymetry drops to −150 m over a distance of 80 km (the modeled continental shelf), to −2000 m over the next 30 km (the continental slope) and to −3800 m in the next 300 km (the continental rise) where it is maintained as the depth of the abyssal plain. These values are consistent with mean values in modern topography. We model a pan-glacial scenario (Hoffman, 2009) in which ice covers all continents. Our ice model at maximum glaciation assumes an equilibrium (parabolic) thickness profile with maximum elevation proportional to the minimum width of the continent as determined from lines passing through the centroid (Fig. 1B). We adopt a proportionality constant between ice elevation and continental width such that the equivalent globally averaged (eustatic) sea-level (ESL) change associated with the maximum ice cover is 1 km, as estimated by Hoffman (2011). We assume that a single ice dome covers the following pairs of continents: B and C, E and F, G and H, and J and K (Fig. 1B). Guided by arguments in Hoffman (2011), the grounding line of each ice sheet extends to the edge of the continental shelf for the smallest ice sheets (as on continents A and L in Fig. 1A) and to 12 km down the continental slope (where local water depths are ∼1000 m) for the largest ice sheets (as on continents E/F and J/K). An equilibrium ice thickness profile is maintained during all modeled deglaciation events. In Fig. 1A, we show the approximate paleo-location of 20 cap dolostones (Li et al., 2013; Hoffman and Li, 2009; Table 1); while not exhaustive, these sites have a rich geological literature documenting the Marinoan Snowball deglaciation and are globally distributed, thus providing reference points for further discussion of sea-level change on a suite of paleo-continental margins. The paleogeography and site locations (Fig. 1A), as well as the maximum ice cover (Fig. 1B), pattern and duration of deglaciation, are subject

76

J.R. Creveling, J.X. Mitrovica / Earth and Planetary Science Letters 399 (2014) 74–85

Fig. 1. (A) Paleogeography based on the 635 Ma reconstruction of Li et al. (2013), where labels refer to the following paleocontinents: A – South China (SC); B – Tarim (plus the Central Asian Fold Belt) and northern Australia (T+NA); C – southern Australia and east Antarctica (SA+EA); D – India (I); E – the Sahara (including the Iranian Block and the Arabian–Nubian Shield), Congo, and Kalahari cratons (S+C+K); F – São Francisco and Rio Plata (Sf+R); G – West Africa and west Avalonia (WA+Aw); H – Amazonia and east Avalonia (A+Ae); I – Baltica (B); J – Laurentia, Greenland, and east Svalbard (L+G+ES); K – Siberia (Si); L – North China (NC). The sites labeled 1–20 represent the approximate paleo-location of 20 cap dolostones (see Table 1) associated with Marinoan Snowball deglaciation that have been described in the literature. (B) Grounded ice thickness (in meters) at maximum glaciation in the model Snowball Earth scenario.

Table 1 The paleo-locations of 20 illustrative Marinoan diamictites and cap carbonates (adapted from the Supplementary File of Li et al., 2013) used as reference sites for modeled deglacial sea-level predictions. The locations of the sites within the adopted palegeography are shown in Fig. 1A. Site

Region

Diamictite (Fm)

Cap carbonate (Fm)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

South China (Yangtze) Mongolia (E. Altaids) Australia (East Kimberly) Australia Australia Australia (King Island) Lesser Himalaya Oman (Jebel Akhdar) Namibia (Gariep Belt) N. Paraguay Belt Bambuí autochthon N. Damara Belt Volta Basin Adrar de Mauritanie Death Valley Mackenzie Mountains Ogilvie Mountains East Svalbard East Greenland Finnmark

Nantuo Khongoryn Moonlight Valley Olympic Elatina Cottons Breccia Blaini Fiq Namaskluft Puga Jequitai/upper S. Lagoas Ghaub Kodjari Jbéliat/Bthaat Ergil (triad) Wildrose Ice Brook (Stelfox) Upper Tindir unit 3b Wilsonbreen Storeelv Smalfjord

Lower Doushantuo (member 1) Tsagaan Oloom (OI member) Moonlight Valley (Cap Dolomicrite Unit) Mount Doreen Nuccaleena Cumberland Creek Dolostone Upper Blaini Hadash Holgat (Dreigratberg member) Mirrasol d’Oeste Sete Lagoas Keilberg Mid Sud-Banboli Amogjar Noonday Ravensthroat and Hayhook Hard Luck Lower Dracoisen Lower Canyon Lower Nyborg

to large uncertainties. Accordingly, we emphasize that our goal in the predictions below is not to reproduce the sea-level history inferred from the study of cap successions at each of these 20 sites,

but rather to use the model geometries (Fig. 1) to explore the range of possible sea-level histories driven by a Snowball deglaciation.

J.R. Creveling, J.X. Mitrovica / Earth and Planetary Science Letters 399 (2014) 74–85

77

Fig. 2. Predicted post-glacial sea-level change (in meters) computed for a globally synchronous deglaciation of the Marinoan Snowball ice cover. Ice volumes are assumed to decrease linearly in time from maximum thickness (shown in Fig. 1B) to zero over a deglaciation phase of duration 2 kyr while maintaining an equilibrium ice thickness profile. (A, B) Computed net sea-level change over the 2 kyr deglaciation phase. The latter frame is masked using contemporaneous continental paleogeography from Fig. 1A. (C, D) As in (A, B), except for the predicted sea-level change over the 10 kyr period following the end of the model deglaciation. The blue dots in frames B and D are the 20 sites shown in Fig. 1A, while the red dot in frame B is the site considered in the inset to Fig. 3A. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Our predictions are based on a gravitationally self-consistent treatment of post-glacial sea-level change developed to compute the sea-level response to the Plio–Pleistocene ice age cycles (Mitrovica and Milne, 2003; Kendall et al., 2005). The model accounts for time-varying shorelines associated with local changes in sea level and the extent of grounded, marine-based ice, as well as the impact on sea level of perturbations in Earth’s rotation. Our calculations employ the pseudo-spectral algorithm discussed by Kendall et al. (2005) with a truncation at spherical harmonic degree 256. These calculations incorporate viscoelastic deformation of a spherically-symmetric Earth in which the elastic and density structure are given by the seismic model PREM (Dziewonski and Anderson, 1981). We adopt a radial viscosity structure characterized by a 70-km-thick elastic lithosphere, upper mantle viscosity of 5 × 1020 Pa s, and lower mantle viscosity of 5 × 1021 Pa s. This model is at the lower bound on mantle viscosity as inferred from observations related to Holocene sealevel changes and mantle convective flow (Lambeck et al., 1998; Mitrovica and Forte, 2004). In the results below we include simulations that explore the sensitivity of the sea-level predictions to an increase in the adopted lower mantle viscosity. 3. Sea-level fingerprints of the Marinoan Snowball deglaciation Theoretical predictions and numerical climate modeling suggest that the collapse of Snowball ice cover occurred over a period as short as a few thousand years (Hoffman et al., 1998; Hyde et al., 2000). Furthermore, inferences of the rate of deposition of constituent sedimentary structures – including largeamplitude, oscillatory wave ripples (Jerolmack and Morhig, 2005; Allen and Hoffman, 2005); vertically aggrading ‘tubestone’ stromatolites (Hoffman et al., 2007); and aragonite crystal fans (Hoffman and Schrag, 2002) – suggest rapid accumulation of cap dolostone lithofacies. (However, Bosak et al. (2013) argued that carbonate saturation states and, therefore, depositional rates were comparable to the modern, and inferred a deglaciation time-scale of >20 kyr.) Accordingly, in our first simulation, we assume a glob-

ally synchronous melt event in which ice volumes on all continents decrease linearly to zero over a 2 kyr period (Fig. 2). The predicted sea-level change at the end of the 2 kyr ice sheet collapse exhibits significant geographic variability (Fig. 2A). In the far field of the ice cover, the computed sea-level rise peaks at 1175 m, ∼20% higher than the ESL change associated with the melt event (1 km). Near field regions once covered by ice are characterized by a predicted sea-level rise significantly lower than the eustatic value, and, near the center of the deglaciation over continents E and F (which encompass the Sahara, Congo, Kalahari, São Francisco and Rio Plata cratons), a predicted sea-level fall of ∼75 m (Fig. 2A). Given the short duration of the adopted deglaciation phase, the predicted syn-deglacial deformation of the Earth model is dominated by elastic effects (i.e., the calculation is insensitive to the adopted viscosity profile). The physics of the sea-level change in this case is straightforward. The disappearance of ice cover removes the gravitational attraction of the ice on the ocean (e.g., Clark and Lingle, 1977; Mitrovica et al., 2001), and water migrates from the near field to the far field, leading to a sea-level change in the latter region that is greater than the eustatic value. In the near field, the migration of water and post-glacial uplift of the unloaded crust both contribute to a sea-level fall that opposes the eustatic rise. How would this range of modeled post-Snowball sea-level variation appear at specific sites along Ediacaran coastlines? To illustrate this, Fig. 2B shows the instantaneous sea-level prediction of Fig. 2A, but with a mask of the adopted paleocontinental configuration (from Fig. 1A), while Fig. 3A shows time-dependent sea-level histories at a set of seven representative sites. Any site close to a continental margin adjacent to an open ocean would experience a significant (∼500 m) transgression during the deglaciation phase (Fig. 2B; see the results from −2 to 0 kyr for sites 7, 9, 10 and 17 in Fig. 3A). Alternatively, sites that were previously under significant ice cover within, for example, continental seaways or rifted margins are subject to significantly greater near field gravitational and deformational effects (Fig. 2B). In the case of site 3, this leads to a more muted sea-level rise (∼200 m) throughout

78

J.R. Creveling, J.X. Mitrovica / Earth and Planetary Science Letters 399 (2014) 74–85

Fig. 3. (A) Sea-level histories for 7 representative sites (as labeled; see Fig. 1A) computed for the case of a globally synchronous deglaciation of duration 2 kyr (as in Fig. 2). The zero on the abscissa scale refers to the end of the deglaciation phase, while the sea-level heights are plotted relative to the sea level obtained at t = 100 kyr (by which time the “background” continental hypsometry described in the text is re-established). Dotted segments of the curves denote times when the site was covered by grounded ice. Computed sea-level curves for the remaining 13 sites shown in Fig. 1A are similar in form to the following representative cases: site 3 – sites 4, 5 and 13; site 7 – site 1; site 17 – sites 2, 15, 16 and 18; site 9 – sites 6 and 14; and site 10 – sites 8, 19 and 20. The inset to the figure is the sea-level history predicted for the site marked by the red dot in Fig. 2B over a 4 kyr period beginning at the start of the deglaciation. (B) Solid line – sea level curve for site 9 computed for the globally synchronous deglaciation of duration 2 kyr considered in Fig. 2 (reproduced from frame A). Dotted line – analogous to the solid, except for the case that the volume of ice cover over the northeast sector of the ice sheet covering continents E and F in Fig. 1B decreases linearly in time, from 100–0% over 2 kyr, while the volume of the remaining global ice cover decreases linearly in the second half of this deglaciation phase (this is the scenario considered in Figs. 4A, B).

the deglaciation phase, while site 12, located under the largest of the model ice sheets, is predicted to experience a sea-level fall of ∼100 m (Fig. 3A). Site 11 is also located under this ice sheet, but nearer to its margin, at a location where near field effects essentially balance the ESL rise in the simulation (Figs. 2B, 3A). The net effect of these signals at this site is an approximately constant sea level through the deglaciation phase (Fig. 3B). Sea level is defined globally, but to observe changes in sea level requires that a site be connected to the ocean and free of grounded ice. The dotted segments of the curves in Fig. 3A indicate the time interval over which a particular site is ice covered. In the case of site 12, for example, this interval covers the entire deglaciation phase (since the site is near the center of the ice sheet covering continents E and F), and, therefore, indicates that the local fall in sea level predicted across this interval would not be reflected in the geological record. As they become ice free, sites 3, 11 and 12 will be subject to an inundation of water into these marine sectors; sites 3 and 11 will then be characterized by a relatively constant sea level until the end of the global deglaciation. The inset to Fig. 3A shows the sea level curve for an additional site located between sites 9 and 12 (red dot in Fig. 2B). The prediction at this site is characterized by a moderate, ∼25 m sea-level

fall over the period between the time when the site becomes ice free to the end of the deglaciation phase. Thus, it is possible for a site to record a sea-level fall across a globally synchronous Snowball deglaciation phase. However, such records would be localized to regions characterized by major ice cover that disappeared before the end of the deglaciation phase. The sea-level records at such sites, as well as any other location subject to large near-field deformational and gravitational effects, will not provide a robust estimate of globally averaged sea-level rise. Sea level during the post-deglacial phase continues to change due to ongoing viscoelastic deformation of the solid Earth (Figs. 2C, 2D and 3A). Fig. 2C focuses, for the purpose of illustration, on the first 10 kyr after deglaciation (Fig. 2D is a masked version of the same field). Post-glacial uplift of the crust drives sea-level fall over regions once covered by ice, while subsidence of the crust at the periphery of these zones of deglaciation results in a zone of sealevel rise surrounding the continents. Both regions have a peak amplitude of ∼500 m over the 10 kyr time window (Fig. 2C). At further distance from the centers of deglaciation, sea level is predicted to fall as water migrates from these regions to fill the accommodation space created by the subsidence peripheral to the ancient ice cover (Mitrovica and Milne, 2002). Sites 9 and 10 are located in a zone of peripheral subsidence, and are characterized by an additional sea-level rise of 250 m and 100 m, respectively, over the 10 kyr time window (Fig. 3A). In contrast, sites 3, 7, 11 and 12, are all in zones of post-glacial uplift and are predicted to experience sea-level regression of 200–500 m over this interval (Fig. 3A). Finally, site 17, located at the hinge point between zones of uplift and subsidence, shows relatively little change in sea level subsequent to the deglaciation phase (Fig. 3A). A Snowball deglacation of duration 2 kyr would suggest active, contemporaneous melting of all continental ice sheets. Nevertheless, even within such a rapid event, it is possible that the deglaciation was not perfectly synchronous and that localized melting began on one or more ice sheets (Hoffman and Macdonald, 2010). To investigate this possibility, we ran a simulation in which the northeast quadrant of the ice sheet covering continents E and F, located near the equator, melted over a 2 kyr deglaciation phase, while the melting of all other ice cover (Fig. 1B) took place only throughout the second half of this interval. (That is, melting of all but the localized zone was delayed relative to the timing adopted in the first simulation.) The total ESL change associated with this scenario is 35 m after 1 kyr (all sourced from the localized melt) and ∼1 km after 2 kyr (Figs. 3B, 4). Fig. 4C shows the total sea-level change 1 kyr (dashed red line) and 2 kyr (solid red line) into the model deglaciation phase along the transect shown in Fig. 4A. Results for the synchronous simulation are also shown for comparison (blue lines). Prior to the onset of global melting, sea level is predicted to fall ∼300 m at the center of the localized zone of melt (Fig. 4), an amplitude that is about an order of magnitude greater than the ESL change associated with the event, and ∼100 m at the edge of the continent (which is close to site 9 on Fig. 1A). At the end of the deglaciation phase, both simulations predict a net sea-level rise along the entire profile. Fig. 3B shows the full time history of relative sea-level change at site 9 for both the synchronous and asynchronous melt simulations. It is clear from this figure that gravitational and deformation effects of the early, localized melt dominate the ESL change associated with the event, and drive an early sea-level fall at site 9. However, once the remaining ice cover begins to disappear in the second half of the deglaciation phase, this early trend is reversed by the global ESL rise. Next, we investigate the sensitivity of the results to the adopted duration of deglaciation. Fig. 5 shows relative sea-level histories at four representative sites for which the duration of synchronous deglaciation varies from 2 kyr (as in the original simulation;

J.R. Creveling, J.X. Mitrovica / Earth and Planetary Science Letters 399 (2014) 74–85

79

Fig. 4. Predicted post-glacial sea-level change (in meters) for a scenario in which the volume of grounded ice cover in the northeast sector of the ice sheet covering continents E and F in Fig. 1A decreases linearly in time (from 100–0%) over the first 2 kyr of the simulation, and the volume of remaining global ice cover in Fig. 1B decreases linearly over the second half of this time window. In all cases ice sheets are constrained to maintain an equilibrium ice thickness profile. (A, B) Computed sea-level change over the first 1 kyr of the simulation in which 50% of the ice cover over the northeast sector of continent E/F in Fig. 1 melts (and the remaining ice cover is unchanged). Frame (A) includes a mask for contemporaneous continental paleogeography (Fig. 1A), while the inset in frame (B) does not include this mask. (C) Computed sea-level change across the transect shown by the dashed white line in frame (A) at the mid-point of the deglaciation phase (dashed red line), as in frames (A) and (B), and at the end of the 2 kyr deglaciation phase (solid red line). The blue lines are analogous predictions for the case of a globally synchronous deglaciation phase of 2 kyr (i.e., the simulation illustrated in Fig. 2). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 3A) to 200 kyr. As the duration of the deglaciation phase increases, the amplitude of post-deglacial sea-level change decreases. For instance, post-deglacial uplift and sea-level fall at site 7, and post-deglacial (peripheral bulge) subsidence and sea-level rise at site 10, decrease ∼95% as the duration of the model deglaciation phase increases from 2 kyr to 200 kyr. The amplitude of post-deglacial sea-level change decreases because more viscoelastic post-glacial adjustment occurs during the deglaciation phase when the duration of that phase is increased; as a result, the system evolves closer to isostatic equilibrium by the end of the deglaciation phase, thus minimizing post-deglacial sea-level change. The predictions at sites 3 and 11 show a somewhat more complex sensitivity to the timescale of deglaciation. In particular, the sign of the predicted post-deglacial sea-level trends changes as the duration of the deglaciation is increased. For rapid deglaciation scenarios, the post-glacial adjustment at these sites is dominated by uplift (and sea-level fall) due to local ice unloading, whereas for slow deglaciation scenarios, the adjustment becomes progressively more dominated by eustatic sea-level rise associated with global melting and (in the case of site 11) subsidence driven by melting of surrounding ice sheets. We further discuss the implication of this sign change in Section 5.

Finally, we consider the sensitivity of our predictions to variations in the adopted radial profile of mantle viscosity. To this point, the simulations have been based on an Earth model with a lower mantle viscosity of 5 × 1021 Pa s, which is at the lower bound of inferences from analyses of ice-age sea-level histories (Lambeck et al., 1998; Mitrovica and Forte, 2004). In Fig. 6 we reproduce the calculations for sites 10 and 11 in Fig. 5, and compare these with results from simulations in which the lower mantle viscosity is increased to 2 × 1022 Pa s. (For clarity, we have omitted from the figure predictions in which the model deglaciation phase has a duration of either 5 kyr or 10 kyr.) As we argued above, the predicted syn-deglacial sea-level change associated with our standard run, in which the deglaciation phase is very rapid (duration of 2 kyr), is insensitive to the adopted viscosity model. However, the predicted post-deglacial sea-level change for the run does show sensitivity to this choice, with a slower return to equilibrium as the adopted lower mantle viscosity is increased. As the duration of the deglaciation phase is increased from 2 kyr to 20 kyr, the sensitivity of both the syn- and post-deglacial sea-level predictions to the lower mantle viscosity increases. In particular, runs with the higher viscosity value show a more muted post-deglacial sea-level variation and longer adjustment time scales. Finally, the predictions become progressively

80

J.R. Creveling, J.X. Mitrovica / Earth and Planetary Science Letters 399 (2014) 74–85

Fig. 5. Computed sea-level histories for four sites (as labeled) for scenarios in which a globally synchronous deglaciation of grounded ice cover in Fig. 1B occurs over a duration of (from left to right on each frame) 200 kyr, 100 kyr, 50 kyr, 20 kyr, 10 kyr, 5 kyr and 2 kyr. The 2 kyr case is identical to the scenario considered in Fig. 2 and the sea-level predictions for the associated sites is given in Fig. 3A. Dotted segments of the curves denote times when the site was covered by grounded ice.

less sensitive to the lower mantle viscosity as the duration of the deglaciation phase is increased above 20 kyr. These various trends reflect the Maxwell time of the two Earth models, ∼2 kyr for a lower mantle viscosity of 5 × 1021 Pa s, and ∼10 kyr for a viscosity of 2 × 1022 Pa s. Syn-deglacial sea-level histories should thus become less sensitive to this range of lower mantle viscosity for deglaciation time scales less than 2 kyr and greater than 20 kyr. In any event, it is clear from Fig. 6 that the general form of the sea-level predictions, and our conclusions based on the results in Figs. 2–5, are not sensitive to plausible variations in the lower mantle viscosity. 4. Comparing fingerprints of Snowball and modern ice sheet collapse

Fig. 6. Computed sea-level histories for two sites (as labeled) for scenarios in which a globally synchronous deglaciation of grounded ice cover in Fig. 1B occurs over a duration of (from left to right on each frame) 200 kyr, 100 kyr, 50 kyr, 20 kyr and 2 kyr. Predictions are based on Earth models with a lower mantle viscosity of (solid lines; as in Fig. 5) 5 × 1021 Pa s, and (dotted lines) 2 × 1022 Pa s.

Our predictions of sea-level change across a Snowball deglaciation may be compared with fingerprints of sea-level change following rapid melting of modern polar ice sheets and glaciers. The latter are characterized by a near-field sea-level fall that is an order of magnitude greater than the eustatic sea-level rise associated with the melt event (e.g., Clark and Lingle, 1977; Mitrovica et al., 2001). Our prediction of sea-level change in the case where Snowball deglaciation is initiated with a phase of highly localized melting shows a similar, order of magnitude difference between the peak near-field sea-level fall and the ESL rise over the first 1 kyr of the model deglaciation (Fig. 4C, dashed red line), prior to the onset of global deglaciation. However, this similarity in the pattern of deglacial sea-level change does not extend to the globally synchronous melt simulation where pronounced sea-level rise occurs along the majority of near-field margins (Fig. 2A). There are two reasons for this difference. First, the volume of any of the individual ice sheets in Fig. 1B is significantly smaller than the global ice volume, thus the eustatic signal arising from a globally synchronous deglaciation is greater than the signal from near-field

J.R. Creveling, J.X. Mitrovica / Earth and Planetary Science Letters 399 (2014) 74–85

81

traction) associated with the deglaciation of a continental-scale ice sheet. Determination of the eustatic signal would depend on the scale of the dynamic signal, but for the size of the modeled ice sheets considered here (Fig. 1B), a eustatic signal greater than ∼500 m would produce syn-deglacial transgression at sites adjacent to an open ocean. Therefore, robust observations of globally distributed syn-deglacial transgression could provide a stratigraphic test for the deglaciation of ‘Snowball’-scale ice volumes. 5. Mechanisms for sea-level regression during and after Snowball deglaciation

Fig. 7. Total sea-level change across the model Snowball deglaciations of duration 2 kyr at the 20 sites shown in Fig. 1A. (A, B) Results for the synchronous (Figs. 2B, 3A) and asynchronous (Fig. 4) deglaciation simulations, respectively.

gravitational and deformational effects (which oppose the eustatic sea-level rise) associated with melting of the nearby ice sheet. Second, the amplitude of the peak sea-level fall driven by near-field effects decreases with the area of the ice sheet, and all but the ice sheet on continent L (Fig. 1A) are much larger than the modern Greenland or West Antarctic Ice Sheets. The net effect is a more muted near-field signal in Fig. 2A than one might expect on the basis of published predictions of modern sea-level fingerprints. However, as we noted in the context of Figs. 2 and 3, even in the case of the synchronous Snowball deglaciation, the net sealevel signal from near-field gravitational and deformation signals is non-negligible. Fig. 7 shows the net sea-level change across the 2 kyr deglaciation at all 20 sites shown in Fig. 1A for both the synchronous and asynchronous deglaciation simulations. The total sea-level change across the deglaciation is nearly the same for both simulations (see also Fig. 4C). At sites previously under ice cover (3–5 and 11–13), near-field effects compensate the eustatic signal at a level of 70–100% (i.e., the predictions show a sea-level rise of 0–30% of the eustatic value). The remaining 14 sites, all at the (open ocean) edge of the ice cover, show consistency in the predicted sea-level rise, regardless of the location of the site relative to the eight ice sheets of various spatial scale included in the simulations. At these sites, local sea-level histories exhibit a sea-level rise of ∼50–70% of the eustatic value. Hoffman and Macdonald (2010) suggest that the near-field gravitational and deformational effects would be an order of magnitude smaller than the ESL value associated with a globally synchronous Snowball deglaciation. Our results, in contrast, indicate that the magnitude of these effects is ∼30–100% of the ESL change. We note that the cap successions described in the geological literature were deposited at near-field sites, and thus using sea-level histories inferred from these records as a measure of the total global (eustatic) sea-level rise (or, equivalently, the total volume of ice) associated with the Snowball deglaciation (i.e., assuming near-field effects are negligible) would lead to a significant underestimate of the eustatic sea-level change. Finally, if robust, observations of globally distributed syndeglacial transgression over a formerly glaciated margin would provide a strong argument for a deglaciation far more extensive than the Pleistocene. Indeed, such observations would require a eustatic sea-level rise of a magnitude sufficient to overcome the dynamic effects (elastic rebound and the loss of gravitational at-

Hoffman et al. (2007) and Hoffman and Macdonald (2010) present geological evidence for an early, regional regression punctuating the overall cap dolostone transgressive sequence tract. If robust, how can we explain this evidence for regional sea-level fall in the context of our modeling? Our model predictions suggest at least three possibilities. First, as originally proposed by Hoffman and Macdonald (2010), near field gravitational and deformational effects during localized (asynchronous) melt would give rise to regional regression. This is the model explored in Figs. 3B (dashed line) and 4. It is important to note that these localized melt events do not have to occur at the initiation of the deglaciation, as in our asynchronous simulation. Regression can occur during any interval in which local melt takes place during a global stasis in melting, that is, during intervals in which the ESL change associated with the melting of distant ice sheets is small compared to the near-field effects associated with the localized zone of collapse. Of course, if the localized melt events occur at some point within the deglaciation, rather than at its start, viscous effects associated with any previous melting may significantly contaminate the predictions in Figs. 4 and 3B (an example of such a viscous signal is shown in Fig. 2C). Our numerical predictions in Fig. 3A suggest a second possible cause for inferred episodes of regional regression in cap dolostones: sea-level fall may reflect the post-glacial uplift of sites in regions of significant ice loss. This behavior is evident, for example, at sites 3, 7, 11 and 12 in Fig. 3A, which are all situated well within the ice margin at maximum glaciation. In this particular simulation, where a globally synchronous deglaciation is assumed, the post-glacial sea-level fall is generally only evident subsequent to the completion of the globally synchronous deglaciation, when the eustatic sea-level rise associated with the global deglaciation has ended. However, we note that a regression driven by post-glacial uplift is not necessarily limited to the period after deglaciation has ended. The only requirement is that local post-glacial uplift be greater than the ESL rise associated with distant melt events. The inset to Fig. 3A provides an example of a site where post-glacial uplift is sufficient to counter an ongoing eustatic rise. A modern analogy to this behavior is the post-glacial emergence of Hudson Bay following the Last Glacial Maximum, as reflected in flights of raised beach terraces which formed from the moment the region became ice free until the present day. All predictions characterized by post-glacial uplift in Fig. 3A (sites 3, 7, 11 and 12) and its inset show sea-level histories that end in regression, but this is not necessarily so. If, during the global Snowball deglaciation phase (i.e., sometime within the 2 kyr deglaciation phase modeled in Figs. 2–4), the global meltwater flux slowed, then these sites would be characterized by a period of sealevel fall. In this case, a re-establishment of global scale melting would then once again initiate a global scale transgression. This explanation replaces the model of Hoffman and Macdonald (2010), which attributes regression to localized melting during an otherwise global-scale hiatus in melt, with a model in which the regression is explained by post-glacial uplift during the global-scale hiatus. A distinction of the latter model is that it would provide

82

J.R. Creveling, J.X. Mitrovica / Earth and Planetary Science Letters 399 (2014) 74–85

Fig. 8. Comparison of Marinoan post-glacial sea-level histories as determined from the geologic record of cap succession to numerical predictions for sea-level change during and after a modeled synchronous deglaciation of duration 2–200 kyr. (A) Reconstructed relative sea-level histories from lithofacies variation in post-glacial ‘cap’ stratigraphic successions from the Mackenzie Mountains, Yukon, Canada (site 16; adapted from the Ravensthroat River section, Hoffman and Halverson, 2011), the Gariep Belt, Kalahari craton, Namibia (site 9; adapted from the Namaskluft Camp section, Hoffman and Macdonald, 2010), and the Adelaide Rift Complex, Flinders Range, Australia (site 5; an idealized composite section from stratigraphic columns in Rose and Maloof, 2010 and Rose et al., 2013). The vertical axis represents meters of cap stratigraphy (younging-upwards). Post-glacial time is likely distributed differently throughout each section. Reconstructed relative sea-level curves (rise towards the ‘+’ symbol, fall towards the ‘–’ symbol) show inundation (sea level rise) across the glacial unit–cap carbonate contact (0 m) and any subsequent change in sea level based on lithofacies variation described in the references above. The dashed portion of the curve indicates the sea-level fall (and resulting sequence boundary) recorded in overlying strata (not pictured, and not to scale). (B) Sea-level histories for model sites 5, 9, and 16 (as labeled; see Fig. 1A) computed for simulations with a globally synchronous deglaciation of duration (from left to right on each frame) 200 kyr, 100 kyr, 50 kyr, 20 kyr and 2 kyr (as in Fig. 6). Dotted segments of the curves denote times when the site was covered by grounded ice. Insets reproduce the reconstructed sea-level history as determined from the cap successions at each site (that is, identical to those in Fig. 8A, but rotated such that time runs in the same direction as the abscissas in Fig. 8B).

a simple explanation for regression events measured at disparate sites. We note that both localized melt or a global cessation in melting may be difficult to reconcile with geochemical signatures for high pCO2 within Marinoan cap carbonates (Bao et al., 2008; Kasemann et al., 2010). However, Sansjofre et al. (2011) and Cao and Bao (2013) discuss the possibility of a post-Snowball atmosphere with a low concentration of carbon dioxide. Finally, the simulations that explored the sensitivity of the predictions to the duration of the deglaciation phase (Fig. 5) provide a third possible explanation for geological evidence of a period of regional regression punctuating an overall transgressive sequence tract. In particular, a slow deglaciation event can lead, at some sites (e.g., 3 and, particularly, 11), to a transition from an early period of sea-level fall dominated by post-glacial uplift to a later phase in which a net sea-level rise is due to the combined effects of a eustatic sea-level rise due to ongoing global deglaciation and either slowing post-glacial uplift or peripheral subsidence associated with melting on nearby continents. The transition from uplift due to local melting to peripheral bulge subsidence due to nearby melting is evident at site 11 on Fig. 5. The predictions for this site show a sea-level rise after the deglaciation phase has ended (i.e., after time 0) in the case where the modeled duration of the deglaciation phase is 50 kyr or longer. How might one discriminate between these various scenarios for a syn-deglacial sea-level regression followed by transgression? As we have noted, whether a sea-level history of this type is in-

ferred at an isolated site or a group of sites may help to distinguish between a local, asynchronous melt event or post-glacial uplift of a suite of previously ice covered regions during a period of slowing global ice melting. Furthermore, a protracted sea-level regression suggests post-glacial uplift rather than gravitational and deformational effects associated with local, asynchronous ice melting, and it provides a measure of the minimum duration of the deglaciation phase. Compare, for example, the dashed line in Fig. 3B, which is characterized by a regression of 1 kyr duration, with the leftmost curve of the site 11 results of Fig. 5, which has a regression that persists (after the site becomes ice-free) for over 20 kyr. 6. Discussion Given the uncertainty in Cryogenian–Ediacaran paleogeography and the adopted ice and Earth models, we cautioned in Section 2 that the sea-level histories computed for the twenty sites shown in Fig. 1 were illustrative, and not meant to reproduce the sea-level histories inferred for specific geologic sites. Nevertheless, it is instructive to explore whether the variability in our predicted sea-level histories encompasses the range of variability in post-glacial sea-level change observed in the geologic record of cap successions. To this end, Fig. 8 shows sea-level reconstructions inferred from lithofacies variation in cap successions from three regions with a well-developed, three-dimensional, post-glacial stratigraphic framework: the Mackenzie Mountains,

J.R. Creveling, J.X. Mitrovica / Earth and Planetary Science Letters 399 (2014) 74–85

Canada (site 16; James et al., 2001; Hoffman and Halverson, 2011; Macdonald et al., 2013), the Gariep Belt, Namibia (site 9; Hoffman and Macdonald, 2010), and the Adelaide Rift Complex, Flinders Range, Australia (site 5; Rose and Maloof, 2010; Rose et al., 2013). The figure also shows numerical predictions of sea-level change at these three locations for synchronous 2–200 kyr duration deglaciation scenarios. Broadly, these localities reveal syn-deglacial transgression (site 16), or syn-deglacial transgression following regression within either the diamictite (site 5) or the basal cap carbonate (site 9). For modeled synchronous Marinoan Snowball deglaciations of 2–200 kyr duration, numerical predictions for syn-deglacial sealevel change at site 16 show sustained, high-amplitude transgression. The longer-duration (>20 kyr) deglaciation scenarios also predict a residual, low-amplitude (up to ∼10 m) sea-level rise that extends ∼20–30 kyr into the post-deglacial phase (Fig. 8B). Beyond this time, all simulations are characterized by a gradual, long-term regression. Hoffman and Halverson (2011) concluded that the Ravensthroat, Hayhook, and lower Sheepbed formations (below the maximum flooding interval; see Macdonald et al., 2013) record syn-deglacial transgression alone (Fig. 8A), an interpretation consistent with the sea-level predictions for the shorter-duration deglaciation scenarios (<20 kyr; Fig. 8B). However, given that the transgression predicted at site 16 need not be confined to the duration of deglaciation, ruling out moderate post-deglacial deposition may be difficult. In this regard, if the Hayhook and/or lower Sheepbed formations were to record post-deglacial deposition, then a longer-duration (>20 kyr) deglaciation scenario could instead explain the Mackenzie Mountain sea-level reconstruction (Fig. 8). To explore this possibility further, we ask how might our model for Snowball deglaciation be revised to yield a continued, highamplitude post-deglacial sea-level rise at site 16? The results in Fig. 3A indicate that site 10, at the northern margin of continent H (Fig. 1A), is predicted to have a total sea-level rise of ∼900 m extending through both the syn-deglacial and post-deglacial phases, where the former is dominated by the global ESL rise while the latter is due to the fact that the site experiences peripheral bulge subsidence associated with melting of large ice sheets to the east and west (Fig. 2C). Therefore, in order to predict a high-amplitude post-deglacial sea-level rise at site 16, the model paleogeography and/or ice geometry must be tuned to bring the site into a region of peripheral subsidence. Test calculations we have performed indicate that this is easily accomplished by: (1) altering the geometry of the northern margin of continent J so that site 16 is placed on this margin, rather than being inland of the margin (as in Fig. 1) and/or (2) moving southward the maximum northern perimeter of the ice sheet that covered this continent. In this regard, we note that the results for site 10 in Fig. 5 indicate that a high-amplitude post-deglacial sea-level transgression, as reconstructed for site 16 in Fig. 8A, suggests a relatively fast (<20 kyr) deglaciation phase. Next, the predictions for synchronous 2–200 kyr sea-level change in Namibia (site 9) agree with a reconstructed largeamplitude syn-deglacial transgression – but do not agree with sedimentological observations for a sea-level fall preceding this transgression – in the Holgat Formation (Fig. 8). In Section 5 we discussed several plausible mechanisms that could explain this regression (Figs. 3B, 4). For example, an asynchronous phase of melting localized to the region close to site 9 would introduce an early, short-lived regression, as shown in Fig. 3B (Hoffman and Macdonald, 2010). Alternatively, we could revise the model of paleogeography and ice cover so as to locate site 9 further inland or under substantially more initial ice cover (Fig. 1). In this case, in analogy with the predictions for site 11 (see Fig. 5), site 9 would experience an initial regression associated with a dominant contribution from post-glacial uplift, followed by a transgression as this

83

uplift was overtaken by the global ESL rise. As we discussed above, these two scenarios could be distinguished by constraints on the duration of the initial regression. Finally, consider the geological record and predictions for site 5 (Fig. 8). The synchronous 50–200 kyr duration deglaciation scenarios predict well the syn-deglacial sea-level fall and rise reconstructed by Rose et al. (2013) for the Elatina–Nuccaleena formations. We note that the same numerical predictions also predict the initial syn-deglacial inundation of the site as it became icefree. In contrast, shorter-duration deglaciation scenarios (<50 kyr) do not reproduce the complexity of sea-level change inferred for this site (Fig. 8). We re-emphasize that tuning model parameters – including paleogeography, ice sheet geometry and magnitude, the duration and synchronicity of melt, and mantle viscosity and rheology (the Earth model) – could improve site-specific predictions of post-glacial sea-level change, as well as capture any sea-level variability not presently reconstructed from a cap succession. Nonetheless, this preliminary comparison indicates that the variability in sea-level histories predicted from our modeling encompasses the range of post-glacial sea-level change inferred from the geological record. Additional regional post-glacial sea-level histories reconstructed from cap successions that preserve a shelf–slope transition would further constrain the geometry, magnitude, and deglaciation sequence of Cryogenian ice sheets. If δ13 Ccarb chemostratigraphic variation within cap carbonates reflects the global, time-evolving chemical composition of postglacial oceans (Hoffman et al., 1998), then the record of δ13 Ccarb variation at a specific site will depend on the post-glacial sea level history (Hoffman et al., 2007). For example, a site that experienced a sea-level stillstand and yet maintained submarine deposition throughout the post-glacial phase should record the full sigmoidal trajectory of post-glacial δ13 Ccarb evolution as defined by Hoffman et al. (2007). Alternatively, delayed transgression (due to early regression or post-glacial uplift) would result in partial capture of the negative limb of the δ13 Ccarb excursion. We note, however, that the δ13 Ccarb chemostratigraphic variation may not reflect steady-state behavior (Hoffman and Schrag, 2002; Higgins and Schrag, 2003), in which case alternative methodologies must be sought to determine the relative timing of deposition of the globally-distributed cap successions. 7. Conclusions If the Marinoan Snowball deglaciation occurred rapidly, then the geological observation of widespread syn-deglacial transgression over glaciated margins speaks to the magnitude of the glaciation: by the end of the deglaciation phase, the eustatic contribution from the melting of far-field ice sheets was (generally) larger than the local sea-level fall driven by near-field effects, a consequence of the immense volume of ice distributed around the Cryogenian globe. However, the results in Figs. 2B and 3–6 demonstrate that syn-deglacial transgression inferred from the geological record at specific sites is neither an accurate measure of the eustatic sealevel change associated with the global deglaciation, nor a global phenomena. Alternative syn-deglacial sea-level histories, including local sea-level fall and sea-level stillstand, may have coincided with transgression elsewhere. In principle, these complex sea-level histories should be manifest in the stratigraphic architecture of cap successions; however such histories may be largely absent (due to erosion, or hiatus) or difficult to read from the rock record, perhaps perpetuating the notion of global transgression. Nonetheless, these more complex syn-deglacial geological records would be entirely consistent with a globally-averaged, high-amplitude eustatic sea-level rise resulting from a pan-glacial Snowball deglaciation.

84

J.R. Creveling, J.X. Mitrovica / Earth and Planetary Science Letters 399 (2014) 74–85

However, we caution that the geological observation of transgression atop a glacial deposit (or its unconformity) need not reflect syn-deglacial deposition. While some cap dolostones record evidence for syn-deglacial deposition such as ice-rafted debris (e.g., Hoffman and Macdonald, 2010; Rose and Maloof, 2010; Petterson et al., 2011), others do not (e.g., James et al., 2001; Hoffman and Halverson, 2011), and in the latter case the supposition of syn-deglacial deposition is a conceptual prediction of the diachronous deglaciation model (Hoffman et al., 2007) rather than a geological observation. We emphasize this point because our model predictions for both short and long duration deglaciation scenarios (2–200 kyr) indicate that at some localities sea level will continue to rise post-deglaciation, allowing for the possibility that cap dolostones (or, more broadly, the transgressive systems tract) record syn-deglacial deposition, post-deglacial deposition (Fig. 3A), or both, and need not be confined to the deglaciation of Snowball ice sheets. This conclusion, based on the physics of sea-level change, supports geochemical models that require a sustained alkalinity flux over timescales longer than the shortest hypothesized duration of Snowball melt (103 –104 yr) in order to deposit cap dolostones (Le Hir et al., 2009). Taken together, our numerical simulations indicate that strata associated with the Marinoan Snowball deglaciation and its aftermath may reflect a more complicated history of transgression and regression than previously considered. (We note that this complexity may be further increased if hydrological processes led to syn-deglacial fluctuations in ice volumes; Halverson et al., 2004.) This conclusion is not dependent on the extent of sea ice during the Snowball glaciation since this would not have had an impact on sea level. Thus, our conclusions are not sensitive to the ongoing debate as to whether the Marinoan glaciation was a global hard “snowball” or, for example, a “Jormungand” state with localized areas of open ocean (Hoffman, 2009; Abbot et al., 2011). These complex post-glacial sea-level histories can inform models of postglacial carbonate and silicate weathering (through the time-history of exposed shelf area), and, thus, post-glacial alkalinity fluxes, atmospheric CO2 concentration, and the δ13 Ccarb isotopic evolution of cap carbonate successions (Higgins and Schrag, 2003). Acknowledgements The Agouron Institute Geobiology Postdoctoral Fellowship (JRC), Harvard University (JXM), and the Earth System Evolution Program of the Canadian Institute for Advanced Research (JXM) provided support for this research. We thank David Evans for a highresolution image of the Ediacaran paleogeography, Eric Morrow for digitizing this image, Paul Hoffman for comments on several versions of the manuscript, and Adam Maloof and an anonymous reviewer for insightful comments that significantly improved this submission. References Abbot, D.S., Voigt, A., Koll, D., 2011. The Jormungand global climate state and implications for Neoproterozoic glaciations. J. Geophys. Res., Atmos. 116. http:// dx.doi.org/10.1029/2011JD015927. Allen, P.A., Hoffman, P.F., 2005. Extreme winds and waves in the aftermath of a Neoproterozoic glaciation. Nature 453, 123–127. http://dx.doi.org/10.1038/ nature03176. Bao, H., Lyons, J.R., Zhou, C., 2008. Triple oxygen isotope evidence for elevated CO2 levels after a Neoproterozoic glaciation. Nature 453, 504–506. http://dx.doi.org/ 10.1038/nature06959. Bosak, T., Mariotti, G., Macdonald, F.A., Perron, J.T., Pruss, S.B., 2013. Microbial sedimentology of stromatolites in Neoproterozoic cap carbonates. In: Bush, A.M., Pruss, S.B., Payne, J.L. (Eds.), Ecosystem Paleobiology and Geobiology. In: Paleontological Special Papers, vol. 19. Paleontological Society, pp. 51–77. Cao, X., Bao, H., 2013. Dynamic model constraints on oxygen-17 depletion in atmospheric O2 after a snowball Earth. Proc. Natl. Acad. Sci. USA 110 (36), 14546–14550. http://dx.doi.org/10.1073/pnas.1302972110.

Clark, J.A., Lingle, C.S., 1977. Future sea-level changes due to West Antarctic ice sheet fluctuations. Nature 269, 206–209. http://dx.doi.org/10.1038/269206a0. Dziewonski, A.M., Anderson, D.L., 1981. Preliminary reference Earth model. Phys. Earth Planet. Inter. 25, 297–356. http://dx.doi.org/10.1016/0031-9201(81)900467. Evans, D.A.D., 2000. Stratigraphic, geochronological, and paleomagnetic constraints upon the Neoproterozoic climatic paradox. Am. J. Sci. 300, 347–433. http:// dx.doi.org/10.2475/ajs.300.5.347. Halverson, G.P., Maloof, A.C., Hoffman, P.F., 2004. The Marinoan glaciation (Neoproterozoic) in northeast Svalbard. Basin Res. 16, 297–324. http://dx.doi.org/ 10.1111/j.1365-2117.2004.00234.x. Higgins, J.A., Schrag, D.P., 2003. Aftermath of a snowball Earth. Geochem. Geophys. Geosyst. 4, 20. http://dx.doi.org/10.1029/2002GC000403. Hoffman, P.F., 2009. Pan-glacial – a third state in the climate system. Geol. Today 25, 107–114. http://dx.doi.org/10.1111/j.1365-2451.2009.00716.x. Hoffman, P.F., 2011. Strange bedfellows: glacial diamictite and cap carbonate from the Marinoan (635 Ma) glaciation in Namibia. Sedimentology 58, 57–119. http:// dx.doi.org/10.1111/j.1365-3091.2010.01206.x. Hoffman, P.F., Halverson, G.P., 2011. Neoproterozoic glacial record in the Mackenzie Mountains, northern Canadian Cordillera. In: Mem. Geol. Soc. Lond., vol. 36, pp. 397–412. Hoffman, P.F., Halverson, G.P., Domack, E.W., Husson, J.M., Higgins, J.A., Schrag, D.P., 2007. Are basal Ediacaran (635 Ma) post-glacial “cap dolostones” diachronous? Earth Planet. Sci. Lett. 258, 114–131. http://dx.doi.org/10.1016/j.epsl.2007.03.032. Hoffman, P.F., Kaufman, A.J., Halverson, G.P., Schrag, D.P., 1998. A Neoproterozoic snowball earth. Science 281, 1342–1346. http://dx.doi.org/10.1126/ science.281.5381.1342. Hoffman, P.F., Li, Z.X., 2009. A palaeogeographic context for Neoproterozoic glaciation. Palaeogeogr. Palaeoclimatol. Palaeoecol. 277, 158–172. http://dx.doi.org/ 10.1016/j.palaeo.2009.03.013. Hoffman, P.F., Macdonald, F.A., 2010. Sheet-crack cements and early regression in Marinoan (635 Ma) cap dolostones: regional benchmarks of vanishing ice-sheets? Earth Planet. Sci. Lett. 300, 374–384. http://dx.doi.org/10.1016/ j.epsl.2010.10.027. Hoffman, P.F., Schrag, D.P., 2002. The snowball Earth hypothesis: testing the limits of global change. Terra Nova 14, 129–155. http://dx.doi.org/10.1046/ j.1365-3121.2002.00408.x. Hyde, W.T., Crowley, T.J., Baum, S.K., Peltier, W.R., 2000. Neoproterozoic ‘snowball Earth’ simulations with a coupled climate/ice-sheet model. Nature 405, 425–429. http://dx.doi.org/10.1038/35013005. James, N.P., Narbonne, G.M., Kyser, T.K., 2001. Late Neoproterozoic cap carbonates: Mackenzie Mountains, northwestern Canada: precipitation and global glacial meltdown. Can. J. Earth Sci. 38, 1229–1262. http://dx.doi.org/10.1139/ cjes-38-8-1229. Jerolmack, D.J., Mohrig, D., 2005. Palaeoclimatology: formation of Precambrian sediment ripples. Nature 436 (7049), E1. http://dx.doi.org/10.1038/nature04025. Kasemann, S.A., Prave, A.R., Fallick, A.E., Hawkesworth, C.J., Hoffmann, K.H., 2010. Neoproterozoic ice ages, boron isotopes, and ocean acidification: implications for a snowball Earth. Geology 38, 775–778. http://dx.doi.org/10.1130/G30851.1. Kendall, R.A., Mitrovica, J.X., Milne, G.A., 2005. On post-glacial sea level – II. Numerical formulation and comparative results on spherically symmetric models. Geophys. J. Int. 161, 679–706. http://dx.doi.org/10.1111/j.1365-246X.2005.02553.x. Kennedy, M.J., 1996. Stratigraphy, sedimentology, and isotopic geochemistry of Australian Neoproterozoic postglacial cap dolostones: deglaciation, δ13 C excursions, and carbonate precipitation. J. Sediment. Res. 66, 1050–1064. Kilner, B., Niocaill, C., Brasier, M., 2005. Low-latitude glaciation in the Neoproterozoic of Oman. Geology 33, 413–416. http://dx.doi.org/10.1130/G21227.1. Lambeck, K., Smither, C., Johnston, P., 1998. Sea-level change, glacial rebound and mantle viscosity for northern Europe. Geophys. J. Int. 134, 102–144. http:// dx.doi.org/10.1046/j.1365-246x.1998.00541.x. Le Hir, G., Donnadieu, Y., Goddéris, Y., Pierrehumbert, R.T., Halverson, G.P., Macouin, M., Nédélec, A., Ramstein, G., 2009. The snowball Earth aftermath: exploring the limits of continental weathering processes. Earth Planet. Sci. Lett. 277, 453–463. http://dx.doi.org/10.1016/j.epsl.2008.11.010. Li, Z.X., Evans, D.A.D., Halverson, G., 2013. Neoproterozoic glaciations in a revised global palaeogeography from the breakup of Rodinia to the assembly of Gondwanaland. Sediment. Geol. 294, 219–232. http://dx.doi.org/10.1016/ j.sedgeo.2013.05.016. Liu, Y., Peltier, W.R., 2013. Sea level variations during snowball Earth formation: 1. A preliminary analysis. J. Geophys. Res. 118, 4410–4424. http://dx.doi.org/10.1002/ jgrb.50293. Macdonald, F.A., Strauss, J.V., Sperling, E.A., Halverson, G.P., Narbonne, G.M., Johnston, D.T., Kunzmann, M., Schrag, D.P., Higgins, J.A., 2013. The stratigraphic relationship between the Shuram carbon isotope excursion, the oxygenation of Neoproterozoic oceans, and the first appearance of the Ediacara biota and bilaterian trace fossils in northwestern Canada. Chem. Geol. 362, 250–272. http:// dx.doi.org/10.1016/j.chemgeo.2013.05.032. Meert, J.G., Torsvik, T.H., 2004. Paleomagnetic constraints on Neoproterozoic ‘Snowball Earth’ continental reconstructions. Geophys. Monogr. 146, 5–11.

J.R. Creveling, J.X. Mitrovica / Earth and Planetary Science Letters 399 (2014) 74–85

Mitrovica, J.X., Forte, A.M., 2004. A new inference of mantle viscosity based upon joint inversion of convection and glacial isostatic adjustment data. Earth Planet. Sci. Lett. 225, 177–189. http://dx.doi.org/10.1016/j.epsl.2004.06.005. Mitrovica, J.X., Milne, G.A., 2002. On the origin of late Holocene sea-level highstands within equatorial ocean basins. Quat. Sci. Rev. 21, 2179–2190. http://dx.doi.org/ 10.1016/S0277-3791(02)00080-x. Mitrovica, J.X., Milne, G.A., 2003. On post-glacial sea level: I. General theory. Geophys. J. Int. 154, 253–267. http://dx.doi.org/10.1046/j.1365-246X.2003.01942.x. Mitrovica, J.X., Tamisiea, M.E., Davis, J.L., Milne, G.A., 2001. Recent mass balance of polar ice sheets inferred from patterns of global sea-level change. Nature 409, 1026–1029. http://dx.doi.org/10.1038/35059054. Petterson, R., Prave, A.R., Wernicke, B.P., Fallick, A.E., 2011. The Neoproterozoic Noonday formation, Death Valley region, California. Geol. Soc. Am. Bull. 123, 1317–1336. http://dx.doi.org/10.1130/B30281.1.

85

Rose, C.V., Maloof, A.C., 2010. Testing models for post-glacial ‘cap dolostone’ deposition: Nuccaleena formation South Australia. Earth Planet. Sci. Lett. 296, 165–180. http://dx.doi.org/10.1016/j.epsl.2010.03.031. Rose, C.V., Maloof, A.C., Schoene, B., Ewing, R.C., Linnemann, U., Hofmann, M., Cottle, J.M., Budnick, A., 2013. The end-Cryogenian glaciation of South Australia. Geosci. Can. 40, 256–293. http://dx.doi.org/10.12789/geocanj.2013.40.019. Sansjofre, P., Ader, M., Trindade, R.I.F., Elie, M., Lyons, J., Cartigny, P., Nogueira, A.C.R., 2011. A carbon isotope challenge to the snowball Earth. Nature 478 (7367), 93–96. http://dx.doi.org/10.1038/nature10499. Torsvik, T.H., 2003. The Rodinia jigsaw puzzle. Science 300, 1379–1381. http:// dx.doi.org/10.1126/science.1083469. Trindade, R.I.F., Font, E., D’Agrella-Filho, M.S., Nogueira, A.C.R., Riccomini, C., 2003. Low-latitude and multiple geomagnetic reversals in the Neoproterozoic Puga cap carbonate, Amazon craton. Terra Nova 15, 441–446. http://dx.doi.org/10.1046/ j.1365-3121.2003.00510.x.