Introduction

The Snowball Earths (SBE) are events in Earth’s history when ice sheets covered most, if not all of the Earth’s surface (e.g. Cox et al. 2016; Donnadieu et al. 2004; Hoffman et al. 2017; Kirschvink 1992; Stern and Miller 2019). Physical geological evidence traces ice sheets, leaving the extent of ocean ice to be inferred from climate models, geochemical proxies, and extreme climatic hysteresis inferred from low-latitude ice sheet longevity. During this period, the surface biosphere survived in cryoconite meltwater-based ecosystems on ablative ice sheets and ice shelf surfaces, as well as terrestrial dry-valley lakes, streams, and soils in the equatorial zone of net ablation (Hawes et al. 2018; Hoffman 2016; Vincent et al. 2000). Eventually, volcanic outgassing causes atmospheric CO2 to reach sufficient levels to create a runaway greenhouse warming effect (Caldeira and Kasting 1992; Halverson et al. 2005; Hoffman et al. 1998; Kirschvink 1992; Walker et al. 1981). Carbonate δ13C (permil deviation of the 13C/12C from the Vienna Pee Dee Belemnite - VPDB) values in contemporaneous cap carbonates deposited following deglaciation reflect the global scale of the event (e.g. Hoffman et al. 2017 and references therein). Isotopic excursions captured in the post-SBE cap dolostones are traditionally considered to begin with a steep decline in δ13C to -6 ± 1 ‰ (Brasier and Shields 2000; Hoffman et al. 1998; Hoffman and Schrag 2002; Hoffman et al. 2017; Kaufman and Knoll 1995; Kaufman et al. 1997). The post SBE carbon isotope excursions (CIEs) likely represent 13C depleted carbonates recycled during SBE ocean buffering, and during deglaciation through weathering of carbonate loess exposed by ice-sheet recession (Higgins and Schrag 2003). As bicarbonate formation rates are accelerated and amplified from CO2 brought into the oceans following the glacial events, low δ13C values persist even post climatic recovery (Hoffman and Maloof 1999; Hoffman and Schrag 2002). The timescales to obtain equilibrium are difficult to quantify (Gernon et al. 2016 and references therein), but have been suggested to be on the order of 5 \(\cdot\) 104 years (Yang et al. 2017), although it is possible full equilibrium was never attained.

Oxygen isotopes provide another useful tool in interpreting the Snowball events, with δ18O (permil deviation of 18O/16O from VPDB) values of carbonates providing temperature, source strengths, reservoir exchange magnitudes and relative sizes of oxidized and reduced carbon (Jaffres et al. 2007; Veizer et al. 1999). Ocean water’s δ18O is modeled to have been as low as -13 ‰ during the Neoproterozoic, with values of δ18O of -5 to -10 ‰ observed in Maieberg Formation samples (Hoffman et al. 2021; Jaffres et al. 2007; Shields 2005). Of key interest with oxygen is the 5 ‰ offset between calcite and dolomite δ18O observed in virtually all Maieberg Formation sections.

The association between deglaciation and large changes in biological productivity (Canfield et al. 2007; Hoffman and Schrag 2000; Kennedy 1996; Peterson et al. 2004) render Zn isotopes as suitable biotracers (Isson et al. 2018). Zinc isotope ratios are significantly modified by biological (metabolic) activity (Dupont et al. 2006; Moynier et al. 2017 and references therein). By coupling Zn isotope ratios with the primary global productivity proxy, carbon isotope measurements, deeper insights to the cumulative biogeochemical processes during deglaciation could be gained (Hoffman et al. 2021).

Zinc isotopes have been previously reported in cap carbonates from the Nuccaleena Formation in the Flinders Ranges, South Australia (John et al. 2017; Kunzmann et al. 2013), and the Doushantuo Formation, China (Lv et al2018; Yan et al. 2019). These studies have generally observed a decrease in δ66Zn values, followed by an increase above initial values. In the Flinders Ranges this is from 0.47 to 0.07 ‰, followed by an increase to as high as 0.87 ‰. The first stage has been explained as a consequence of the influx of isotopically light Zn from weathering of isotopically light source rocks, while the second stage is attributed to the effects of the reactivation of export production (John et al. 2017; Kunzman et al. 2013). These trends have been suggested to reflect the consequence of global precipitation of carbonate enriched in heavier Zn isotopes, with a positive shift due to the removal of lighter Zn by sulfides (Yan et al. 2019). Of note is that these Zn isotope data are from platformal settings, which were flooded only after substantial ice-sheet melting had occurred.

In this work, we report isotopic data of Zn, C, and O from cap carbonates of the Otavi Group in Namibia and include multiple sample sites to assess the global or regional character of Zn isotopic patterns within the basal sections of the distal foreslope section, representing a much earlier stage of deglaciation than previous studies. This could date before the low-latitude ice sheet on the Otavi platform had disappeared (Hoffman and Macdonald 2010). Our emphasis is on the scale of variation in Zn isotope ratios in the immediate aftermath of the snowball event with δ13C as an excursion benchmark.

Materials and methods

Samples overview

For our study we focused on strata representing the transition from the late Cryogenian Ghaub Formation to the early Ediacaran Maieberg Formation of the Otatvi Group of northern Namibia. The samples were collected from near the base of the cap dolostone (Keilberg member (Mb) / Maieberg Formation (Fm)) representing by definition the onset of Ediacaran time (Knoll et al. 2007; Narbonne et al. 2012). These samples have been extensively described previously, and sample names are in accord with these studies (Gyollai 2014; Gyollai et al. 20142017).

The sample stratigraphy includes the Ghaub Formation, a fragmentarily preserved section of glacial deposits, rich in ice-rafted debris derived primarily from the underlying upper Ombaatjie Formation (Hoffmann and Prave 1996). It is separated from the overlying cap carbonates by a boundary, cm-scale reddish clay layer, developed from erosional surface conditions and tectonic overprinting. The Bethanis Member of the Ghaub Formation is a meters thick glacial debrite, fining upwards until its contact with the Maieberg Formation (Hoffman 2011). The Maieberg Formation is a ca. 400 m thick expanded transgressive-recessive deposition sequence, from which a multi-million year glacial duration was inferred (Hoffman et al. 1998). The Keilberg Member of the Maieberg Formation is a meters thick pale-grey to pale pink dolostone (Hoffman 2011). Samples representing the changeover between the uppermost Ghaub glacial diamictite and the basal post-glacial cap dolomite (Keilberg Mb / Maieberg Fm) were collected from proximal slope facies deposits (Bethanis Mb/Ghaub Fm) abutting the Otavi platform at outcrops near Fransfontein and near Narachaamspos and from outcrops representing inner Otavi platform facies deposits, located within the Khowarib Fold Belt near Ongongo as well as from an outcrop at the head of a valley extending from Entrance (near Devede village) towards the SW (Figs. 1, 2 and 3).

Fig. 1
figure 1

Geological map of the Otavi fold belt showing the Neoproterozoic carbonate dominated succession (Otavi group) and the foreslope-platform facies change. Fransfontain (sample C1) and Naraachamspos (C2, C3, C17) samples are from the foreslope. Ongongo (K2 and DC samples) and Entrance South Valley (C14) are from the Ghaub-Maieberg transition. (Map based on Hoffman, Snowballearth.org)

Fig. 2
figure 2

© 2020 Maxar Technologies

Location of our Ongongo Sample Locality (red Arrow, Profile 2). Annotated satellite image (from Hoffman et al. 2021, Figs. 112) looking obliquely NW toward Ongongo. Magenta lines indicate composite section through entire Otavi Group. Formations in stratigraphic order: Ob, Beesvlakte Fm; Od, Devede Fm; Ok, Okakuyu Fm; Ac, Chuos Fm; Ar, Rasthof Fm; Ag, Gruis Fm; Ab, Ombaatjie Fm; Tm, Maieberg Fm; Te, Elandshoek Fm; Th, Hüttenberg Fm; Ms. Sesfontein Fm. Rasthof Fm members: 1, rhythmite-turbidite; 2a, lobate stromatolite; 2b, roll-up microbialaminite; 2c, cryptic microbialite; 3, highstand grainstone. Google Earth: Image

Fig. 3
figure 3

Schematic representation of our sample site at the Ongongo locality. A tectonic hiatus exists between the Ombaatjie and Maieberg Formations overprinting a subglacial erosional surface of the Late Cryogenian period. Our measured sample K2 is a reddish tectonized clayely dolostone at the base of the cap carbonate. This sample was deposited above an erosional surface wherein Ghaub Fm diamictites of the Ombaatjie Fm were eroded and replaced with postglacial cap carbonate sediments. The measured drillcore (DC) samples cover 320 cm of the basal Keilberg cap carbonate deposits (DC5-49 = 140 cm above the base, DC5-50 = 240 cm above the base, DC5-54 = 320 cm above the base)

The samples from Naraachamspos and other foreslope locations represent the lithologically distinct basal turbidite interval (≤ 1.0 m) of the Keilberg Member (Hoffman et al. 2017). This is chronologically followed by a brief fall in relative sea level, at the start of a much larger relative sea level rise, in the Keilberg to middle Maieberg transgression, representing the “initial regression” observed in Greenland at the start of the last Quaternary deglaciation (Clark 1976). The initial shoaling at Naraachamspos, in turn, was interpreted as a sea-level fall resulting from the loss of gravitational attraction generated by the waning Otavi platform and Congo craton ice sheet on the regional ocean (Hoffman and Macdonald 2010). This would mean that the basal dolomite turbidites represent a very early stage of global deglaciation, following the collapse of the tropical sea glacier but with the ice sheets still remaining. At this point, the meltwater lid would be relatively shallow as it would derive from sea-glacier melting, allowing SBE brine to penetrate it given enhanced vertical mixing at the continental margins (Phillips 1970; Wunsch 1970), i.e. the foreslope locations.

Foreslope facies sample

The sample (C1C1) from the Otavi platform´s proximal slope facies is taken from the immediate transitional layers between glacial Ghaub deposits and post-glacial Keilberg cap carbonates. Further information is available in the electronic supplementary material (ESM).

Narachaamspos locality (Fig. S1 in the ESM)

The Fransfontein slope at Narachaamspos (“Pip´s rock”) exhibits stratified and finally ice-rafted debris flows of the Bethanis Member (upper Ghaub Fm), which is overlain with a sharp contact by stratified dolomicrite of the basal Keilberg Member. Here the C2 and C17 samples span the basal 20 cm of these capping dolomite strata, whereas sample C3 is a collection of samples from the direct contact of the two formations.

Fransfontein locality (Fig. S2)

Collected at outcrops NE from Fransfontein village, sample C1 comes from the top layer of a sequence composed of diamictite and well-bedded, fining upward dolostone choked throughout with outsized clasts of ice-rafted debris. These represent the upper Ghaub Formation (Bethanis Member). Paraconformably overlying the Bethanis Mb follows the cm-layered, micropeloidal Keilberg cap dolostone featuring giant wave ripple-structures at its base and transitioning into limestone rhythmite in the hanging wall. The Keilberg cap dolostone is characterized by its content of marly silt-mudstone.

Platform facies samples (Fig. S3)

Our samples from the Otavi Group´s inner platform facies region focus on Marinoan postglacial transition layers but are diachronous with respect to our slope facies samples. A) At the Ongongo locality, our sample profile 2 consists of 3 drill core-samples covering a 3.5 meter profile of the brown-gray dolomite from the distal foreslope of the Ghaub-Maieberg transition (Keilberg Mb/Maieberg Fm). Its basal contact (Samples K2) is a glacigenic erosional surface situated atop the Ombaatjie Formation. We collected these samples from the tectonized, clay boundary layers with respect to the underlying Ombaatjie Formation. Due to Pan-African folding and thrusting, this contact is tectonically overprinted and layers in the Ongongo area are vertically dipping. Concerning the subjacent upper Ombaatjie strata our sample profile 2 transgresses are an erosional surface incised into its basement due to glacial erosion that formed a Marinoan Trough cut into Ombaatjie deposits and afterwards rapid deglaciation caused removal of glacial recession deposits (e.g. Fig. 118 B in: Hoffman et al. 2021). B) Within the Khowarib fold belt our sample C14 (Fig. S4, Entrance South Valley) is collected from the thin red-green marly lutite marking the boundary between diamictite and cap dolomite, the latter being reddish stained by oxidized iron at its basal layer.

Analytical methods

Multi-isotopic analysis of cap carbonates (dolostone) samples dated to the Marinoan glaciation event (635 Myr) in Northern Namibia was performed for C, O, and Zn isotopic compositions. Samples were chosen based on previous X-ray diffraction, X-ray fluorescence, mineralogical analyses, Raman spectroscopy, and instrumental neutron activation analyses work (Gyollai 2014). The samples (Fig. 1) were collected seeking representative, pristine (unweathered) samples to best fit with the detailed stratigraphic setting, including localities and names from studies by Hoffman (20022005). As cap carbonates from this area are some of the most studied samples from the Neoproterozoic, they are appropriate for exploration through new isotope system measurements.

For the zinc isotopic measurements, powdered samples of 0.2 g were dissolved in cold 2N acetic acid to extract calcite (as per Pichat et al. 2003). Samples were alternatively sonicated for two hours and mechanically shaken to ensure full dissolution, until no CO2 was observed to be released. The acetic acid phase containing the carbonates was centrifuged and dried under a heat lamp, then redissolved in 1.5N HBr with the addition of 0.001 vol% H2O2. Separation of the Zn from its potentially interfering isobars and sample matrix was achieved through anion-exchange chromatography (Moynier et al. 2006, 2009). The total Zn yield of the chemistry is >99%, monitored through the use of United States Geological Survey standard rocks BCR-2 (Columbia river basalt, USA) and BHVO-2 (Hawaiian basalt, USA). The typical blank was less than 10 ng. Isotopic measurements of the isolated Zn were performed on a Thermo Scientific Neptune plus MC-ICP-MS following Paniello et al. (2012) and Chen et al. (2013). Replicate analyses of the same samples as well as standards define our external reproducibility (2σ deviation) of ±0.04 ‰ for δ66Zn, ±0.05 ‰ (δ67Zn), and ±0.07 ‰ (δ68Zn).

It is worth noting that in our study we use the same bulk carbonate dissolution of carbonates with acetic acid as was described in Pichat et al. (2003). It has been previously suggested that a “true” carbonate value can only be found through stepwise leaching of the sample (e.g. Liu et al. 2013; Liu et al. 2014; Bellefroid et al. 2018). While these studies likely find the correct isotopic value of the carbonate during the leaching, by the nature of leaching it is unclear what the actual value of the carbonate should be given the wide range of values that have been reported. Regardless, for the purposes of comparison between post-Snowball sample sites, bulk dissolution allows for a more robust comparison as this is the standard method for zinc isotope measurements in cap carbonates (John et al. 2017; Kunzmann et al. 2013; Yan et al. 2019). In this case one also analyzes other elements via wet chemistry, stepwise leaching may prove necessary, at which point a comprehensive methodological analysis of the behavior of Zn isotopes during stepwise leaching will become necessary.

For carbon and oxygen isotopic measurements, fine powdered samples were added to a quartz reaction vessel containing 100 % phosphoric acid (H3PO4) in the side arm. The sample tube was pumped overnight to remove trapped air in the acid. The acid was added to the evacuated carbonate powder at room temperature and allowed to react for 8h at 25±2 °C. The CO2 gas was cryogenically separated using two glass spiral traps (trap I at -65 °C to collect water and trap II at -196 °C to collect CO2) and non-condensable gas was removed from the system. Trap II was replaced with ethanol slush at -65 °C and CO2 was cryogenically transferred to the sample tube. After CO2 extraction, the sample tube was placed in a silicone bath at 150°c for 4 hours and CO2 released at this temperature collected as described above. This CO2 extraction at two different temperatures allowed separation of the Ca-rich (calcite) carbonate phase from the Fe-Mg-Mn-rich carbonate phase (Al-Aasm et al. 1990). Al-Aasm et al. (1990) have experimentally shown that this is one of the most effective ways to chemically isolate calcite from dolomite, siderite, and magnesite. At 25 degrees calcite is the only reactive species and is nearly 100% reacted, while all other carbonate species are removed at 50 degrees, as determined in their control pure component extractions. The present work employs their optimized extraction protocol. The C and O isotopic composition of purified CO2 was measured on a Finnigan MAT 253 Isotope Ratio Mass Spectrometer, as described in detail by Shaheen et al. (2010). Control samples of carbonate of known isotopic composition were routinely analyzed to insure isotopic precision and accuracy. For replicate mass spectrometric analysis 1 sigma standard deviations (SD) typically are 0.01‰ for δ13C and δ18O. Overall reproducibility for acid dissolution and gas chromatography steps for CO2 gas (n = 5) is ± 0.05 ‰ for δ13C and ± 0.1‰ for δ18O (Shaheen et al. 2015).

Results

The Zn, C, and O isotopic data from the cap carbonates are presented in Table S1 and are combined with previous data from Gyollai (2014). All δ13C and δ18O measurements are standardized to VPDB. The δ18O values vary between (n = 23) -3 and 12 ‰, with the calcite and low-Ca carbonates offset by 3 ‰ (Fig. 4) and δ13C (n=34) between 0 and -10 ‰, with most data points above -6 ‰. There is little correlation between the δ13C and δ18O values, though the Naraachamspos calcites show minor correlation. The data are compared to other measurements in the discussion section.

Fig. 4
figure 4

δ13C vs. δ18O: δ18O values lie primarily between -1‰ and 5‰. The calcites generally cluster around δ18O values of 3‰, while non-calcite (primarily dolomite) carbonates lie between -1 and 0‰. Values are given in standard per mil notation

Global temporal equivalents for the Maieberg δ13C excursion show depletion in 13C relative to typical seawater values of 0.0 ‰ (Halverson et al. 2010). The evidence of the excursion is also found in C isotopic data reaching an average nadir between -5 to -7 ‰, which agree with the ranges reported in previous studies of the Marinoan (Halverson et al. 2002; Hoffman et al. 1998; Hoffman and Maloof 1999; Kasemann et al. 2005; Kaufman et al. 1991).

Measurements of δ13C in this study were performed on chemically separated calcium-poor carbonates (primarily dolomite) and a calcium-rich phase (calcite). The goal was to determine the difference in the chemical record between the two minerals described in the experimental section. The extraction protocol removes calcites from other carbonates, following the experimental studies of Al Aasm et al. (1990). The δ13C value exhibits a general correlation between the chemically separated low-Ca carbonate and calcite. In the majority of the samples the non-calcite component (high temperature extract) has a concentration of released CO2 ~10 times that of calcite. At the sampling location of the C1 samples1 δ13C values remain constant while δ66Zn values vary widely. At C2 and C3, the δ66Zn values do not fluctuate, while the δ13C values reach the low of -6 ‰. At δ13C values of -3 ‰, δ66Zn undergoes greater variation. At C14 the δ66Zn value reaches its heaviest ratio, while the δ13C value is isotopically lightest.

The greatest range of variation of δ66Zn values coincide with δ13C values near zero (Fig. 5). These δ66Zn values exceed those typically measured in clastic sediments and igneous rocks. In turn, the δ13C excursions occur almost exclusively within the range of δ66Zn values of 0.3 and 0.5 ‰. The variation of Zn isotopic ratios beyond standard terrestrial ranges occurs in association with δ18O values of 2 to 5 ‰. In turn, δ18O variations occur primarily within the relatively narrow range of δ66Zn values 0.4 ± 0.1 ‰.

Fig. 5
figure 5

δ66Zn (WR) vs δ13C (both phases) in per mil: In samples from both the foreslope and platform, the zinc isotopes are decoupled from the carbon (R2 of 0.0931 for foreslope, 0.1575 for platform). The biggest variation in δ66Zn occurs in samples which do not show the negative carbon isotope values associated with post snowball samples, while the largest variations in δ13C show δ66Zn values which fall towards the upper range of most modern terrestrial samples (here titled Earth crust, spanning the range of modern clastic sediments and igneous samples), although almost all samples are far below modern carbonates (dark grey). In blue, we highlight the mixing of the C and Zn isotopic reservoirs, showing the two as largely decoupled

Discussion

Overview

The SBE ended with high pCO2, giving a distinct isotopic signature in cap carbonates. The initial post-glacial weathering and nutrient fluxing was likely dominated by highly reactive glacial rock powder (loess) exposed by the receding ice sheets, on account of its high surface area to volume and unweathered composition. Intense chemical weathering on the continents led to a high flux of nutrients and high primary productivity. However, given the low primary productivity during the Precambrian, C isotopes cannot directly reflect primary productivity while Zn potentially can. While previous studies discussed this possibility by reporting Zn isotopic data at low spatial resolution (ca. 1 m between samples), here we analyze at far higher (1 cm) resolution to better constrain the small-scale variations (Fig. 6).

Fig. 6
figure 6

Stratigraphic Profile. Section of the Naraachamspos distal foreslope. While the profile is the first 17cm, the member extends for a further 20m. The δ18O (VPDB) broadly increases through time. The extremely negative δ13C values are clustered towards the base of the section, in the bottom few centimeters of the profile. These are the most likely section to be altered through alteration at the dolomite contact. Finally, the δ66Zn varies wildly through time, likely a reflection of ocean mixing times (between 100s to 1000s of years for modern oceans)

As most modern terrestrial inorganic zinc has a relatively narrow range of δ66Zn values (Fig. 5), the wide range observed in the dataset suggest different Zn sources control the isotopic signatures. These sources provide estimates to contributions which drive the secular change previously reported for Zn isotopes. We will discuss these processes alongside the variations in δ13C and δ66Zn values and the decoupling of the two isotope systems.

Oxygen, carbon, and phases of carbonate

δ18O and δ13C ranges (Fig. 4) are consistent with the composite carbonates data reported by Jaffres et al. (2007). The modest variability in δ18O compared to δ13C suggest that both calcitic and non-calcitic carbonates have undergone mild diagenesis (Bao et al. 2008; Zhou et al. 2010). Furthermore, little evidence for diagenesis based on trace element and petrographic analysis have been observed in these samples (Gyollai 2014). Our δ18O data in calcite are similar to those of non-calcite carbonate, but do not correlate with the δ13C variations. This decoupling may indicate a partial reset in the oxygen record given that secondary alteration processes conventionally smooth out isotopic structure (Swart 2015). Of note is the particularly high values of δ18O (up to 12) reached in some samples, as such high oxygen isotopic values could indicate evaporitic conditions. However, this is by no means a “smoking gun,” simply the most likely explanation for the magnitude of effect. Macroscopic evidence suggests that calcite in the Keilberg Member cap dolomite is a late-stage (possibly post-tectonic) fracture- and void-filling sparry cement, the origin of which is unrelated to Neoproterozoic seawater other than being buffered by Neoproterozoic carbonate rock along its possibly meteoric fluid pathway. Accordingly, the calcite in the cap dolomite tends to be more 13C depleted than the dolomite because the pore-fluid from which it precipitated was buffered by the thick overlying middle Maieberg Formation limestones, which are consistently 2−3 permil more 13C depleted than the cap dolomite. However, in these older samples, there is little evidence of those same events.

Previous studies focused on bulk carbonates through micro-drilling sawed surfaces to avoid macroscopic late-stage fracture and void-filling calcite, though this would not exclude microscopic calcite. In contrast, our study employs chemical separation and finds that the calcite tends to exhibit a more negative δ13C value than the non-calcite. The difference (average 0.5 ‰, n=34) between the calcite and dolomite dominated non-calcite could derive from different sourcing of materials rather than a depositional selective isotopic effect where one would expect a transition between the carbon isotopes rather than clusters. While calcite easily forms from dissolved CO2, dolomite may form in areas richer in organic matter (Bramlette 1946, Kelts and McKenzie 1982, Swart 1993). Therefore, the different regions may capture redox/pH changes, however one cannot conclusively state that this is the inherent mechanism for dolomite formation during diagenesis. In the post-glacial depositional environments, this means calcite incorporates atmosphere CO2 more directly, while non-calcite carbonates more extensively sampled the diagenetic environment in their formation (Swart 2015). This could explain why the non-calcite carbonates are isotopically lighter than the "typical" carbon signature following deglaciation (Hoffman et al. 2017). It is likely that more than one fractionation factor (e.g. pH, temperature, and molality of solute species – especially electrolytes and organics) is operative, and interlinked with environmental changes. These effects would include single stage fractionation factor variations, chemical reaction mechanisms, and position of isotopic chemical equilibria. Interestingly, at the isotopically lightest values the calcite is lighter than the non-calcites. This reversal effect could be the result of a change in the composition of the source region or of the chemical formation mechanism of the non-calcite carbonates, which would change the concomitant relevant isotopic fractionation factors (Swart 2015). Determination of the specific compositional/chemical differences could be found through trace element proxies or alternate isotopic results (Abanda and Hannigan 2006). Our measured range of δ13C values (between 0 and -10 ‰) is similar to that found in bulk carbonates by Hoffman et al. (2007). These variations have many potential sources, as previously discussed (Bechstädt et al. 2018). These observations thus lead to a next step in quantifying the relative contributions of different proportions of carbonates to chemical/depositional and source strength variations over time.

Δ13Ccalcite – non-calcite (the difference between the δ13C of calcite and non-calcite) is variable but typically negative. Differences in δ13C between coexisting dolomite and calcite indicates an absence of homogenizing alteration that post-dates the calcite. If the calcite is a late-stage (burial or post-burial) product, the earlier dolomite might still be altered from its composition on the seafloor (e.g. Ahm et al. 2019; Hoffman and Lamothe 2019). Similarly, δ18O values show offset between calcite and non-calcite. Secondary or random scattering could account for a diagenetic reset affecting only one phase. As this is not observed, we rule this process out. This is probable, as Naraachamspos is paleobathymetrically below the level of intense fluid-buffered diagenesis (Hoffman and Lamothe 2019). The extremely high Sr/Ca ratios (as reported in the ESM) can also indicate a primary signature, although this is not a reliable metric.

Comparison of zinc and carbon isotope ratios

We observe little correlation between the δ66Zn and δ13C values amongst the samples (Fig. 5). The δ13C values are low, consistent with models invoking a greater role of biological organics in the upper ocean, as well as upwelling of deep water. Upwelling deep water would be enriched in nutrients, leading to algal blooms in the shallower parts of the ocean, where the samples were deposited. The observed δ66Zn falls primarily between 0.2 ‰ and 0.6 ‰, with some values reaching as low as -0.3 ‰. These ranges could be the result of crustal contributions or upwelled deep water, with the negative excursions the result of enhanced biological activity. Uniformity in δ13C could reflect equilibration of meltwater-lid dissolved inorganic carbon with the larger atmospheric CO2 reservoir, given that carbonates in the turbidites formed upslope from their site of final deposition. The massive carbon pool of this dissolved inorganic material would also potentially mask smaller changes in δ13C values. δ13C values of -0.5 ± 0.5 ‰ fall within the δ66Zn range of -0.30 to 0.45 ‰. δ13C values below -1 ‰ uniformly track with a δ66Zn of 0.45 ‰. The absence of correlation between the (possible) bio-signatures of primary productivity in zinc and carbon isotope systematics demonstrates that biological processes may not exclusively account for the isotopic values. The lack of correlation also demonstrates the value of simultaneous isotopic measurements, as it allows for elimination of potential causes of isotope fractionation.

Differences in δ66Zn between the calcite and non-calcite carbonates (Fig. 5) would likely be due to their different structures and chemical compositions (i.e. Mg and Ca). The larger ionic radius of Ca compared to Mg (and Zn) affects the bonding environment of Zn and therefore would result in distinct isotopic compositions of Zn between calcite and dolomite. However, the coupled carbon (calcite and non-calcite) isotope variations along the profiles indicates that there must also be some source compositional variations as C mineral structure changes would decouple between phases. This argument is qualitative, pointing to the need for extended measurements of coupled isotopic systems to parse out ambient changes and their role in weathering versus redox conditions. This would allow for a comprehensive model that folds in future experimental measurements of biological fractionation factors for more species and over a greater range of pH and temperature. More extended trace element and multi-isotope measurements in differing depositional environments are important. Ultimately, the roles of biology, climate, and geochemical weathering will be better understood. Nevertheless, our data, by coupling for the first time Zn and C isotopes in calcite and non-calcite minerals, suggest that the change of mineral structure is not the only control on carbonate Zn isotopic composition. Instead, variations in the Zn source must also exert a strong control on the final observed isotopic composition.

Zinc isotopes

An important facet in interpreting the Zn isotopic results is to determine the different Zn sources. Distinct Zn isotope compositions are associated with changes in Zn speciation, metabolic processes, hydrothermal activity, and geochemical weathering. Changes in the speciation of Zn (particularly when sulfides are involved), affect the Zn isotopic composition of the available pool of Zn in the seawater, and therefore could affect the isotopic composition of the carbonates (Fujii et al. 2011). The proliferation of life into a newly sunlit, nutrient-rich, and rapidly warming ocean following Snowball deglaciation strongly suggests biology as a potential source. The long period of liquid ocean isolation with exclusively hydrothermal inputs would suggest that hydrothermal values should appear in the record.

The speciation of zinc within the ocean is primarily pH controlled, which can affect its residence time and concentrations (Bruland et al. 1994; Fujii et al. 2011, 2014). Following global deglaciation, the oceanic pH may have fluctuated due to high variations in dissolved CO2 and seawater composition, as well as delivery of acidic and alkaline materials through pathways such as weathering. The pH of standard ocean water during the Cryogenian is modeled as 8.1, the same as present, though in the post-glacial periods the pH was suggested to have dropped as low as 7.0 (Kasemann et al. 2005) or lower (Ohnemueller et al. 2014). At both pH 7 and 7.8, as in the modern ocean, the dominant species of Zn in the modern ocean is Zn2+ (Zirino and Yamato 1972). However, Zinc hydroxide (Zn(OH)2) overtakes Zn2+ as the dominant Zn species at pH 7.8 and Zn will have a far shorter residence time in more alkaline conditions. The fractionation factor for calcitesolid-liquid is pH dependent with an 0.6 ‰ range of fractionation between pH of 8.5 and 6.1 (Mavromatis et al. 2019). This fractionation factor remains constant, irrespective of ambient pH, rendering Zn isotopes an ideal pH proxy. Furthermore, the Zn isotopic fractionation process is highly sensitive to ionic strength, and a near 0.4 ‰ variation in Zn isotope ratios may be achieved by variation from high to low ionic strength (Dong and Wasylenski 2016). The presence of anoxic iron rich deep ocean environments along with enhanced hydrothermal zinc concentrations would provide high levels of oceanic zinc, while post-glaciation ice melt or alkaline detrital inputs would provide zinc deficits. Overall, changes in the redox states of the ocean would affect biochemical access to zinc as a consequence of the change in its chemical speciation and bioavailability (John et al. 2007; Scott et al. 2013; Little et al. 2014).

The large detrital content of the post-glacial upper ocean lifts the restraints on limiting nutrients (such as bioavailable Zn), creating an ideal environment for reproduction and global primary productivity expansion (Brand et al. 1983). The δ66Zn of the modern ocean (0.5‰) is higher than crustal values (~0.3‰), in part the consequence of biological activity (Bermin et al. 2006; Little et al. 2014; Vance et al. 2016). δ66Zn values vary with oceanic depth, the consequence of biological processes (Bermin et al. 2006; Conway and John 2014; John et al. 2008; Lohan et al. 2002; Zhao et al. 2014) or reversible scavenging of Zn onto organic matter by ion-surface electrostatic attachment (Conway and John 2015; John and Conway 2014; Weber et al. 2018). Zinc isotope composition follows the depth profile of concentration, with surface values as high as 0.8‰, decreasing to -0.1‰ at 80 m depth before increasing back towards a more globally homogeneous deep-water δ66Zn value of 0.5‰. Isotopically light Zn is utilized in the productive zone, while dead organic matter is remineralized at ca. 80 m, and stable values are reached at depth. This light Zn remains in biologically active use by other organisms, potentially producing two distinct isotopic zinc reservoirs as heavier Zn is removed from the biosphere (John et al. 2007; John and Conway 2014; Moynier et al. 2017). This could account for our observed variations in δ66Zn. The present limit is that there are finite culturing experiments (species wise) and the precise zinc concentration, pH, and zinc chemical state are not sufficiently constrained.

Isotopically heavy Zn observed in our carbonates therefore also may be explained as the combination of biological uptake and reversible scavenging of Zn onto sinking particles (Gelabert et al. 2006; John et al. 2007; John and Conway 2014; Weber et al. 2018). As shown by John et al (2007) the adsorption of heavy zinc onto cell surfaces undergoes a desorption reabsorption process (John and Conway 2014). The unstable climate post-snowball could lead to largescale biosphere die-offs, allowing for different reservoirs (i.e. hot or cold) to be captured in the carbonate formation process.

Many of our observed δ66Zn values beyond the range of typical crustal values may relate to the addition of hydrothermal Zn, as hydrothermal Zn can bear some of the lightest δ66Zn values (-0.51 to 0.01 ‰) (Lemaitre et al. 2020). Although the influence of hydrothermal Zn is relatively small compared to terrestrial inputs in the modern ocean (Little et al. 2014), in the absence of continent – ocean interactions during the Snowball phase, hydrothermally sourced Zn will become the dominant input of Zn to the ocean. Support for SBE brine-meltwater mixing from the top of Marinoan cap dolomites exists in the form of sea-floor barite, in which Ba-isotopes indicate a globally homogenous source indistinguishable from modern pelagic Ba (i.e. the well-mixed residual SBE brine), and sulfate Δ17O carrying the SBE atmospheric O2 signature indicating that the sulfate was delivered in post-SBE run-off from oxidative weathering of sulfides in loess, tills, and rocks (Crockford et al. 2018). The expected range of δ66Zn in the upwelled water would fall below crustal values (0.3 ‰). This is one of the most likely sources of the negative δ66Zn isotope signatures. Variations may be attributable to cycling of the hydrothermally dominated deep waters and deglaciation derived water reservoirs (Yang et al. 2017). However, it is difficult to assert whether these values are exclusively hydrothermal and biological fractionation processes may also be significant.

Other sources of variations of δ66Zn within our own samples beyond 0.3 ± 0.1 ‰ could be the result of zinc-containing sediments being deposited in the ocean by continental erosion during the wet, warm period following deglaciation. This is unlikely, however, as the southwest cape of Congo craton, if not the entire craton, was blanketed by carbonate rock before, during, and after the Cryogenian SBEs, and thus would not be a source for non-carbonate terrigenous Zn. Furthermore, our methods exclusively extract zinc from the carbonates, so potential heterogeneities introduced directly via clastic non-carbonate components are unlikely. The variations of Zn isotopes related to mixed signals from overprinting, the carbonates, and the clastic non-carbonates are best determined using step-leaching (e.g. Liu et al. 2013; Liu et al. 2014; Bellefroid et al. 2018). Even if these components were dissolved to the extent that the Zn was incorporated into the carbonate phase, which is chemically very unlikely due to the electron volume of Zn, the seawater record would overprint it given the differences in reservoir sizes. As such, we suggest that the primary δ66Zn sources were likely upwelled hydrothermally dominated deep-ocean water in the aftermath of glaciation, with variations derived from biological activity. As oceanic Zn mixing and residence times were roughly the same (ca. 51 kyr) (Yang et al. 2017), regional Zn heterogeneity is unsurprising.

Relevance to previous studies

Comparison of our results with those of previous studies (John et al. 2017; Kunzmann et al. 2013; Lv et al. 2018; Yan et al. 2019) allows for an assessment of the global or regional character of the Zn isotope evidence. However, our samples are sourced from the distal foreslope, while these other studies have analyzed platform settings. As such, direct comparison is not possible. In both cases, the initial δ66Zn may be representative of continental, hydrothermal, and deep-water inputs. Kunzmann et al. (2013) and John et al. (2017) results from the Nuccaleena section of the Flinders Ranges trend towards lighter values than typical carbonates, which we also observed in some of the transition layers. These δ66Zn values are in the range of mantle values (0.2 to 0.4 ‰), and those of hydrothermally dominated systems (e.g. John et al. 2008), as would be expected following resumption of oceanic mixing following deglaciation given the dominance of volcanic inputs during the Snowball period. While the general trend is mirrored amongst the samples in the Jiulongwan Section of the Doushantuo Formation (Yan et al. 2019), rather than a descent from ca. 0.4 to 0.2 ‰, these samples fall from 0.6 to 0.4 ‰. Thus, there is a similar magnitude of change between the Jiulongwan and Nuccaleena sections, indicative of a global oceanic trend. It is noteworthy that the range of variation in δ66Zn, which we observe in the 17 cm span of one profile from the basal turbidite interval in a distal foreslope setting (Fig. 6) is mirrored in the meters long profiles from full-thickness cap dolomites in platform settings shown in Fig. 5 of Yan et al. (2019) (containing data from John et al. 2017; Kunzmann et al. 2013; Lv et al. 2018). It is worth noting that the stratigraphic span of our data (ca. 20 cm) falls between any two data points of the previous studies (50 cm – 1 m between samples). As such, determining the larger and true extent of Zn isotopic trends can be improved with smaller distances between samples, to avoid the risk of sample bias creating a potential false signature.

While it is tempting to arrive at a lifetime and magnitude of the zinc reservoirs, there are significant restrictions. As discussed, the light isotope effect from cellular uptake is highly dependent on source Zn concentrations which, unfortunately, is not known. The same is true for temperature dependence of the fractionation process. The chemistry of surface adsorption, desorption, and immediate re-adsorption is known to produce heavy isotope enrichment, but this has been determined for one species at one temperature (John et al. 2007). The molecular level details are not known, and thus we cannot currently address the quantitative role of the oceanic speciation of zinc (elemental vs organic ligand complexed). The hydrothermal process is restricted in quantification of the process as it is a temperature dependent fractionation process and we do not know the global extent or temperature of these sources following the Marinoan. The actual chemistry of the subsurface process also requires a more detailed mechanistic description.

It is plausible that the biospheric signal overrides the hydrothermal system, as has been suggested based on the increasing dominance of eukaryotes (Isson et al. 2018), but given the restrictions listed above, it is not possible to generalize this. Simultaneous sulfur and boron isotopes as well as rare earth elements and iron isotopes (e.g. Busigny et al. 2018) may provide better evidence of hydrothermal waters’ modulation of the carbonate isotopic signature as sulfur captures both source and redox state parameters with the multiplicity of available sulfur valence states (Thiemens 2006; Thiemens and Lin 2019), and serve as an expansion on the work of Ohnemueller et al. (2014). Boron records the pH, capturing another parameter for the post-glacial oceans. The precipitation of sphalerite (Zn sulfide) could cause hydrothermal sources of Zn to trend towards becoming isotopically heavy. This, in turn, leads to a high potential for measuring sulfur and zinc isotopes together in future experiments. This would be especially helpful for redox consideration as the four isotope system of sulfur and its possession of multiple valence states (plus, minus, and zero) render it a sensitive measure of chemical mechanisms and redox conditions, including sedimentation and biological interactions.

An important focal point of the work has been to address the issue of local versus non local sources of environmental change in this time period. The Zn isotopic variations are measured in several regions of Namibia, suggesting that these signals may record localized events (Site 1, C14, C15), without ruling out global events with impact of local geochemical signatures. These isotopic variations are more likely to be the result of biological fractionation of Zn isotopes (Isson et al. 2018), with isotope ratios below the continental average being non-uniquely characteristic of biological processes. These data further support the role of local biological processes, as one would not expect a homogeneous global recovery of the biosphere, and variations in a smaller local reservoir may produce an amplified signal due to mass conservation requirement for a smaller environment.

Conclusions and future work

Multi-site, multi-isotope analyses of the Marinoan cap carbonates from Namibia reveal multiple contributing Zn isotopic processes in the aftermath of the SBE glaciations. Non-crustal zinc isotope ratios (δ66Zn values beyond the range of 0.3 to 0.5 ‰) are unlikely to be the consequence of detrital inputs, leaving biological activity as the most likely explanation for these values. However, mixing reservoirs, such as upwelled hydrothermally dominated deepwater, could be another potential source. Local variations in δ66Zn values are found, with the isotopically light excursions best explained as being derived from local resumption of biological activity, likely in the form of algal blooms.

We find that the most negative δ13C values, from -4 ‰ and lower are associated with δ66Zn values of approximately 0.45 ‰. The primary source of the carbon signature is likely the consequence of enhanced primary productivity during post-glacial melting. However, this leaves variations in the carbon signature unaccounted for, which require more detailed analysis.

The key finding is that we observe the same magnitude of variation in Zn isotope compositions (ca. 0.8 ‰ variation in δ66Zn) for an extremely short timescale following deglaciation as seen in other profiles of the cap carbonates (ca. 0.6 ‰). The full-thickness cap dolomites in South Australia (0.2 – 0.8 ‰) and South China (0.4 – 1.0 ‰) have different ranges of variation as this study (-0.3 to 0.5 ‰), though the ranges are demonstrably heavier. The lack of homogeneity amongst the δ66Zn values in the samples from this study is the likely consequence of regionalized export production, as Zn taps multiple sources. A key next step would be to obtain and analyze full-thickness sample sets spanning the Keilberg Mb cap dolomite on the distal foreslope (6-10 m) and inner platform (ca. 20 m). This would allow comparison of δ66Zn with records from other paleocontinents. Most critically, though, it would construct a more thorough image of an interval of time and paleoenvironment not recognized elsewhere.