The 1.85 Ga Belomorian Belt, Karelia, Russia, hosts ultralow δ18O and δD (as low as −27.3‰ and −235‰ standard mean ocean water [SMOW], respectively), high-Al gneisses and amphibolites that we attribute to the Paleoproterozoic “Slushball Earth” glaciation. They now occur in at least 11 localities spanning 450 km. To constrain distribution of 18O-depleted rocks, we performed detailed field mapping in Khitostrov, where δ18O values are the lowest. Using 430 new and previously published laser fluorination isotope analyses, we show that the elongated, concentrically zoned area of δ18O depletion is greater than 6 × 2 km in areal extent, ∼10 times larger than previously thought. Relationships between δ17O versus δ18O strictly adhere to the equilibrium terrestrial mass-dependent fractionation with a slope of 0.527. We also report the results of ion microprobe U-Pb geochronology of zircons coupled with co-registered oxygen isotope spot analyses for mafic intrusions and host gneisses in six localities. The 2.9–2.7 Ga gneiss zircon cores are normal in δ18O (5‰–7‰). They show truncated oscillatory cathodoluminescence (CL) patterns and rounded shape indicative of original igneous crystallization with subsequent detrital overprinting. A younger 2.6–2.55 Ga metamorphic zircon domain with normal δ18O, low Th/U, dark cathodoluminescence, and also with rounded crystal morphology is commonly preserved. Cores are surrounded by ubiquitous rims highly depleted in δ18O (re-)crystallized with Svecofennian (1.85–1.89 Ga) ages. Rims are interpreted as metamorphic due to bright and uniform CL and Th/U <0.05. Mafic intrusions preserve few igneous zircon crystals between ca. 2.23 and 2.4 Ga in age, but neoblastic zircon in these intrusions originated mostly during 1.85 Ga Svecofennian metamorphism. The δ18O-age relationship for metamorphic rims in zircon and corundum grains suggests that δ18O values of fluids were subtly increasing with time during metamorphism. Large metamorphic corundum grains have ∼3‰ intracrystalline δ18O isotope zonation from –24 to –21‰, which likely developed during interaction with metamorphic fluids. The Zr-in-rutile geothermometer temperatures are in the range of 760 to 720 °C, in accordance with mineral assemblages and amphibolite metamorphic grade. High and irregular rare-earth element (REE) abundance in cores and rims of many zircons correlates with high phosphorus content and is explained by nanometer-scale xenotime and monazite inclusions, likely in metamict zones during 1.85 Ga Svecofennian metamorphism. A survey of oxygen isotopes in ultrahigh-pressure diamond and coesite-bearing metamorphic terrains around the world reveals normal to high-δ18O values, suggesting that the low δ18O in metamorphic rocks of Dabie Shan, Kokchetav, and in Karelia, are genetically unrelated to metamorphism. We discuss alternative ways to achieve extreme δ18O depletion by kinetic, Rayleigh, and thermal diffusion processes, and by metamorphism. We prefer an interpretation where the low-δ18O and high-Al signature of the rocks predates metamorphism, and is caused by shallow hydrothermal alteration and partial dissolution of the protolith surrounding shallow mafic intrusions by glacial meltwaters during pan-global Paleoproterozoic “Slushball Earth” glaciations between ca. 2.4 and ca. 2.23 Ga.

There are now 11 known localities spanning over 450 km across the Belomorian Belt, which formed in the Late Archean–Early Proterozoic, where we have observed remarkably 18O depleted rocks (as low as −27.3‰ standard mean ocean water [SMOW], Fig. 1)1. Ultralow-δ18O rocks crop out over many tens to hundreds of meters, and comprise chiefly four lithologies: Al-rich paragneisses, 2.4 Ga high-Mg gabbro-noritic intrusions, amphibolites at the contact between the two, and high-Fe intrusions with a tentative age of 2.3–2.1 Ga. The δ18O depletions often display a “bulls-eye” concentric pattern with progressively greater 18O depletion in the proximity of the intrusions (Fig. 2), which are similar to trends observed around shallow intrusions in modern meteoric-hydrothermal systems (Taylor, 1974). A feature of aluminoferrous Karelian rocks (Fig. 3) is that the maximum level of 18O depletion is also characterized by the highest concentration of Al, insoluble major and trace elements in the paragneisses, most notably in their high-Al/Si, Ti, and Zr concentrations. These features can be explained by residuum enrichment during partial dissolution of silicate minerals in hydrothermal process at large water/rock ratios (Bindeman and Serebryakov, 2011). It is also undoubtedly important to resolve the influence of ca. 1.85–1.89 Ga metamorphic processes, and in particular the effects of synmetamorphic fluid flow, on the protolith of these rocks. Preservation of diverse δ18O protolithic values during and after high-grade metamorphism has been described in other terranes (e.g., Valley and O’Neil, 1984; Fu et al., 2012), but synmetamorphic origin of ultralow δ18O Karelia rocks is advocated by some researchers (e.g., Ustinov et al., 2008). Karelian rocks represent the current world record in the level of 18O depletion. The only known terrestrial oxygen reservoir that could conceivably cause rocks and minerals to attain such low δ18O values is glacial meteoric water (e.g., Bindeman, 2011). We thus envision high-temperature water-rock interaction in a subglacial rift zone where ca. 2.4 Ga and perhaps younger Paleoproterozoic mafic intrusions have caused the depletion in the pre–1.85 Ga metamorphic protoliths. Because Karelia was located at near equatorial latitudes during most of the Paleoproterozoic (Evans and Pisarevsky, 2008), these ultralow-δ18O, high-Al paragneisses were interpreted to represent the first direct evidence for pan-global “Slushball Earth” glaciation (Bindeman et al., 2010). Dating of the intrusions can also serve as a novel method to constrain the individual glaciations, or the total duration of a pan-global freeze. The current geologic and geochronologic data suggest three to four individual glaciations between 2.5 and 2.26 Ga (Young, 2004; Hoffman, 2009; Bekker, 2011; Hoffman, 2013; Rasmussen et al., 2013). Which of the three events caused the appearance of atmospheric oxygen and disappearance of mass independent sulfur-isotope fractionation is a matter of current debate.

Here, we expand our previous isotopic mapping of the Belomorian Belt and present 430 new analyses of individual minerals, mostly garnet, by laser fluorination (Fig. 2; see Table A1 in the Supplemental File2). We also report new analyses for the southernmost low-δ18O locality, Shueretskoye (Fig. 1B), which expands the zone of known low-δ18O localities to 450 km across the Belomorian Belt. We also present 216 high-spatial resolution geochronologic and oxygen isotope spot analyses from 11 rock units at six localities, with auxiliary trace-element analyses of zircon. Zircon is unique for these mineral assemblages in that it retains normal-δ18O cores, which survived hydrothermal alteration and subsequent metamorphism. Metamorphic recrystallization during ca. 1.85 Ga Svecofennian metamorphism formed some neoblastic low-δ18O zircons in equilibrium with other metamorphic minerals but mostly triggered the epitaxial crystallization of low-δ18O rims onto normal-δ18O cores. In order to determine the age of depletion, we use the age of the youngest detrital zircon core in the paragneisses to constrain the age of last deposition prior to hydrothermal alteration and metamorphism. We further have determined the ages and δ18O values of zircons in mafic intrusions that are present in the vicinity of the ultralow-δ18O rock halo. In an attempt to constrain hydrothermal alteration in nonmetamorphic equivalents of the gabbronorite and gneiss complex, we investigated Sumian and Sariolian ca. 2.4–2.2 Ga Karelian supracrustal rocks for their δ18O composition: basaltic pillow rims, amygdaloidal vugs, varves from lacustrine periglacial lakes, as well as Tertiary basalts and their alteration products (vugs) from Antarctica for comparison purposes.

Sampling and Laser Fluorination Analyses

For isotopic mapping of the Belomorian Belt, hand specimens were processed in the field to extract 1–2 mg of garnet and other major phases. The separates were then transferred to the University of Oregon Stable Isotope Laboratory. For the majority of mapping samples, a single crystal of garnet was analyzed (Table A1 in the Supplemental File [see footnote 2]). Bulk oxygen isotope analyses of 0.5–2 mg aliquots of garnet and/or mineral separates of plagioclase, ruby corundum, kyanite, biotite, amphibole, zircon, monazite, and rutile were conducted using laser fluorination (e.g., Bindeman, 2008). Samples were heated with a 35W NewWave Research infrared laser in the presence of purified BrF5 reagent to liberate oxygen. The O2 gas generated in the laser chamber was purified through a series of cryogenic traps held at liquid nitrogen temperature and a mercury diffusion pump to remove traces of fluorine-bearing waste gases. The oxygen was converted to CO2 gas in a small, heated platinum-graphite converter, the yields were measured, and then the CO2 gas was analyzed using a Thermo Scientific MAT 253 mass spectrometer in a dual inlet mode. Four to seven garnet standard aliquots (UOG, δ18O = 6.52‰ and GMG, δ18O = 5.75‰) were analyzed together with the unknowns during each of seven analytical sessions. Measurements of unknowns were adjusted to correct for day-to-day variability, and precision of the standards was typically <0.1‰ (1 standard deviation).

Because Karelian samples span the largest yet measured terrestrial range in δ18O, we additionally performed three isotope (16, 17, 18) oxygen isotope measurements including 17O of several Khitostrov samples spanning a 37‰ range in δ18O, by measuring the isotopes in O2 gas directly without conversion to CO2, to check for mass independent isotope fractionation. The procedure includes twice collecting the gas on chilled 13 Å molecular sieve to prevent 14NF+ (mass/charge = 33) contamination. This interference was monitored by scanning for 14NF2+ (mass/charge = 52), which was below detection, suggesting the absence of contaminants in the O2 gas.

Ion Microprobe U-Pb Dating, Oxygen Isotope, and Trace-Element Analysis

Zircons were extracted from crushed samples using standard density separation procedures involving heavy liquids and magnetic separation. Extracted zircons were mounted in the center of a round (2.54 cm in diameter) epoxy mount along with standards and polished to expose crystal interiors. The crystals were imaged by backscatter electron (BSE) and cathodoluminescence (CL) methods prior to analysis. Ion microprobe U-Pb dating of zircons (Table 1) was performed at University of California, Los Angeles (UCLA) using the CAMECA IMS 1270 large magnet-radius ion microprobe, applying routine instrumental and calibration procedures based on zircon standard AS3 (Schmitt et al., 2003; Bindeman et al., 2010). A subset of zircons was dated using the Stanford–U.S. Geological Survey (USGS) sensitive high-resolution ion microprobe reverse-geometry (SHRIMP-RG) ion microprobe at Stanford University using a U-Pb dating protocol that is calibrated to zircon standard R33 (419 Ma, Black et al., 2004) and includes trace-element analysis. Sample AB3513 was dated by SHRIMPII in VSEGEI, St. Petersburg, Russia, and used TEMORA (Black et al., 2004) zircon as age standard.

Oxygen isotope analyses of zircons were performed at UCLA using methods described in Trail et al. (2007). After repolishing the mount to level the crystal topography and remove all traces of oxygen implanted during the U-Pb dating analysis, a ∼3 nA Cs+ primary beam at 25 µm spot diameter was targeted directly onto the same crystal domains used for dating. Beam size and repolishing likely resulted in core/rim overlap in some analyses, especially for small zircons such as in sample K5 (Bindeman et al., 2010). Instrumental fractionation was calibrated using bracketing and interspersed analyses of zircon standards mounted together with the unknowns. Standard values were as follows: AS3, δ18O = 5.31‰ (Trail et al., 2007), KIM5, δ18O = 5.09‰ (Valley, 2003), and TEMORA, δ18O = 8.20‰ (Valley, 2003). Instrumental mass fractionation factors (IMFs) varied between ∼1 and 3‰, with shifts occurring after exchanging sample mounts, but δ18O drift was absent for runs of individual mounts. Uncertainties for individual spot analyses (Table 1) are based on the external reproducibility of the standards in each analytical session, and average ∼0.2‰. In addition, we analyzed zircon rims by pressing euhedral grains into indium metal along with standards. Because zircon prismatic growth domains are thus laterally extensive perpendicular to the direction of ion beam penetration, the composition of the outermost rims of zircon can thus be analyzed to a depth of ∼1 µm. In another experiment, we applied a ∼0.5 µm Ga+ beam to detect oxygen isotope variations in an ∼25 µm lateral profile. This technique (a novel method, developed for this work by AKS using the UCLA CAMECA 1270 ion microprobe) used pre-implantation of the analysis area with Cs+ to enhance secondary ion yields, and simultaneous Faraday cup measurements of 16O with ion counting of 18O using a Hamamatsu Mark III electron multiplier. The instrumental fractionation value determined on 91500 zircon was ∼−40‰, with a spot-to-spot reproducibility of ∼1‰ (1 standard deviation [s.d.]). Analyses of zircons are presented in Table 1.

Analysis of oxygen isotopes within a large corundum crystal from sample K1 (Fig. 4) was conducted with a CAMECA 7fGEO ion microprobe at California Institute of Technology (Caltech). Given the large crystal size (∼1–2 cm), synthetic corundum (−6.37‰ SMOW, determined by laser fluorination at the University of Oregon) was mounted around the margins and Instrumental Mass Fractionation (IMF) corrections were applied using standards in close proximity to the unknown. We estimate the overall error associated with these corrections to be ∼0.5–1‰ (1 s.d.). The range of observed core-to-rim variability was independently confirmed by laser fluorination analysis (Fig. 4).

Analyses of trace-element concentrations in zircons (Table A4 in the Supplemental File [see footnote 2]) were performed with the Stanford-USGS SHRIMP-RG using a negative O2 primary beam and a mass resolution of 8000–8500 in order to resolve potential interferences for rare-earth elements. Concentrations were calculated from Zr2O+-normalized secondary ion yields relative to those from an in-house concentration standard (MAD, Barth and Wooden, 2010), which was calibrated to zircon standard SL13 (Mattinson et al., 2006) as well as synthetic zircon measured by electron microprobe (Claiborne et al., 2006). Zirconium concentrations in rutiles were measured by laser ablation–inductively coupled plasma mass spectrometry (LA-ICP MS) at ETH-Zurich (Marcus Walle, analyst) and are presented in Figure 5.

Mapping of Isotopic Anomalies

The spatial distribution of the most 18O-depleted rocks at the Khitostrov locality was mapped over four field seasons between 2009 and 2013 and is based on 639 mineral δ18O analyses (430 new in Table A1 in the Supplemental File [see footnote 2] and 209 analyses published in Bindeman and Serebryakov, 2011,) that form the basis for an isotope contour map (Fig. 2C). In order to determine rock δ18O values, we relied on garnet because it (1) is present in nearly all of these rocks and is an alteration-resistant mineral with a high closure temperature; (2) is chemically and isotopically inert under retrograde metamorphism, and thus best suited to record peak-metamorphic and protolithic compositions; (3) has a comparatively simple stoichiometry with little (∼0.2‰) variation in the oxygen isotope fractionation factor α among common types of garnet (Kohn and Valley, 1998); and (4) serves as a good proxy for the whole-rock oxygen isotope composition (with calculated 1000lnαWR-Grt = 0‰ and +0.5‰ for mafic and silicic Karelian lithologies, respectively; Bindeman and Serebryakov, 2011).

The new data suggest that the extent of the Khitostrov isotope anomaly is much larger than first described (Fig. 2A). Our current estimates indicate that the isotopic depletion zone now covers a >6 × 2 km area. Most 18O depletions of Khitostrov and the area inside and beyond Upper Pulongskoye Lake trace the elongated high-Fe mafic body but also extend to a wider area around it. Significant (−7 to −10‰) 18O depletion also characterizes rather ordinary looking gneiss and amphibolite without or with very limited desilication (Table A1 in the Supplemental File [see footnote 2]). Zones of maximal 18O depletion are correlated in the field with desilication trends leading to disappearance of quartz and enrichment in kyanite and corundum (e.g., Figs. 3 and 4), and these trends are used to identify and target low-δ18O areas in the field. Based on our present investigations, we suspect that expanding isotope mapping may require random sampling of ordinary looking gneisses and amphibolites, which lack visible evidence of desilication.

The second-lowest values of δ18O and δD ever documented (−20‰ and −232‰, respectively) were identified at Varatskoye, ∼20 km south of Khitostrov (Fig. 1). There as well, the isotopically anomalous rocks are associated with high-Mg amphibolites, gneisses, and a contact between gneiss and amphibolite contacts, and are accompanied by desilication leading to the enrichment in Al, Ti, and other aqueous-fluid insoluble trace elements. For Varatskoye and other localities, we previously presented isotope profiles (Bindeman and Serebryakov, 2011).

Here we also report low-δ18O values for Shueretskoye locality (Table A1 in the Supplemental File [see footnote 2], samples SH-) measured in four samples with garnets (up to 10 cm) inside of amphibolites of the southernmost Belomorian Belt, earlier described by Glebovitsky and Bushmin (1983). The Shueretskoye garnet deposit is located 150 km to the southeast along the Belomorian Belt (Fig. 1B). The discovery of low-δ18O values there expands the zone of known δ18O depletion of Karelia to 450 km. We also report δ18O analyses of eclogites of the Belomorian Belt (e.g., Shipansky et al., 2012), which represent higher-metamorphic grade, devolatilized equivalents of gneisses and gabbros of the Belomorian Belt (Table A1 in the Supplemental File [see footnote 2]).

Isotope Diversity and Zoning of Minerals

Individual minerals from single hand specimens from the Khitostrov locality are heterogeneous and zoned in δ18O by 1–2.5‰, and were interpreted to be caused by cm-scale oxygen isotope heterogeneity in the protolith, which is typical for modern hydrothermal systems (Bindeman et al., 2010). However, isotope mapping at the outcrop scale also revealed that some isotope diversity is synmetamorphic, related to fluid activity of variable δ18O (but generally ultralow δ18O with values of ∼–27 to –10‰) isotopic composition (Bindeman and Serebryakov, 2011). In particular, this diversity is observed on the contact of plagioclase-rich leucosomes (called plagioclasites) containing large corundum crystals (e.g., Figs. 3D and 4 and Table A2 in the Supplemental File [see footnote 2]). These plagioclasites are thought to have formed under the presence of synmetamorphic fluids and be metasomatic in origin (e.g., Serebryakov and Aristov, 2004; Serebryakov and Rusinov, 2004). We analyzed an individual large crystal of corundum extracted from the contact zone between melanosome and leucosome in the paragneiss using an ion microprobe, and it shows isotopic zonation with a low-δ18O core (−24‰), and heavier δ18O values closer to the rim (−21 to −22‰), next to the plagioclase leucosome (Fig. 4). Serebryakov (2004) and several other researchers consider leucosome to be metasomatic in origin rather than the result of a true partial melting. We performed Zr-in-rutile thermometry on rocks in the corundum-bearing assemblages (Fig. 5), and the temperatures range from ∼765 °C for rutiles enclosed in corundum to 720 °C for rutiles in the matrix, in agreement with upper to mid-amphibolite–facies metamorphism in these samples, and above wet solidus temperature for these rocks.

We performed melting-crystallization simulations with the MELTS (Ghiorso and Sack, 1995) program using high-Al lithologies as the protolith and 6 kbar pressure and 3–6 wt% water (Table A3 in the Supplemental File [see footnote 2]). Low degree, 15%–25%, partial melts have 12–20 wt% water and a general composition of plagioclasite, suggesting that hydrous partial melting of high-Al lithologies can generate a plagioclase-rich (up to 90%) “metasomatic” leucosome with corundum (Fig. 3D and Table A2 in the Supplemental File [see footnote 2]). Moreover, water and temperature variations in the near-solidus hydrous melting process can generate a wide range of plagioclasite and quartz-muscovite plagioclasite compositions. These are evident in crosscutting, late-stage metasomatic overprinting of preexisting lithologies and are higher in δ18O than the rocks they cut through (see Fig. 3; Serebryakov, 2004; Bindeman and Serebryakov, 2011). The increase in δ18O is then simply explained by involving normal δ18O fluids drawn into the zones of partial melting from outside the low-δ18O zones during metamorphism. This water-rich partial melt (rather than a water-dominated fluid) is also capable of concentrating REE-bearing phosphates (monazite and xenotime), which impregnate zircon rims (see below) and whole-rock assemblages with high light rare-earth element (LREE) concentrations (see below; also Terekhov, 2007).

Triple Oxygen Isotope Analysis of Ultradepleted-δ18O Rocks

The Belomorian Belt comprises rocks with the largest known terrestrial range in δ18O values. These rocks range from +9 to 10‰ in the starting Chupa gneiss to –27.3‰ in corundum-bearing rocks from Khitostrov, for which the Chupa gneiss was the protolith. The triple oxygen isotope analysis was conducted to check for the mass independent, photolytic, or extraterrestrial origin of the ultralow-δ18O rocks (e.g., oxygen having meteoritic or cometary origin with excesses or depletions in 17O; Boss, 2011). We have determined δ17O = 17O/16O − 0.5X × 18O/16O at a precision on X <0.01–0.02‰ (Fig. 6 and Table A3 in the Supplemental File [see footnote 2]), and observe strict adherence to the terrestrial mass-dependent fractionation line (Fig. 5) precluding the possibility of extraterrestrial origins. Furthermore, the triple oxygen analysis is helpful to precisely determine the nature of the mass-dependent process that caused 18O depletion and large-scale isotopic fractionation, because equilibrium or kinetic fractionation will result in different slopes (e.g., Clayton and Mayeda, 2009). The samples selected for δ17O analyses strictly adhere to the equilibrium mass-fractionation line with a slope of 0.527, in agreement with earlier work regarding hydrothermal alteration elsewhere (e.g., Young et al., 2002). Evidence of equilibrium and solid-state processes has recently been precisely determined by laser fluorination and shows 17O-18O slopes in a narrow 0.526 to 0.528 range (Rumble et al., 2007; Spicuzza et al., 2007). Kinetic isotope fractionations involving gas have shallow slopes of ∼0.516 (Young et al., 2002) or less (∼0.503; Clayton and Mayeda, 2009). Thus, the 0.527 slope that we determined for extremely diverse Karelian rocks further suggests that the extreme hydrothermal process with significant mass dissolution and loss that was responsible for formation of these rocks is characterized by equilibrium mass-dependent oxygen fractionation.

δ18O in Other Ultrahigh-Pressure Terrains

The Karelian amphibolite-grade metamorphic rocks host the world’s lowest-δ18O signature, but there are other low-δ18O localities in metamorphic terrains, including the diamond-bearing ultrahigh-pressure (UHP) rocks from Dabie-Shan Sulu (Rumble and Yui, 1998; Zheng et al., 2004; Fu et al., 2012) and Kokchetav (Kazakhstan, Masago et al., 2003). In order to explore more deeply this counterintuitive connection between isotopic depletion and UHP metamorphism, to identify causes and effects, and to test whether ultrahigh-pressure rocks are generally associated with anomalously low-δ18O values, we analyzed δ18O in minerals within diamond- and coesite-bearing samples from nine other localities (Table 2). The δ18O values for UHP rock range between ∼5‰ (Alpe Arami, Alps) and +12.5‰ (Sederonero, Greece), characteristic for hydrothermally unaltered igneous and metasedimentary sources. We thus propose that association of low-δ18O values with UHP metamorphism is purely coincidental. Presumably, these low-δ18O anomalies were discovered fortuitously because of studying UHP minerals (Dabie Shan or Kokchetav) or nearly gem-quality ruby corundum (Karelia). This survey also suggests that the ultralow-δ18O values in metamorphic and other inconspicuous rocks may be much more abundant than previously thought.

Zircon Age-δ18O Relationships for Individual Localities

The descriptions below are from localities shown in Figures 1 and 2, with data presented in Table A1 in the Supplemental File (see footnote 2) and in Figures 7–11.


We extracted zircons from corundum-bearing rocks from inside the Chupa gneiss (samples K5 and X424 in Bindeman et al., 2010 and Bindeman and Serebryakov, 2011, respectively). Backscatter electron and CL imaging reveals rounded igneous and metamorphic detrital cores sometimes overgrown by thin (<20 μm) metamorphic rims, similar to what was observed in another sample from these rocks previously imaged and dated by Serebryakov et al. (2007). New high-spatial resolution analyses of additional spots on cores and rims of zircons confirm previously published trends in δ18O versus age (Fig. 8). Specifically, our combined data set shows that (1) zircon cores are exclusively older than 2.55 Ga and are normal δ18O (∼5.3–7‰; typical for metapelitic crustal rocks, e.g., Lackey et al., 2008) with the exception of one 2.72‰ grain; (2) cores are likely detrital because of their abraded shape and range of ages; (3) a minor population of zircon cores is oscillatory-zoned with ages of ca. 2.7–2.9 Ga mantled by uniformly gray CL ca. 2.6 Ga domains; both are normal in δ18O (crystal KV10_1,2, Fig. 7); (4) low-δ18O zircon rims are of Svecofennian metamorphic age (ca. 1.85 Ga) with δ18O ranging from –27.3‰ to values identical with normal δ18O cores; (5) most Khitostrov zircons (1.8 Ga rims and 2.6 Ga internal zones) have Th/U ratios <0.1 indicative of metamorphic growth (Hoskin and Schaltegger, 2003; Rubatto et al., 2009), with the exception of the ca. 2.7–2.9 Ga oscillatory-zoned zircon cores (Fig. 8B). Mostly low Th/U zircon crystals from corundum-bearing paragneiss contrast with higher Th/U zircon from neighboring mafic intrusions. These yield mostly ca. 1.85 Ga Svecofennian metamorphic ages (Table 1). We attribute the strong variability of U, Th, and REE in metamorphic zircon crystals to co-crystallization of coeval monazite (1.87 Ga; Bindeman et al., 2010) competing with zircon for Th and light REE (see further discussion below).


Zircons in Varatskoye (Tables 1 and 2 and Fig. 9) were extracted from corundum-bearing rocks surrounded by Chupa paragneiss and from amphibolites that represent a metamorphosed ca. 2.4 Ga mafic intrusion. Our amphibolite sample is compositionally similar to a low-grade metamorphic mafic intrusion at 2 km distance from the Varatskoye location, which yielded U-Pb zircon ages of 2.40 Ga (Bibikova et al., 2004). Large zircons in Varatskoye rocks permitted co-registered U-Pb age and δ18O analysis of clearly separated core and rim domains. Zircons in corundum-bearing rocks with gneiss as protolith display straightforward relationships between δ18O and age: cores are exclusively 2.8–2.55 Ga with δ18O = 5–7‰, and rim ages are 1.9–1.8 Ga with δ18O = −19‰. The youngest cores interpreted to be metamorphic are ca. 2.55 Ga, and, similar to Khitostrov, are characterized by dark, featureless CL (Fig. 7). Zircons in the Varatskoye amphibolites were exclusively of Svecofennian age and display subtle rimward increases in δ18O (Figs. 9B and 9C). Within the age range between 1.89 and 1.85 Ga, we observe a nearly linear increase in δ18O of ∼1–1.5‰. This temporal variation is too large to be caused by increasing zircon–whole-rock fractionation (1000lnα18O(zircon-WR) during retrograde metamorphism and instead reflects a subtle increase in δ18O for the metamorphic or intergranular fluids (e.g., Lackey et al., 2008; Peck et al., 2010; D’Errico et al., 2012). Unlike predominantly high-U concentrations and low-Th/U ratios at Khitostrov, at Varatskoye, U concentrations are mostly low (several to tens of ppm, resulting in comparatively large errors in the U-Pb ages), whereas Th/U ratios are highly variable.

Height 128 and Lyagkomina

Zircons in these two localities were extracted from corundum-bearing paragneiss samples. Height 128 zircons have dark CL, rounded cores with some oscillatory zoning and uniformly bright-CL rims that yielded 1.85 Ga ages (Figs. 7 and 10). The thickness of these Svecofennian rims varies from grain to grain, and they are extremely U and Th poor (<10 ppm) causing comparatively large U-Pb age uncertainties (Table 1). There are also few distinct intermediate zircon domains with dark CL and low Th/U that are metamorphic with 2.59–2.6 Ga ages. Zircon rims with dark CL and low Th/U have δ18O = −9‰ (Figs. 7 and 10), in equilibrium with garnet and other minerals of the host rock. Lyagkomina (sample L-1) zircons have oscillatory-zoned igneous cores with 2.7–2.75 Ga ages and δ18O = 6.5–7.7‰ (n = 22) and the δ18O at the rims of −4‰, also in equilibrium with the host assemblage. The Th/U ratios of interiors are typical for igneous crystallization, whereas rims are bright in CL and are low in Th/U and δ18O values. In contrast to Height 128, crystal domains of 2.55–2.65 Ga age are absent. Intermediate age and δ18O values (Fig. 9) are the result of beam overlap onto core and rim domains (marked C/R in Table 1), and are thus not considered further.

Mafic Intrusions

Two main types of mafic intrusions regionally occur in Fennoscandia (Figs. 1 and 2, and Fig. A2 in the Supplemental File [see footnote 2]) and are abundantly represented in the area: high-Mg gabbros (druzites) attributed to the 2.4–2.5 Ga global rifting episode (Amelin et al., 1995; Puchtel et al., 1997; Sharkov et al., 1997; Stepanova and Stepanov, 2010) and less abundant, high-Fe tholeiitic dikes tentatively dated at 2.1 Ga (Stepanov, 1981; Stepanova et al., 2003). The respective ages of 2.4–2.5 Ga and 2.1 Ga are thus typically assigned in the field on the basis of high-Mg and high-Fe composition. In the Belomorian Belt, both intrusion types are variably overprinted by 1.85 Ga Svecofennian metamorphism, and localities with extreme 18O depletions also have intense chemical modification trends of desilication and aluminum enrichment. However, the composition-age designation, although more difficult to see, is still recognizable (e.g., sample X451, Table 1; Kulezhma KY21, Fig. 11). Both types of high-Mg and high-Fe intrusives sometimes occur in close proximity in the studied localities (e.g., at Kulezhma and Khitostrov).

In an attempt to determine the U-Pb age of mafic intrusions that host, or are in close proximity to, low-δ18O localities, we extracted zircons from hand samples collected near the center of these bodies (Fig. 11) where normal, mantle-like δ18O values for all major minerals prevail. The majority of zircon interior ages from Dyadina Gora, Khitostrov, and Varatskoye (Figs. 9 and 11) yielded 1.85 Ga Svecofennian ages with the exception of a single concordant zircon interior from Dyadina Gora, which returned a 2.39 ± 0.022 Ga (2σ) 207Pb/206Pb age. The least metamorphically modified high-Mg intrusive body at the Keret’ River has been dated at 2.40 Ga in close proximity to the Varatskoye locality (Bibikova et al., 2004). A single slightly discordant zircon age from inside of high-Fe intrusion at Khitostrov gave a 207Pb/206Pb age of 2.228 ± 0.014 (2σ) Ga (sample X245, Fig. 10). Unfortunately, most zircons from this sample and all zircons extracted from high-Fe amphibolite (sample X451, Table 1) in close proximity to the contact with lowest-δ18O corundum-bearing rocks yielded ca. 1.85 Ga Svecofennian ages indicative of overprinting and metamorphic zircon recrystallization. Zircon interiors in a high-Fe intrusion at Kulezhma were better preserved and yielded 207Pb/206Pb zircon core ages ranging between 2.15 and 2.04 Ga with highly variable Th/U ratios, and minor 1.85 Ga zircons, which suggest that some of these zircon ages represent metamorphic recrystallization. However, a recent compilation of intrusion ages determined by various zircon dating methods (mostly by ion microprobe) suggests that high-Fe dikes could also be older than 2.1 Ga (see Hanski and Melezhik, 2013, figs. 3.8 and 3.9 and references therein). Our single-zircon core age for Khitostrov at 2.23 Ga could potentially also be a minimum age (Fig. 11B).

High-Spatial Resolution Core-Rim δ18O Relationships in Zircons

The drastic isotopic gradients between normal δ18O (detrital) zircon cores and ultralow-δ18O rims in nearly all gneiss samples (Fig. 7) detected by conventional spot analyses are limited in spatial resolution to the lateral beam dimensions used (∼25 µm). To constrain the spatial scale over which the δ18O transition occurs, and to assess if this correlates with the sharp (<1 µm) CL contrast between cores and rims, we performed a high-spatial resolution isotope profiling using a ∼0.5 µm diameter Ga+ primary ion beam (Fig. 12). A zircon crystal from Varatskoye corundum-bearing rocks was selected because it had one of the thickest rims detected by CL imaging (Fig. 12). Over the length of the profile (∼25 µm), the entire ∼23‰ shift in δ18O values occurs over a 5–6 µm interval. The shape of the δ18O data array is symmetric, resembling a “Fickian” diffusion profile that is centered at the dark-light CL boundary between core and rim.

Rare-Earth Element Analyses of Zircons

Zircon rare-earth element (REE) abundances were determined synchronously with U-Pb ages for a subset of zircons (Fig. 13). In contrast to “normal” magmatic or metamorphic zircon REE patterns that are characteristically depleted in LREE, enriched in high HREE, and commonly have pronounced positive Ce and negative Eu anomalies (e.g., Rubatto and Hermann, 2007), many zircons from the low-δ18O gneisses show a relative overabundance of LREEs and medium rare-earth elements (MREEs) (Fig. 13), sometimes with a positive Eu anomaly and a lack of a Ce anomaly (e.g., samples K5, KV10, and L1). Light rare-earth element–enriched zircons were also recently reported by Krylov et al. (2012) for Khitostrov locality using different instrumentation. Zircons from gabbros (samples X245 and DG150) also display concave-up REE patterns. To our knowledge, these are the most anomalous distributions observed for zircons anywhere in the world, even exceeding similar trends in detrital Hadean zircons from Western Australia (Cavosie et al., 2006). The rare exceptions to this anomalous behavior are zircons in sample B51 and several zircon cores (e.g., in sample K5; Table 1), which match the zircon-typical REE patterns (Hoskin and Schaltegger, 2003). Excess LREEs in zircon were interpreted to indicate “hydrothermal” or “hydrothermal-pegmatitic” origins (Hoskin and Schaltegger, 2003) because such compositions deviate from magmatic zircon-silicate melt partitioning (e.g., Rubatto and Hermann, 2007; Reid et al., 2011). Pettke et al. (2005), however, observed similar abundances of LREE in magmatic and hydrothermal zircons but found that hydrothermal zircons are characterized by strongly negative Eu anomalies, a feature entirely absent in zircons studied here. Co-crystallization of zircon with garnet cannot explain the observed relationships because garnet depletes the HREE.

At face value, chondrite-normalized LREE abundances of 10–500 suggest a “super”-pegmatitic parental solution or melt containing weight percent REE abundances. We consider such an interpretation improbable because: (1) no such melt has been described in nature and is uncharacteristic at 600–720 °C (e.g., Hermann et al., 2013), and (2) the anomalous distribution of LREE and MREE characterizes not only synmetamorphic zircon rims of Svecofennian 1.85 Ga age but also variably affects 2.9–2.5 Ga detrital zircon cores (see an example of diversity in REE profiles in sample K5).

The best explanation of this unusual phenomenon is the presence of an LREE- and MREE-enriched phase such as monazite (CePO4) or xenotime (YPO4), which contains 500,000 ppm of LREE and Y, and thus even tiny amounts can contaminate the zircon analysis. Monazite is present in analyzed rocks, and it has been dated specifically in sample K5 to yield Svecofennian 1.89 Ga age (Bindeman et al., 2010). We examined zircon crystals and analysis craters using electron beam imaging (secondary electrons [SE] and backscatter electrons [BSE]) to test whether the unusually large LREE enrichments are due to monazite or xenotime inclusions in zircon, but we found no recognizable monazite. If such inclusions existed, they would be smaller than detectable by SE or BSE imaging (<100 nm). To further test the hypothesis that accumulation of nanoinclusions occurred inside metamict, radioactively damaged zones, we analyzed selected spots in sample K5 using a CAMECA SX100 electron probe with wavelength dispersive spectrometers (calibrated relative to synthetic zircons, Fig. 13C; see Fig. A5 in the Supplemental File [see footnote 2]), exciting X-rays with high beam current (300 nA) and voltage (30 kV) to achieve high precision for trace elements (calculated detection limits are: P = 20 ppm, Ce = 90 ppm, and Y = 130 ppm). Because monazite and xenotime are both phosphate minerals, we concentrated on detecting phosphorous and correlations between REE and P. Analyzed zircons had up to 900 ppm P, and a positive correlation between Y and P exists. However, correlations between P and especially Ce are subtle. Although the amount of P in K5 zircons is abnormally high compared to common zircon (e.g., Hoskin and Schaltegger, 2003), it is still too low to account for the total abundance of REE (Fig. 13). This requires non-phosphate REE-rich phases in these zircons.

These measurements support our inference that zircons are “infected” by nanometer-scale inclusions, at least some of which would be monazite or xenotime like. Monazite recrystallizes and is dissolved in fluids at lower temperatures than zircons during regional metamorphism (Rubatto et al., 2001). The hydrous silicate melts that form plagioclasites, the leucosome in studied localities, serve as an appropriate transport media for these REE-rich phases. Our interpretations are in line with those of Cavosie et al. (2004, 2006), who attributed more modest LREE enrichment in >3.9 Ga Hadean detrital zircons to microinclusions formed at very low water/rock ratio inside of metamict zircon zones. These authors also found positive albeit scattered correlations between P and REE abundances.

High LREE abundances in zircons correspond to elevated and highly fractionated LREEs/heavy rare-earth elements (HREE) in corundum-bearing rocks across Karelia, which were reported by Terekhov (2007). Presence of just 0.01 wt% monazite (as detected in these rocks) is capable of explaining the elevated LREE budget.

The δ18O Analyses of Nonmetamorphic Protoliths

In a reconnaissance study of coeval synrift sedimentary and igneous rocks that may preserve direct evidence for interaction with ultra-depleted glacial melt waters, we have sampled low-grade metamorphic rocks to the southeast (Vetreny Belt) and southwest (north of Petrozavodsk) of the Belomorian Belt (Fig. 1). We analyzed materials that reasonably could have interacted with low-δ18O glacial meltwaters: glacial tillite, secondary minerals (amygdaloids and quartz vugs) inside basaltic lavas intercalated with glacial deposits, alteration minerals formed between margins of basaltic pillows, and quenched pillow basalt rinds (Table 3 and Fig. 14). For comparative purposes, we have also analyzed secondary quartz from Antarctica’s Minna Bluff area, which represents subglacial and proximal alteration at 12–9 Ma (Antibus et al., 2012). Alteration history and isotope thermometry of these rocks have been studied in detail and confirm a 10–100 °C temperature window for deposition of quartz (Antibus, 2012), which constrains quartz-water isotope fractionation.

If alteration and secondary products in Karelian rocks were formed at low-temperatures (50° ± 30 °C) as suggested by comparison with Minna Bluff Antarctic quartz data, then the altering waters must have been ultra-low-δ18O, comparable to Antarctic ice. The Karelian quartz has δ18O values that are 10–15‰ lower than subaerial Oligocene quartz vugs in Oregon, formed in temperate near-coastal conditions (Fig. 14).

Summary of Zircon Age-δ18O Relationships

There are clear spatial and temporal patterns for δ18O in zircon from corundum-bearing rocks in the Chupa gneiss that are consistent for all sampled localities: normal δ18O cores with ages ranging from >2.9 to 2.55 Ga and ultralow-δ18O rims with younger ca. 1.85 Ga Svecofennian ages. Zircon-rim δ18O values are in isotopic equilibrium with the host-rock mineral assemblage, ranging from δ18O = −4.1‰ in Lyagkomina to –27.3‰ in Khitostrov, which to our knowledge is the lowest value ever reported (samples X424 and X425). Rims also have universally low Th/U ratios characteristic for metamorphic zircon, whereas cores have variable but generally higher “igneous” Th/U. Some zircon cores have an internal mantle with dark CL and low Th/U values that date to a ca. 2.6 Ga metamorphic event previously reported for the Belomorian Belt (Bibikova et al., 1994, 2001, 2004). Importantly, these ca. 2.55–2.6 Ga internal mantle domains are normal in δ18O, and no unambiguously igneous zircon of this age is present in our data. The morphology of 2.6–2.55 Ga zircon cores is rounded and non-euhedral (Fig. 7) and represents metamorphic recrystallization of older cores. The sharp termination of 2.55–2.6 Ga metamorphic zircon growth zones suggests cooling potentially due to uplift and unroofing of the Chupa gneiss. The amount of denudation after 2.55 Ga but before regional rifting at 2.4 Ga remains unknown, but rocks must have resided at this time at depths shallow enough to enable open fractures and attendant alteration by meteoric water at large water/rock ratios to imprint the observed ultralow-δ18O glacial meltwater signature on them shortly thereafter.

The next episodes of zircon growth that reflect a regional geologic event are recorded exclusively in gabbroic intrusions but not in the gneisses they intrude. The zircon record of the high-Mg gabbros yields only a single magmatic core of 2.4 Ga (i.e., in DG150-11, Fig. 11), which is in agreement with similarly aged zircon crystals described for mafic intrusions at locations that are not surrounded by low-δ18O depletions (Puchtel et al., 1997; Bibikova et al., 2004). Rare zircons in high-Fe intrusions (e.g., 2.23Ga in Khitostrov) characterize less abundant high-Fe magmatism with nominal 2.1Ga age but may be locally older (ca 2.2–2.3 Ga, Hanski and Melezhik, 2013). Pervasive 1.85 Ga Svecofennian metamorphism appears to have obliterated much of the pre–1.85 Ga zircon record in the mafic intrusives immediately adjacent to the low-δ18O localities. This is not surprising because the odds for survival and detection of any older igneous or metamorphic zircons in such mafic intrusions, if they were ever present, are low because of low Zr abundance and small grain sizes typical for such rocks. It should be mentioned, however, that the cores are better preserved (and dated) in the least metamorphosed mafic intrusions 1–2 km away from localities of interest (2.4 Ga at Varatskoye; Bibikova et al., 1994, 2004).

The last episode of zircon growth is during the 1.85 Ga Svecofennian metamorphism, affecting zircons in both gneisses and mafic intrusions. This stage is represented by zircon overgrowths on existing cores and crystallization of new zircons (Fig. 7).

Zircon geochronology reveals a rather short (∼50–100 m.y.) time gap between the 2.55 Ga youngest metamorphic crystallization of zircons in the gneiss and their inferred shallow hydrothermal alteration by glacial meltwaters near the surface by 2.45–2.4 Ga in the vicinity of superplume-related (Evans and Pisarevsky, 2008), high-Mg gabbroic intrusions. Two scenarios are possible. (1) Cessation of zircon growth after 2.55 Ga but before 2.45 Ga can be reasonably achieved at comparatively slow metamorphic unroofing (and cooling) rates of ∼1 mm/yr: uplift from amphibolite-grade depths (15–18 km) to the near surface would take ∼15–36 m.y. Unroofing in response to postcollisional erosion-driven isostatic uplift (removal of the mountain top overburden) or by lateral gravitational flow of the thickened crust, as is observed in the Alps (Ruppel et al., 1988; Champagnac et al., 2009), is conceivable. At such rates, uplift and cooling (Fig. 14) will be completed ∼50 m.y. before the superplume event and the inferred initiation of rifting at 2.5 to 2.45 Ga. (2) Uplift caused by extension during the 2.5–2.4 Ga rifting episode, a known regional (Rybakov et al., 2000) and global event. In this scenario, rifting resulted in rapid exhumation from mid-crustal depths. For this tectonically driven uplift history, exhumation would have been completed in less than 10 m.y.

The Khitostrov locality may require hydrothermal alteration to take place at 2.23 Ga, coincident with high-Fe gabbro intrusions.

Origin of Sharp Isotope Boundaries inside Zircon

Karelian zircons show sharp isotopic shifts across the core-rim boundaries, with the rims being in isotopic equilibrium with the host assemblage (Fig. 12). For a ∼3 µm half-thickness of the δ18O diffusion profile in the Varatskoye zircon, and assuming that the zircons spent at least 10 m.y. at high peak metamorphic conditions of ∼700 °C indicated by the metamorphic paragenesis and Zr-in-rutile geothermometry (Fig. 5), the estimated oxygen-diffusion coefficients are ∼10−23 to 10−24 m2/s, in agreement with Watson and Cherniak’s (1997) “dry” diffusion coefficients for zircon. For “wet” diffusion coefficients, complete annealing of the zircon crystals would occur over the same timescale, even at a lower metamorphic temperature of ∼600 °C, which is contrary to our observations. “Wet” diffusion would produce the observed diffusion profile in ∼0.1 m.y., a duration that is too brief to be reconciled with reasonable cooling rates for regional retrograde metamorphism. There are several considerations here for explaining the comparatively sharp δ18O diffusion profile we have observed: (1) experimental uncertainties in wet versus dry diffusion coefficients, in particular the “wet” diffusion coefficients of Watson and Cherniak (1997) being too large (e.g., Page et al., 2007; Bowman et al., 2011); (2) the metamorphism occurred under drier conditions than usually assumed for such rocks (Kohn, 1999); (3) metamorphic rims crystallized at much lower “metasomatic” temperatures (<500 °C), where diffusion was vanishingly small3; and/or (4) the metamorphic event was unusually short. Under these conditions, zircon-rim crystallization should have occurred at a maximum temperature 510–485 °C for only a few million years. However, such low temperatures contradict metamorphic grade, Zr-in-rutile temperatures (Fig. 5), and Ti in zircon temperatures recently reported by Krylov et al. (2012) for the Khitostrov locality, all suggesting 650–770 °C for zircon-rim crystallization. The age difference between amphibolitic zircon and monazite crystallization age of 1.89 Ga (Bindeman et al., 2010) and inferred ca. 1.75 Ga ages of exhumation of these rocks (Terekhov, 2007) suggests a normal duration of metamorphism and exhumation, lasting tens of millions of years. We thus prefer explanation (1) because it is in agreement with geological evidence (e.g., Page et al., 2007; Bowman et al., 2011).

Search for a Nonmetamorphic Low-δ18O Protolith

The Karelian ultralow-δ18O anomaly is exclusively hosted by amphibolite-grade metamorphic rocks of the 1.85 Ga Belomorian Belt extending >450 km, but oxygen isotope depletion of the protolith is likely associated with a Paleoproterozoic “Slushball Earth” episode of subglacial alteration around rift-related mafic intrusions (Fig. 1). No ultralow-δ18O rocks have yet been reported in coeval low-grade volcanic and sedimentary rocks (Table 3 and Fig. 14). Although ultralow-δ18O supracrustal rocks would be the ultimate proof for extremely 18O-depleted glacial waters (which, if present at low latitudes would indicate a “Slushball Earth” condition; Hoffman, 2009; Bindeman et al., 2010), we have reason to surmise that such evidence would be difficult to find, or be uncommon. First, this is because water-rock interaction is kinetically restricted in supracrustal rocks, and isotope fractionation factors are large at low temperatures, favoring higher-δ18O solids. In contrast, isotopic exchange is rapid and extensive in hydrothermal systems, and 1000lnαrock-water whole-rock–water fractionation is close to 0‰, thus more reliably recording water δ18O values at large water/rock ratios. Second, there are uncertainties regarding the geological position of sediments and lavas as representing subglacial or lacustrine environments (low-δ18O water), or marine environments (∼0‰ δ18O SMOW seawater values). Nonetheless, despite the unknown temperature of formation, the δ18O values in supracrustal samples from Karelia (Fig. 14) are at the lower end of geologically more recent analogues, overlapping with the Antarctic quartz data.

Protolithic versus Synmetamorphic Low-δ18O Signatures

Our isotope mapping for the Khitostrov zone of ultralow-δ18O rocks shows concentric zonation in the proximity of elongated mafic intrusion, sheared during 1.85 Ga Svecofennian metamorphism (Fig. 2). If the elongated concentric isotopic pattern were to be explained by infiltration of ultralow-δ18O synmetamorphic fluids at 1.85 Ga (as suggested by Terekhov, 2007; Ustinov et al., 2008), the following very specific conditions must be met.

  1. In the absence of any recognized reservoirs of mantle- or crustal-derived, low-δ18O and low-δD, high-temperature fluids, the only conceivable source of such ultralow-δ18O fluids is devolatilization of (previously surface-exposed) buried metamorphic rocks.

  2. Because typical gneisses and amphibolites contain only 1–2 wt% H2O (Fig. 15A) and because at amphibolite-grade temperatures of 600–700 °C, 1000lnαwater-rock water–whole-rock fractionation is close to zero (Fig. 15B), the devolatilizing protolith would have to be: (a) essentially the same ultralow-δ18O composition to yield ultralow-δ18O fluids that are inferred for Karelia; and (b) it would take 50–100 times (Fig. 15D) the mass of the devolatilizing ultralow-δ18O rock to produce cubic kilometers of –10 to –27‰ rocks documented in outcrops (Fig. 2).

  3. Isotopic effects are maximized only if fluids were escaping through the same metamorphic shear zone, thus integrating fluid/rock ratios (Figs. 14C and 14D). On an equimolar-oxygen basis, and assuming that fluids are following the same path with 100% exchange efficiency, it will take equal amounts of fluids to rocks to bring δ18O value of rocks closer to that of the fluid, again requiring 100 times the volume of the known low-δ18O exposures (6 × 2 km, Fig. 2C). If the isotope exchange reaction is only 50% efficient, and/or if increments of fluids are following different (non-integrating) paths, the total isotopic effects are significantly reduced, and thus the required amount of devolatilizing rocks at depth is significantly increased.

  4. The flushing fluids will induce progressively lower (but still significant) isotopic effects on rocks down the path of exchange (up the fault). Thus the exposed –27‰ rocks would require an even lower δ18O devolatilizing source (Fig. 15C).

  5. If metamorphic fluids were CO2 rich (e.g., Valley, 1986; Lackey and Valley, 2004), they should have been derived from decarbonation of ultralow-δ18O carbonate, unknown in the area or worldwide. While mass balance would require less rock on an equimolar-oxygen basis (∼10–20 times the amount of decarbonating marble versus 50–100 times the amount of dehydrating amphibole-bearing rock), the fractionation factors required are still problematic: δ18O of CO2 is heavier than coexisting calcite by 3.3‰ (900 °C) to 6‰ (700 °C; Zheng, 1993; Rosenbaum, 1994), thus diminishing the overall fluxing effect of CO2 and requiring an even lighter devolatilizing carbonate source of –30 to −35‰ δ18O SMOW.

In summary, a synmetamorphic fluid-fluxing hypothesis would require >100 times 18O-depleted rocks underneath the Khitostrov and other Karelian localities, making an even stronger case for widespread pre–1.85 Ga surface altered rocks. Additionally, investigated eclogites of the Belomorian Belt (Fig. 1A and Table A1 in the Supplemental File [see footnote 2]), devolatilized equivalents of studied rocks, are structurally below the examined localities, but these eclogites lack low-δ18O minerals.

We consider shallow heated glacial meltwater interaction with rocks as the most parsimonious explanation. This would result in a circular symmetric depletion pattern formed around mafic intrusions prior to 1.85 Ga metamorphism, which, in turn, deformed these depletion zones into elongated and thus likely fault-controlled, ultralow-δ18O localities. Later local devolatilization of the low-δ18O protolith during metamorphism generated zones of equally low-δ18O synmetamorphic fluids, plagioclasite leucosomes (Fig. 3), which caused localized formation of zoned metamorphic minerals, zircon rims, smoothed and/or obscured preexisting δ18O zonation in outcrops, enriching metamict zones in zircons, and whole rocks with LREEs.

Alternative Isotopic Ways to Produce Ultralow-δ18O Fluids?

Are there other conceivable processes that can produce the isotopic ranges observed? We are aware of only three processes other than hydrothermal alteration that are capable of producing large (>20‰) depletion of δ18O: (1) kinetic isotope fractionation, e.g., during devolatilization or evaporation with significant mass loss (Clayton and Mayeda, 2009; Mendybayev et al., 2010); (2) thermal diffusion in water-bearing rocks (Bindeman et al., 2013); and (3) extreme Rayleigh distillation.

Rapid thermal decomposition of hydrous phases (e.g., brucite and serpentine; Clayton and Mayeda, 2009) at low pressure or vacuum involves kinetic removal of increments of low-δ18O water. The process can be viewed as isotope disproportionation into high-δ18O residue and low-δ18O fluid. However, contrary to the Karelian data, such a process generates shallow “kinetic” slopes of δ17O versus δ18O fractionation of 0.503–0.516, which are different from the equilibrium mass-dependent slopes that we measured (0.527, Fig. 6). The kinetic devolatilization process is also unlikely to explain the genesis of the Karelia rocks because high-δ18O anhydrous residues are absent. Moreover, evidence is lacking that high-pressure metamorphic devolatilization would be kinetically equivalent to low-pressure devolatilization.

For thermal isotope redistribution (or “thermal migration”; e.g., Lundstrom, 2009), the sense of O- and H-isotope redistribution in a temperature gradient is low δ18O at a hotter end and high δ18O at the colder end, in spatial agreement to the observed contact relationships between gabbro and gneiss. Covariant mass-dependent fractionations as much as 28‰ for δ18O and 144‰ for δD are possible in a ∼500 °C temperature gradient (Bindeman et al., 2013) using natural, normal-δ18O rocks with 2–4 wt% H2O. As is the case with kinetic devolatilization above, both low- and high-δ18O values are generated in the process, while Karelia rocks universally show a decreasing trend from normal-δ18O gneiss or gabbro to ultralow-δ18O values in corundum-bearing rocks. However, further tests can be applied to explore both of these possibilities: oxygen and hydrogen isotopes should be covariant with other isotopic species (e.g., Si, Mg, and Fe) in a predictable mass-dependent way (Mendybayev et al., 2010; Bindeman et al., 2013).

Finally, an extreme Rayleigh fractionation process of incremental fluid removal, which we can call the “underground distillery model,” can accentuate isotopic differences. However, at high temperatures, isotope fractionation factors are small (Fig. 15B), the amount of ultra-distilled fluid is measured at a few percent of the original, and thus the possibility of any significant “underground distillation” is also unlikely. We suggest that interaction between normal δ18O rocks and ultralow-δ18O glacial meltwater presents the simplest and most realistic explanation to the observed results.

Insights into Duration of Paleoproterozoic Slushball Earth Glaciations

The Paleoproterozoic glaciation on different continents is inferred to have lasted from ca. 2.45 Ga to as late as 2.26 Ga based on existing and newly appearing geochronologic data on diamictites worldwide (see Hoffman, 2013, for review). Three to four individual glacial episodes, each of unknown but likely global or near-global extent, characterized the Earth during this time interval, each lasting multiple millions of years. Coeval to glaciation was the appearance of oxygen in the Earth’s atmosphere and disappearance of the mass-independent sulfur photolysis reactions between 2.4 and 2.26 Ga (Bekker et al., 2004). When exactly the Great Oxidation Event (GOE) occurred remains a matter of debate.

Our work contributes to this discussion because ultralow-δ18O values in Karelian rocks provide direct evidence for low-δ18O glacial meltwaters and thus terrestrial glaciation at low latitudes. The majority of Karelia localities record a 2.4 Ga episode, during coeval rifting and high-Mg plume magmatism (Fig. 16; also Bindeman and Serebryakov, 2011). The new evidence from the most depleted locality at Khitostrov is adjacent to high-Fe gabbro that yielded a single-zircon core of 2.23 Ga, suggesting that hydrothermal alteration could have been caused by the youngest glaciation dated at 2.26 Ga by Rasmussen et al. (2013) in South Africa and Canada. Thus there is a need to date post–2.4 Ga high-Fe intrusions in the Belomorian Belt and their unmetamorphosed variety in Karelia to confirm the age designation. If this evidence is confirmed, then the Karelian low-δ18O rocks record the oldest and the youngest of glaciations during the Paleoproterozoic, so that every shallow intrusion likely underwent subglacial or near-glacial meteoric-hydrothermal alteration.

  1. Eleven newly discovered low- and ultralow-δ18O Paleoproterozoic Karelian localities extend the previously known geographical range of such rocks in the Belomorian Belt to 450 km.

  2. At the Khitostrov locality, which hosts the world’s lowest δ18O rocks (−27.3‰), the mapped zone of depletion is now extended to ∼6 × 2 km, tracing the high-Fe gabbro body in an exposed regional fault.

  3. Isotopic mass balance supports the idea of near-surface alteration by glacial meltwaters at large water-rock ratio as the most likely mechanistic interpretation of the observed depletion patterns, prior to metamorphism.

  4. Zircon crystals in corundum-bearing rocks associated with the Chupa paragneiss display systematic zoning patterns with normal-δ18O, 2.9–2.7 Ga igneous cores that are commonly mantled by normal-δ18O 2.55–2.6 Ga metamorphic domains and ultralow-δ18O Svecofennian 1.85 Ga metamorphic rims. Zircons in metamorphosed mafic intrusions are predominantly of the younger metamorphic Svecofennian age.

  5. Zircon ages of gabbroic intrusions in most low-δ18O localities indicate intrusive ages coincident with 2.45 Ga rifting and the oldest Paleoproterozoic glaciation. Hydrothermal alteration in the rift zones involved heated glacial meltwaters at large water-rock ratios, implying shallow residence of the studied localities at that time.

  6. Whether the depletion at Khitostrov is younger than 2.4 Ga, and is associated with ca. 2.23 Ga high-Fe intrusions, needs to be further tested. If proven, Khitostrov will document a hydrothermal alteration event during the youngest of Paleoproterozoic glaciations.

  7. Svecofennian 1.85 Ga burial and metamorphism have resulted in very limited devolatilization of low-δ18O rocks diluted by the fluids from ambient normal-δ18O rocks, leading to increasing δ18O in metamorphic minerals and zircon rims. This small fluid-rock ratio or hydrous low-degree melting may explain unusually enriched REE and P concentrations in zircons. We invoke dispersed LREE-rich nanometer-sized REE-phosphate inclusions in metamict zones within these zircons.

  8. Geochronologic and geologic evidence indicates depletion in δ18O happened between 2.55 and 1.85 Ga, which is broadly coincidental with glacial deposits of the Sariolian and Sumian age of Karelia. These coeval volcanic and sedimentary rocks and their alteration products have δ18O values broadly comparable to those in modern Antarctica but lack ultralow-δ18O values due to lower (and uncertain) temperatures of alteration as compared to high-temperature hydrothermal alteration around magmatic intrusions.

  9. The redefined range of δ18O in Karelian rocks (from +10 to −27.3‰) allows a more precise determination of exponents of the δ17O versus δ18O fractionation exponent, as 0.527, in strict adherence to the equilibrium distribution of isotopes in the process, suggesting that intense hydrothermal alteration obeys equilibrium mass dependency.

We thank the National Science Foundation for funding of grants EAR-1049351 and EAR-CAREER 0805972; IF-EAR funding of the CAMECA 1270 ion microprobe facility; the University of California, Los Angeles; and the Russian Foundation for Basic Research for grant 12-05-00706a. We also thank Pavel Medvedev, Sergei Svetov, and Denis Korpechkov for help during fieldwork of 2011 and 2012, and Markus Walle for help with LA-ICP MS. V. Kulikov kindly supplied samples from Vetreny Belt, and Kurt Panter is thanked for collaboration on Antarctic secondary quartz. J. Eiler and Caltech Ion microprobe facility are thanked for hosting I. Bindeman during sabbatical, and Cliff Dax is thanked for electronics support. Associate Editor Aaron Cavosie and two anonymous reviewers are thanked for their careful reviews.

1Here and below, “normal-δ18O” rocks are defined as being in the +5.5 to 6.5‰ range characteristic of δ18O values for mantle-derived basic rocks and silicic products of their differentiation; differentiation causes subtle increase of δ18O with increasing SiO2 (“mantle array” of Bindeman, 2008). Normal-δ18O minerals and fluids are in high-temperature (T) isotopic equilibrium with (and within) these rocks; their δ18O values are different from the bulk rock by small (typically less than 1–2‰ at high T) mineral-specific, temperature-dependent fractionation factors (e.g., Taylor, 1974). Physical weathering and isochemical metamorphism do not change bulk δ18O value of rocks despite changes in modal mineral identity. This paper also refers to high-δ18O gneisses, with values >7.5‰, formed by metamorphism of high-δ18O supracrustal sedimentary protoliths. Low-δ18O rocks (0 to 5.5‰) represent the result of hydrothermal interaction of low-δ18O meteoric water with normal- and high-δ18O rocks and minerals, leading to heterogeneous lowering of δ18O values in rocks. The “ultralow”-δ18O rocks featured in this paper (δ18O <<0‰ standard mean ocean water [SMOW]) are formed by high-T interaction of “ultralow”-δ18O glacial meltwater with rocks.
2Supplemental File. Laser fluorination analyses and supplementary tables and figures. If you are viewing the PDF of this paper or reading it offline, please visit or the full-text article on to view the supplemental file.
3To quantify the maximum temperature at which the zircon rims could have crystallized without significantly reequilibrating with the zircon cores, we modeled metamorphic cooling using fast grain boundary diffusion model (Eiler et al., 1992; Peck et al., 2003) over reasonable metamorphic temperatures and cooling duration of 10 m.y. using “wet” diffusion coefficients (Watson and Cherniak, 1997); see figure A1 in Bindeman and Serebryakov (2011).