Perturbation of the deep-Earth carbon cycle in response to the Cambrian Explosion.

Earth's carbon cycle is strongly influenced by subduction of sedimentary material into the mantle. The composition of the sedimentary subduction flux has changed considerably over Earth's history, but the impact of these changes on the mantle carbon cycle is unclear. Here, we show that the carbon isotopes of kimberlite magmas record a fundamental change in their deep-mantle source compositions during the Phanerozoic Eon. The 13C/12C of kimberlites before ~250 Ma preserves typical mantle values, whereas younger kimberlites exhibit lower and more variable ratios-a switch coincident with a recognized surge in kimberlite magmatism. We attribute these changes to increased deep subduction of organic carbon with low 13C/12C following the Cambrian Explosion when organic carbon deposition in marine sediments increased significantly. These observations demonstrate that biogeochemical processes at Earth's surface have a profound influence on the deep mantle, revealing an integral link between the deep and shallow carbon cycles.

INTRODUCTION

Earth’s carbon cycle provides a first-order control on the concentration of oxygen and carbon dioxide in the atmosphere and oceans and, as such, is essential in producing and maintaining a habitable planet. The carbon cycle operates at many different levels, and on a planetary scale, the process of subduction mediates the transfer of surface material, including carbon, into Earth’s mantle. Major changes in the physical, chemical, and biological conditions at Earth’s surface may thus be expected to exert a profound influence on the planet’s interior (). Similarly, temporal variations in the carbon content or isotopic composition of the mantle could be used to track changes in Earth’s surface carbon cycle.

Knowledge of the deep carbon cycle and its evolution through time, however, is presently incomplete. Diamonds are the primary source of information on mantle carbon, and the carbon in most diamonds shows a remarkably consistent isotopic composition through time, that is, δ13C ~ −5 ± 1‰ (, ) , where R = 13C/12C expressed relative to the Vienna Pee Dee Belemnite (VPDB) standard]. This isotopic composition is indistinguishable from that of carbon in most mantle-derived magmas including mid-ocean ridge basalts (MORBs) (–), most CO2-rich kimberlites (, ) and carbonatites (), thus providing a robust baseline for the isotopic signature of mantle carbon. In contrast, the carbon-isotope composition of diamonds containing eclogitic inclusions, i.e., remnants of crustal material subducted into the mantle, is skewed toward lower δ13C values (, ). Organic matter is characterized by distinctly lower δ13C of ~−20 to −30‰ () compared to the ambient mantle (−5 ± 1‰) and marine carbonates (presently +2 ± 2‰) (). The combination of low δ13C values in diamonds with higher-than-mantle , where R = 18O/16O expressed relative to the Vienna Standard Mean Ocean Water (VSMOW) standard] values in their silicate mineral inclusions (, ) indicates recycling of subducted material containing organic carbon into the mantle. This process appears to have occurred episodically throughout most of Earth’s history, with the oldest diamonds that contain eclogitic inclusions forming at ~3 billion years (Ga) ().

Diamonds provide a limited view of the deep carbon cycle because the vast majority are sourced from the subcontinental lithospheric mantle (), which only extends to depths of ~200 to 250 km. Rare sublithospheric diamonds, which are limited to a few localities worldwide, do indicate that surface-derived, isotopically light carbon of organic origin can reach the asthenosphere and mantle transition zone (). Ocean-island basalts (OIBs) could potentially provide a more comprehensive picture of the carbon-isotope variability in the sublithospheric mantle, as their origin is commonly linked to mantle plumes (i.e., solid-state upwellings), many of which are probably sourced in the lower mantle. The low δ13C values of some OIBs (, ) provide evidence for transport of subducted organic carbon into the deep mantle. However, CO2 degassing profoundly alters the carbon budget and isotopic composition of these low-CO2 magmas (, ) because CO2 has low solubility in silicate magmas at crustal pressure. This issue, together with the limited temporal coverage of OIBs [≤150 million years (Ma)], suggests that an alternative approach is required to examine the evolution of the deep (sublithospheric) carbon cycle through time.

To address this knowledge gap, we have assessed the carbon-isotope compositions of hypabyssal (i.e., subvolcanic) kimberlites and related carbonate-rich ultramafic lamprophyres (i.e., aillikites) from localities worldwide with ages between ~2060 and 0.012 Ma (fig. S1). For simplicity and because the results are dominated by kimberlites, we collectively refer to both rock types as kimberlites. All samples used in this study are carbonate-rich magmatic rocks derived from low-degree partial melting of the sublithospheric (convective) mantle beneath thick continental regions (–). An association of most kimberlites with plumes from the lower mantle is supported by the geographic correspondence between Phanerozoic kimberlites at the time of their eruption and seismically anomalous regions (large low shear wave velocity provinces) above the core-mantle boundary, where most mantle plumes probably originate (). The genetic link between kimberlites and deep mantle plumes is reinforced by the distribution of some kimberlite fields along age-progressive corridors corresponding to the continental portions of hot spot tracks (). Furthermore, the relatively homogeneous Nd and Hf isotope composition and the occurrence of negative 182W anomalies in kimberlites older than ~200 Ma requires tapping of a deep source largely isolated from mantle convection and associated recycled crustal components (, ), and hence probably located in the lowermost mantle (). This is supported by the entrainment of superdeep diamonds, which contain inclusions of minerals stable in the mantle transition zone and lower mantle [e.g., ringwoodite () and CaSi-perovskite ()], and the He and Ne isotope systematics of olivine in some kimberlites, which require deep-mantle plume contributions ().

Although kimberlites represent complex mixtures of magmatic phases, mantle and crustal xenocrysts, and hydrothermal components (), their enrichment in magmatic carbonates makes bulk-carbonate carbon-isotope analyses of fresh samples representative of their magmatic values. This was recently demonstrated by comparing the results of in situ analyses of different textural types of carbonates by secondary-ion mass spectrometry (SIMS) with conventional bulk-sample measurements for the same kimberlites (see Materials and Methods) (). As carbon isotopes do not significantly fractionate between source and melt during partial melting of the mantle owing to the high incompatibility of carbon (), and are only marginally affected by kimberlite melt degassing (see below), bulk-carbonate carbon-isotope analyses can be used to trace the composition of the deep-mantle sources of kimberlites through time.

RESULTS

We analyzed 161 samples from 69 localities for their bulk CO2 concentrations and bulk-carbonate carbon and oxygen-isotope compositions and compiled existing data (tables S1 and S2). To minimize the effects of isotopic heterogeneity because of contributions from different textural types of carbonate in each sample (, , ), we calculated the average carbon- and oxygen-isotope composition of bulk-carbonate samples for each pipe or, where insufficient data were available, each cluster of kimberlites (see Materials and Methods). The samples analyzed in this and previous studies are associated with reliable emplacement ages, and for the few kimberlites that have not been dated, the age was estimated on the basis of neighboring kimberlites. When the δ13C values of kimberlites are plotted as a function of time (Fig. 1), a remarkable pattern emerges: All (except Pipes 10 and 14, Finland) examined kimberlites older than 350 Ma (n = 32) exhibit carbon-isotope compositions within the mantle range (δ13C ~ −4 to −6‰). Conversely, 8 of the 32 kimberlites younger than 250 Ma show δ13C (± 1σ) values lower than −6‰. Excluding Pipes 10 and 14, the weighted mean of δ13C for the >350-Ma kimberlites is −5.0 ± 0.6‰, while for the <250-Ma kimberlites the mean value is −6.1 ± 1.2‰. Kolmogorov-Smirnov tests indicate that both kimberlite populations are normally distributed, and a two-sample t test confirms that the difference in mean δ13C of the two populations does not occur by chance (at 95% confidence level, P < 0.05). These results demonstrate that younger kimberlites have, on average, lower δ13C values and a larger spread in carbon-isotope compositions than older kimberlites.

Fig. 1.

Carbon-isotope compositions of kimberlites and aillikites through time.

Each data point represents the average of multiple analyses of ≥3 samples, and error bars indicate the SD of the mean (see table S2). The green bar shows the carbon-isotope composition of the ambient mantle (δ13C = −5 ± 1‰).

DISCUSSION

Kimberlites younger than 250 Ma and with carbon-isotope excursions below the mantle range occur in southern Africa, western and eastern Canada, Brazil, and South Australia (Fig. 1). This isotopic signature is therefore unlikely to stem from local phenomena and implies a global process—although biases related to incomplete preservation of kimberlites in the geological record and sampling limited to fresh hypabyssal rocks cannot be completely ruled out. Exsolution (or degassing) of CO2-rich fluids during kimberlite magma ascent can fractionate carbon isotopes in the residual melt toward lower 13C/12C, whereas oxygen isotope variations are negligible (Fig. 2A) because of the high abundance of oxygen in the melt phase (see Materials and Methods). A degassing process could therefore explain the low δ13C values observed in some kimberlites containing less than 10 weight % (wt %) CO2, as shown in Fig. 2B. In addition, kimberlite samples from some localities (Orapa, Renard; fig. S4) show a direct correlation between δ13C and bulk CO2 concentrations, which is qualitatively consistent with carbon-isotope fractionation because of CO2 exsolution. However, multiple lines of evidence exclude a prominent role of CO2 loss in the carbon-isotope systematics of most kimberlites. (i) Kimberlites from most areas show limited variability in δ13C values that are not correlated to CO2 contents—including localities where kimberlites show a large spread in CO2 contents potentially related to variable CO2 degassing or fluid exsolution (fig. S4). (ii) Most kimberlites with low CO2 contents (e.g., <5 wt %) that might be attributed to CO2 degassing do not show appreciable deviations from carbon-isotope values typical of the mantle (n = 12 of 17; table S2). Lack of correlation between carbon-isotope ratios and CO2 concentrations in kimberlites worldwide probably stems from the fact that CO2 contents in bulk kimberlite samples are controlled by several processes including primary melt composition and CO2 loss as well as variable entrainment and assimilation of mantle and crustal debris, magmatic differentiation, hydrothermal alteration, and attendant carbonate replacement (, , , –). (iii) The broad correlation observed between carbon and hafnium or strontium isotopes (Fig. 3 and fig. S3) essentially excludes CO2 exsolution as a cause of the low δ13C values, because Hf and Sr do not quantitatively partition from melt into exsolved fluid or gas phases. (iv) Attribution of low δ13C values in younger (<250 Ma) kimberlites to CO2 loss is contrary to the absence of light carbon-isotope compositions in older kimberlites, which share similar bulk compositions and experienced similar ascent and emplacement processes, including CO2 loss by degassing/fluid exsolution as recently documented (). In summary, CO2 loss can generate localized low carbon-isotope values in kimberlite magmas, but there is no evidence that this process operated in all the low-δ13C kimberlites documented in this study.

Fig. 2.

Bulk-carbonate carbon and oxygen-isotope composition and bulk-rock CO2 concentrations in kimberlites and aillikites worldwide.

Each data point represents the average of multiple analyses, and error bars indicate the SD of the mean (see table S2). The green bars show the isotopic composition of the ambient mantle (δ13C = −5 ± 1‰; δ18Ocarbonate = 7.5 ± 1.5‰). The blue curve shows how the composition of kimberlite melts changes with increasing loss of CO2 at high temperature (≥1200°C) under open system conditions (Rayleigh fractionation) assuming that the primary kimberlite melt has a mantle-like carbon and oxygen isotope composition and a CO2 content of 20 wt % (, , ); each horizontal dash represents a 1 wt % increment of lost CO2. The gray dots represent the results of Monte Carlo simulations assuming δ13C = −4 to −6‰ in the mantle source region of kimberlites and variable CO2 contents between 15 and 25 wt % in primary kimberlite melts (see Materials and Methods for details). Note the similar spread in δ18O values of pre–250-Ma kimberlites (square symbols, warm colors) and post–250-Ma kimberlites (circles, cold colors) compared to the larger spread in δ13C toward lower values shown by post–250-Ma kimberlites.

Fig. 3.

Comparison of carbon and hafnium isotope compositions of kimberlites and aillikites worldwide.

εHfi represents the deviation of age-corrected 176Hf/177Hf from CHUR (chondrite uniform reservoir) at the time of kimberlite emplacement. Each data point represents the average of multiple analyses, and error bars indicate the SD of the mean (see table S2). The dotted blue circles indicate Pipes 10 and 14 (Finland), the only kimberlites older than 250 Ma, which show δ13C values lower than the mantle range. The green bar shows the carbon-isotope composition of the ambient mantle (δ13C = −5 ± 1‰). The red dashed curves represent mixing lines between the depleted kimberlite-source mantle and partially devolatilized, subducted sediments at ~150 Ma; the red numbers next to each curve indicate the carbon concentration in the sediments as μg/g, while the blue values next to the horizontal red dashes show the relative amounts of sediments in the source. The gray dots are the results of Monte Carlo simulations of the mixing calculations assuming 500 ± 200 μg/g of carbon in the recycled sedimentary component. Note the broad yet statistically robust direct correlation between δ13C and εHfi (R2 = 0.35). Details of calculations and statistical methods are reported in Materials and Methods.

Crustal contamination by CO2-rich fluids sourced from low-δ13C shales or similar wall-rock lithologies could also lower bulk-carbonate 13C/12C while simultaneously increasing δ18O values to above the mantle range (). The <250-Ma kimberlites with δ13C lower than the mantle, however, exhibit both mantle-like δ18O (i.e., Lac de Gras) and very heavy δ18O values (i.e., South Australia; Fig. 2A), and there is no statistically significant correlation between δ13C and δ18O in these kimberlites (see Materials and Methods). These observations suggest that crustal contamination mediated by high-δ18O fluids could have an impact locally but does not control bulk-carbonate δ13C compositions in kimberlites globally. Some kimberlites in both age groups were emplaced in crustal sequences containing shales or their metamorphic equivalents (e.g., the ~1.1-Ga Premier kimberlite), yet low-δ13C kimberlites are essentially limited to the last 250 Ma. Assimilation of metasomatized lithospheric mantle hosting low-δ13C graphite, diamond, carbides, or carbonates (, , ) could decrease the 13C/12C of kimberlite magmas. However, given the antiquity (mostly >1 Ga and occasionally >3 Ga) of low-δ13C lithospheric diamonds (, ), it is unlikely that low-δ13C lithospheric mantle material was quantitatively assimilated only in kimberlites younger than 250 Ma—unless mantle metasomatism introduced additional isotopically light carbon in the lithospheric mantle during the late Phanerozoic.

Recycling of subducted crustal material containing low-δ13C carbon of organic origin into the deep mantle source of some kimberlites younger than 250 Ma arguably represents the most plausible process to explain the low δ13C compositions. This is consistent with the Sr-Nd-Hf isotope compositions of <200-Ma kimberlites from southern Africa, Brazil, western Canada, and South Australia, which have more geochemically enriched signatures (i.e., lower initial 143Nd/144Nd and 176Hf/177Hf, and higher initial 87Sr/86Sr) compared to coeval kimberlites elsewhere (fig. S2). These characteristics have been attributed to contributions from deeply subducted crustal material (). Statistically meaningful correlations are observed between δ13C, bulk-sample initial 176Hf/177Hf, and perovskite initial 87Sr/86Sr (Fig. 3 and fig. S3; see Materials and Methods). The isotopic correlations vary from time-integrated, moderately depleted compositions characterized by suprachondritic 176Hf/177Hf, low 87Sr/86Sr (~0.703), and δ13C approaching −4‰ [i.e., similar to carbon-isotope ratios in global MORBs ()] to geochemically enriched compositions of likely subducted sedimentary origin, i.e., subchondritic 176Hf/177Hf, moderately high 87Sr/86Sr (~0.705), and δ13C as low as ~−8‰ (Fig. 3). Mixing models and mass-balance calculations indicate that incorporation of up to ~10 to 15% of metamorphosed and partially devolatilized marine sediments containing <1000 μg/g of organic carbon into the mantle source of kimberlites can generate the correlated carbon and Hf isotope compositions observed in <250-Ma kimberlites (Fig. 3 and Materials and Methods). These results are broadly consistent with independent modeling of the highest extent of recycled material contribution in the sources of Cretaceous kimberlites from southern Africa, Brazil, and western Canada based on the oxygen-isotope compositions of olivine ().

The shift toward lower carbon-isotope compositions in kimberlites coincides with a remarkable increase in the frequency of kimberlite eruptions after 250 Ma ago (Fig. 4). Increased kimberlite activity in the Phanerozoic was previously linked to the onset of colder subduction regimes in the Neoproterozoic (), which is broadly consistent with a recent statistical analysis of the pressure and temperature conditions of metamorphic rocks globally (). A global change in the thermal regime of subduction zones alone, however, cannot explain the combination of increased kimberlite activity and carbon-isotope perturbation, which also requires an increased contribution from deeply subducted organic carbon. The flux of carbon into subduction zones is presently dominated by carbonates in sediments and altered oceanic crust (), with organic carbon representing less than 20% of the total carbon input (). On the other hand, heavy carbon-isotope compositions in arc magmas and related gases combined with thermodynamic modeling of carbonate stability at the pressure and temperature conditions of subducted plates indicate that carbonates are efficiently stripped from subducted slabs at fore-arc and sub-arc depths (, , , ). In addition, partial melting experiments at high pressure and temperature indicate that the mantle transition zone represents an efficient thermal barrier to deeper carbonate subduction (). Conversely, the low solubility of reduced organic carbon (i.e., graphite and diamond) at high pressure in hydrous fluids () and slab melts () limits the extraction of organic carbon from subducted sediments. Therefore, the combination of carbon-isotope data for kimberlites with geochemical, theoretical, and experimental modeling of the behavior of oxidized and reduced carbon in subducted lithologies (, , , –) supports a higher flux of sedimentary organic (reduced) carbon to the lower mantle in the Phanerozoic. The variable δ13C values in kimberlites younger than 250 Ma might reflect global variability in this influx as well as the heterogeneous distribution of recycled organic carbon in the deep mantle source of kimberlites. An important question remains—what Earth surface processes could have generated such a change?

Fig. 4.

Comparison of kimberlite carbon-isotope compositions, frequency of kimberlite events, and total organic carbon contents in shales in the last 1000 Ma.

(A) Kimberlite and aillikite δ13C versus time (taken from Fig. 1). (B) Number of kimberlite events every 50 Ma [compilation of ()]. (C) Boxplot showing the variability of total organic carbon (TOC) in shales binned into 50-Ma intervals [compilation of ()]. Each box is delimited by the first and third quartile, and the horizontal line represents the median value. The ages of major events of geological significance for this work are highlighted, i.e., start of present-day cold subduction after ~850 Ma (), the Cambrian explosion at ~542 Ma, and the kimberlite “bloom” at ~250 Ma. The horizontal black bar shows the time (≥260 Ma) required for subducted material to return to the surface via magmatism related to deep-mantle plumes (). Entrainment of sedimentary material containing isotopically light organic carbon, which was subducted after the Cambrian Explosion, is consistent with a wider spread in kimberlite δ13C values and increased frequency of kimberlite events after 250 Ma.

The traditional view, based on the long-term stability of the carbon-isotope record of marine carbonates, was that the extent of organic carbon burial was relatively constant on 100-Ma time scales (). More recent work, however, has questioned this conclusion. Specifically, there is increasing appreciation that the extent of global organic carbon burial based on the carbon-isotope record of marine carbonates is dependent on surface oxygen levels (, ). With lower atmospheric oxygen levels in the Precambrian (), oxidation of sedimentary organic carbon would have been limited, questioning one of the fundamental assumptions of the traditional view of the global carbon-isotope mass balance. In this framework, it is possible to have a significant increase in organic carbon burial—linked to an increase in marine primary productivity—without a major shift in the marine carbon-isotope record (, ). Furthermore, there is increasing evidence that the most likely means of stabilizing a low-oxygen system is to have strong nutrient limitation and reduced marine primary productivity, which leads to less burial of organic carbon even with anoxic oceans (, ). Last, a recent compilation of organic carbon content in the sedimentary record suggests that there was a significant increase in the deposition of organic carbon at the Proterozoic-Phanerozoic boundary (Fig. 4) (). Critically, this compilation, in contrast to most sedimentary-rock geochemical databases, focuses on sections with stratigraphically continuous sampling, minimizing the effect of preferential analysis of organic-rich units and providing a more representative record of the sedimentary organic carbon record. Although a large proportion of organic matter is currently deposited on shelves in oxic oceans [~50% at present-day conditions; ()], anoxic deep marine conditions would have favored more efficient burial, preservation, and potentially subduction of organic carbon—e.g., (). In summary, available constraints from the sedimentary record are consistent with a shift in organic carbon burial roughly coincident with the Precambrian-Cambrian boundary.

While increased deposition and subduction of organic carbon following the Cambrian Explosion could introduce isotopically light carbon into the mantle, there is a ≥300-Ma lag between the onset of the Cambrian Explosion and changes in the carbon-isotope geochemistry of kimberlites after ~250 Ma. Geodynamic modeling shows that ≥250 to 300 Ma are required for subducted plates to reach the core-mantle boundary and return to Earth’s surface via plume-related magmatism (). This can reconcile the temporal gap between biogeochemical changes at Earth’s surface and perturbation of the deep carbon cycle recorded by kimberlites and might suggest a short residence time for subducted crustal material in the convective mantle, at least locally. For example, the oldest Phanerozoic kimberlite with δ13C values significantly lower than the mantle range (Jacare in Brazil; δ13C = −7.3 ± 0.1‰) has an age of ~242 Ma (table S2), which requires a minimum subduction age of surface-derived organic carbon of ~500 Ma, i.e., ~40 Ma after the Cambrian explosion. Development of land plants after ~450 Ma and their later rise in the Devonian () could have provided a further boost in the delivery of organic carbon to marine basins, but there is currently no strong evidence from sedimentary records for an increase in total organic carbon contents associated with the rise of plants (). If the kimberlite source is located in the upper mantle as advocated by some (), the carbon-isotope perturbation recorded by kimberlites still requires a fundamental change in their mantle sources during the Paleozoic because of subduction of organic carbon, and, hence, a likely (though more loosely constrained) link to the Cambrian Explosion. Regardless, within this framework, the carbon-isotope record of kimberlites provides independent support for a significant increase in the extent of primary productivity broadly related to the Cambrian Explosion.

In summary, episodic transport of reduced organic carbon into the mantle is demonstrated by the carbon-isotope systematics of eclogitic diamonds as old as 3 Ga () and, potentially, some Precambrian kimberlites (Pipes 10 and 14, Finland; Fig. 1). During the Phanerozoic, however, colder subduction zones and increased rates of organic carbon burial associated with a more productive marine biosphere (, ) promoted a radical change characterized by a higher flux of organic carbon into the mantle that perturbed the deep carbon cycle. This is reflected in the abundance and carbon-isotope values of kimberlites younger than 250 Ma. The higher frequency of kimberlite and carbonatite eruptions after 250 Ma [this work; ()] also suggests that more carbon might have been released from the deep mantle via volcanism after this time. This hypothesis should be tested by exploring the temporal distribution of more widespread silica-undersaturated mafic magmas (e.g., alkali basalts, nephelinites, and basanites), which are sourced in mantle domains that may contain deeply subducted carbon. This work demonstrates that kimberlites provide a previously overlooked perspective on the deep-mantle cycling of surface-derived material and the links between the deep and shallow carbon cycles including their Phanerozoic biogeochemical perturbations. It also provides a novel perspective on the evolution of the global carbon cycle and supports the emerging view that the extent of marine primary productivity has markedly changed through Earth’s history.

MATERIALS AND METHODS

Carbonates in kimberlites

Kimberlites are hybrid rocks containing components of magmatic, xenocrystic, and hydrothermal origin (, , ). Different textural and mineralogical types of carbonates occur in kimberlites: microphenocrysts, interstitial groundmass grains, segregations with or without serpentine, pseudomorphs after olivine or other silicates, veins of variable size (up to centimeters wide), and rarer country-rock xenocrysts (, , , , ). While the dominant carbonate is calcite, with less abundant to scarce dolomite, different carbonate textures can be associated with variations in oxygen, strontium, and carbon isotopes (, , , , , ). Calcite microphenocrysts, interstitial groundmass grains, and serpentine-free carbonate segregations are commonly considered to be magmatic based on their textures, elevated concentrations of Sr, Ba, and LREE (light rare earth elements), unradiogenic Sr isotope ratios, and mantle-like carbon and oxygen isotope compositions (, , , , –). Carbonate segregations containing serpentine and veins crosscutting the igneous texture are interpreted as hydrothermal based on higher 87Sr/86Sr than associated perovskite (, ). Therefore, bulk-carbonate analyses of carbon and oxygen isotopes can potentially provide mixed signals because of the mixture of these components ().

Previous work (, , ), however, has demonstrated that bulk-carbonate δ13C compositions of fresh kimberlites are largely indistinguishable from those of their magmatic carbonates, and therefore, bulk-carbonate analyses can be used to trace the carbon-isotope composition of kimberlite magmas and their sources. In addition, modeling of the effects of hydrothermal fluids on the carbon-isotope composition of carbonates in kimberlites () shows that carbonate δ13C values are minimally affected (~≤1‰) at the typical conditions [i.e., T ~ 200° to 400°C ()] of hydrothermal alteration. Therefore, bulk-carbonate analyses of kimberlites provide robust estimates of kimberlite melt δ13C even for mildly altered samples, where olivine and part of the groundmass are serpentinized. Conversely, hydrothermal alteration can produce bulk-carbonate δ18O values, which deviate considerably from mantle values (e.g., δ18Ocarbonate = +6 to +9‰), typically to higher values (, ). This is commonly observed in kimberlites worldwide as noted in this (Fig. 2) and previous studies (, , , ). Therefore, in this work, we focus on bulk-carbonate carbon-isotope compositions of hypabyssal (i.e., subvolcanic) kimberlites.

Bulk-carbonate carbon and oxygen-isotope analyses

We report the results of 231 bulk-carbonate carbon and oxygen isotope analyses of 144 kimberlite and aillikite samples from 60 localities (table S1). Most samples are either fresh or only mildly altered as defined above, but we have included some altered samples to assess the effects of extensive hydrothermal and supergene alteration. Clean chips of each sample were finely crushed in an agate mortar, and up to three different splits were analyzed separately to test for isotopic heterogeneity, which was found to be negligible. Subsamples were placed in glass vials sealed with septum caps and then placed on a temperature-controlled block set to 70°C, where the ambient air was purged using ultrahigh-purity helium. The powders were then reacted for 30 min with ~0.05 to ~0.3 ml of orthophosphoric acid, depending on sample mass. This is essentially a bulk-carbonate extraction because both calcite and any dolomite are fully dissolved. Subsample mass varied according to the (CO3)2− content of the bulk samples such that ~0.001 mol of CO2 gas from the acid-sample reaction was available for measurement. CO2 was measured using an Analytical Precision AP2003 continuous-flow stable isotope ratio mass spectrometer at the University of Melbourne. The sample gas 13C/12C and 18O/16O were converted to the conventional “delta notation” (δ13C and δ18O) and normalized to the VPDB and VSMOW scale, respectively, using two in-house reference standards (NEW1 and NEW12) previously calibrated against the international reference standards NBS18 and NBS19 (table S3) (). The 1σ reproducibilities for carbon and oxygen isotope compositions of individual measurements were 0.07 and 0.08‰, respectively, based on replicate analyses of the NEW1 standard (n = 98). A very minor scale correction was applied to the oxygen-isotope data based on mean measured versus known values of the three standards—NEW1 (n = 98), NEW12 (n = 41), and NBS18 (n = 18). No correction was necessary for the carbon-isotope data.

Determinations of bulk-sample CO2 concentrations

Bulk-rock CO2 contents were measured in most (n = 79) of the samples analyzed for their bulk-carbonate carbon and oxygen isotope compositions (table S1). For two subsets of samples (those from India and Finland), bulk-rock CO2 contents had previously been reported (table S1) (–). CO2 concentrations in finely ground powders were measured by infrared (IR) spectroscopy at ETH Zurich using a LECO CS844 carbon and sulfur analyzer. Approximately 1 g of copper accelerator and 10 mg of sample powder were placed in a ceramic crucible and combusted in a stream of purified oxygen. CO2 concentrations were determined by two nondispersive IR cells calibrated relative to LECO reference material 502-950 (synthetic carbon, C = 0.12 ± 0.01 wt %). Average CO2 concentrations in Fig. 2B were calculated by combining these new results and previously published data (table S2).

Compilation and assessment of carbon-isotope data

To minimize the effects of isotopic heterogeneity due to contributions from different textural types of carbonate in each sample (, , ), we have averaged the isotopic composition of multiple samples from the same kimberlite pipe or, where insufficient data were available, cluster or field of kimberlites. For each locality, at least three samples were available to calculate an average δ13C value. The exceptions are Jacare (Brazil), Frank Smith (South Africa), and Mir (Siberia) for which only two samples were available (table S2). We excluded samples that are highly weathered based on an inspection of hand samples and thin sections, because these can preserve δ13C compositions distinct from those of fresher samples from the same pipe or area. For example, two altered samples from the Alto Paranaiba field (Brazil), selected to test the impact of alteration, exhibit δ13C of −0.8‰ (Successo-08) and +0.4‰ (Santa Clara-01), respectively, compared to a mean of −7.7 ± 1.6‰ for the fresher samples (table S1). The samples from some regions (e.g., South Australia) exhibit very high bulk-carbonate δ18O (Fig. 2), associated with δ13C below the mantle range. These elevated δ18O values could indicate overprinting of the carbonate composition by crustal fluids. These samples are, however, fresh or only mildly altered, and they contain abundant interstitial carbonates of apparently magmatic origin with carbonate replacing olivine in the South Australian aillikites. We have opted to retain these results; rejecting them would not alter our conclusions. We also note that there is no statistically meaningful correlation between δ13C and δ18O values in <250-Ma kimberlites (R2 = 0.21; Student’s t test: tcalc = 2.8 ~ tcrit (0.01; 32) = 2.8), which further excludes a role of crustal contamination in the low carbon-isotope compositions documented in these kimberlites.

New analyses of volcaniclastic kimberlites from Alto Paranaiba (Brazil), Somerset Island (Canada), and Orapa (Botswana; table S1) show δ13C values that can differ substantially from hypabyssal kimberlites from the same area (e.g., Somerset Island: −3.1 ± 1.1‰ for volcaniclastic kimberlites versus −5.8 ± 0.3‰ for the hypabyssal samples; Orapa: −3.8 ± 2.9‰ versus −6.3 ± 1.2‰). These higher δ13C values probably relate to the effects of alteration and entrainment of country rock material rather than degassing, which can only lower δ13C (see the “CO2 loss model” section), because volcaniclastic kimberlites are easily altered and generally contain abundant crustal xenocrysts. Four previously studied (), apparently coherent but pyroclastic kimberlites () from the Lac de Gras field (Canada) were included in this study because they are fresh and have compositions similar to those of coeval hypabyssal kimberlites, i.e., they experienced mild CO2 degassing that did not affect their carbon-isotope compositions. None of these pyroclastic kimberlites shows δ13C values below the mantle range. We also avoided magmas emplaced into limestone sequences to reduce the effects of crustal contamination. The only exception is the Udachnaya-East kimberlite (Siberia), which shows robust evidence of limestone contribution based on a direct correlation between bulk-carbonate δ13C (−5.9 to −1.6‰) and bulk-sample CO2 concentrations (fig. S4). The average δ13C (−3.5 ± 0.8‰; n = 49) of this kimberlite, however, overlaps the upper end of the mantle range.

The mean δ13C values reported in table S2 and used throughout this work do not include outliers, defined as values beyond ±3σ of the arithmetic mean. An independent test of our approach is provided by the coeval kimberlites and olivine-lamproites in the Wajrakarur field (India), which show indistinguishable Sr, Nd, and Hf isotope compositions (, ); our data extend this similarity to carbon isotopes (kimberlite δ13C = −4.7 ± 1.6‰, n = 4; olivine-lamproite δ13C = −4.9 ± 1.7‰, n = 10; table S2). The oxygen-isotope dataset was not screened in this way because the δ18O of carbonates can be modified by interaction with fluids of deuteric (i.e., magmatic) and/or crustal origin due to fast diffusion of oxygen in carbonates under hydrous conditions, while leaving carbon isotopes largely unaffected (, , ), and the focus of this work is on carbon isotopes.

Carbon isotopes versus Sr-Nd-Hf isotopes

To compare the carbon-isotope composition of kimberlites with other geochemical tracers of kimberlite source evolution through time, we calculated average Sr, Nd, and Hf isotope compositions for kimberlites from the same pipe, cluster, or field using the recent compilation of Giuliani et al. (). For Nd and Hf isotope ratios, we used bulk-kimberlite analyses because these magmas are so enriched in incompatible trace elements () that crustal contamination generally has minimal effect on Nd and Hf isotope ratios of fresh samples. For Sr isotopes, we only used analyses of perovskite, a magmatic phase that is minimally affected by alteration and crustal contamination (). When insufficient radiogenic isotope data were available (i.e., <2 bulk-rock Nd-Hf isotope results; no perovskite Sr isotope data), we used Sr-Nd-Hf isotope data for the entire cluster or field of kimberlites (see table S2 for details). This approach is validated by the limited spread in radiogenic isotopic compositions of kimberlites within the same cluster or field, i.e., very low standard deviations associated with average initial 143Nd/144Nd, 176Hf/177Hf, and 87Sr/86Sr. There are notable exceptions, e.g., the Lac de Gras (Canada) and Kaavi-Kuopio (Finland) fields (, ), which are not considered here. Even adopting this approach, for some kimberlites, only one analysis of perovskite 87Sr/86Sr was available (table S2) to examine correlations between carbon and Sr-Nd-Hf isotopes. It is important to note that this comparison does not require restriction to samples for which both carbon and radiogenic isotopes have been measured in the same sample, because averaging the isotopic composition of multiple samples from the same pipe or cluster of kimberlites yields representative isotopic compositions for the kimberlite mantle source.

Statistically significant correlations are observed between bulk-carbonate δ13C and both initial εHf (R2 = 0.35, n = 43; Fig. 3) and initial 87Sr/86Sr (R2 = 0.34, n = 27; fig. S3). εHf is defined as the deviation of age-corrected 176Hf/177Hf from CHUR (chondrite uniform reservoir) at the time of kimberlite emplacement with CHUR values of Bouvier et al. (). Standard two-tailed Student’s t tests show that none of the correlations examined is likely to have occurred by chance (99% confidence level) because t values are significantly higher than critical values, i.e., C versus Hf isotopes: tcalc = 4.7 > tcrit (0.01; 43) = 2.7; C versus Sr isotopes: tcalc = 3.6 > tcrit (0.01; 27) = 2.8. Conversely, there is no significant correlation between bulk-carbonate δ13C and initial εNd (R2 = 0.11, n = 48, tcalc = 2.3 < tcrit (0.01; 48) = 2.7; fig. S3). This last observation is surprising given the well-established correlation between Nd and Hf isotopes in terrestrial rocks, including mantle-derived magmas (). The lack of correlation between bulk-carbonate δ13C and initial εNd, however, largely stems from the anomalously low initial 143Nd/144Nd of the Lac de Gras kimberlites compared to archetypal kimberlites worldwide ().

The correlations between δ13C, 176Hf/177Hf, and 87Sr/86Sr indicate variable mixing of moderately trace-element depleted and geochemically enriched components in the source of kimberlites. The moderately depleted endmember, which is characterized by a suprachondritic 176Hf/177Hf and mantle-like carbon-isotope composition, might represent the depleted component that ubiquitously occurs in kimberlites worldwide (). Our approach employs the Sr-isotope compositions of perovskite, which are insensitive to crustal contamination (). Hence, the geochemically enriched component(s) could be either assimilated lithospheric mantle or deeply subducted sedimentary material. We reject the hypothesis that the enriched component represents basaltic oceanic crust because >2 Ga of aging is required to develop the required low εHf values (, ), at odds with an origin of this component via Phanerozoic subduction.

Previous studies of southern African carbonate-bearing olivine-lamproites (previously known as orangeites or Group II kimberlites) (), mantle xenoliths (–), and eclogitic diamonds (, ) have shown that metasomatism can introduce isotopically light carbon in the lithospheric mantle. The compositions of the Wajrakarur olivine-lamproites (India) and Torngat aillikites (Canada), however, indicate that this process is not ubiquitous. The Wajrakarur olivine-lamproites (India) have carbon-isotope compositions within the mantle range and are indistinguishable from those of coeval kimberlites, which are dominated by sublithospheric components (table S2). An extensive contribution by enriched lithospheric mantle was identified in the source of the Torngat aillikites based on 143Nd/144Nd and 176Hf/177Hf extending from moderately suprachondritic to subchondritic values (). Despite the isotopic variability, these aillikites show mantle-like δ13C of −5.3 ± 0.8‰ (n = 27; table S2). If kimberlites and aillikites younger than 250 Ma inherited their low-δ13C component from assimilated lithospheric mantle material, this would require metasomatic addition of low-δ13C material during the Phanerozoic. While this process cannot be completely ruled out, it is not favored here because surface-derived, high-δ13C oxidized carbon can be efficiently recycled into the upper mantle as shown by some diamonds () and mantle xenoliths (). Conversely, the transition zone provides an effective thermal barrier to the subduction of oxidized carbon into the lower mantle (), from where most kimberlites are probably sourced (, , ).

Mixing calculations of sediment recycling in the mantle

To address the effects of deep recycling of sedimentary material in the source of kimberlites, we performed mixing calculations using the depleted kimberlite mantle source of Giuliani et al. () and partially devolatilized marine sediments as endmember components. We only modeled the relation between carbon and Hf isotopes in kimberlites due to the large uncertainties associated with Rb/Sr fractionation in marine sediments during subduction due to the effects of high-pressure metamorphism and partial melting. In this model, the carbon concentration of the mantle component is 117 μg/g, which Bekaert et al. () consider representative of the source of lower mantle plumes. Its δ13C is −4.5‰, which is in the carbon-isotope compositional range of kimberlites emplaced before 350 Ma (δ13C ~ −5.0 ± 0.6‰) and intermediate between kimberlites and the MORB source [δ13C ~ −4.0‰ ()]. The depleted source of kimberlites has a Lu/Hf ratio () that is remarkably similar to that of the Early Depleted Reservoir (EDR) modeled by Jackson and Jellinek (). Hence, we use the EDR Hf concentration (0.24 μg/g) to constrain the composition of the kimberlite source before mixing, and estimate the Lu concentration (0.061 μg/g) by using the Lu/Hf ratio of the depleted kimberlite source (). The initial εHf value of the kimberlite source before mixing is calculated using the linear regression exhibited by kimberlites in 176Hf/177Hf versus time ().

The composition of the recycled sedimentary component is based on the GLOSS-II composition () assuming that Cambrian marine sediments had similar compositions to those in modern ocean basins. Subducted sediments contain, on average, 3.1 wt % of CO2 (), of which ~1/10 is of organic origin (). This corresponds to 837 μg/g of organic carbon. Organic carbon is partially lost during progressive graphitization associated with diagenesis and low-grade metamorphism, with estimates ranging between ≥10% and >90% of carbon loss (, ). On the basis of this low-pressure loss plus additional devolatilization during transit through the mantle, we varied the concentrations of organic carbon (Corganic) in the recycled component between 300 and 700 μg/g. During diagenesis and metamorphism, the isotopic composition of organic carbon increases from the canonical δ13C value of ~−25‰ due to methane loss () and/or high-temperature equilibration with coexisting carbonates (). Hence, for the purpose of modeling, a δ13C value of −15‰ was used, which is in the upper range of the carbon-isotope composition previously measured for high-grade metapelites (, ). The Hf content (5.13 μg/g) of the sedimentary component was assumed to be 50% higher than that of GLOSS-II, to account for partial devolatilization and loss of fluid-mobile elements. To estimate εHf at the time of mixing, we calculated the present-day value (−11.5), which is based on the Nd isotope composition of GLOSS-II and the well-established correlation between Nd and Hf isotopes in marine sediments (), and assumed no Lu/Hf fractionation in subducted sediments ().

We assumed a single mixing event at 150 Ma, approximately centered between the oldest (Jacare in Brazil, ~242 Ma) and youngest (Rattler in Canada, ~60 Ma) kimberlites with low δ13C compositions (table S2). Increasing the time of mixing to 250 Ma has a negligible effect on the results (e.g., εHfmix increases from −3.2 to −2.6 assuming 10% of sediments in the kimberlite source). Radiogenic ingrowth in the source after mixing generates minimal variations of εHf (e.g., ~1 unit for 10% of sediment contribution after 200 Ma) unless the source includes more than 20% of sediments, because the Lu/Hf of the source is similar to that of CHUR. To address the variability in mixing age (i.e., 150 ± 100 Ma), compositions of depleted kimberlite mantle source, and subducted sedimentary components, we used Monte Carlo simulations. Each model parameter was randomized assuming an arbitrary relative uncertainty of 20%, except for Corganic, which was set at 500 ± 200 μg/g (i.e., ± 40%), and blocks of 300 simulations were constructed for each parameter.

The results of two mixing calculations (i.e., Corganic = 300 and 700 μg/g) and these Monte Carlo simulations are shown in Fig. 3. These calculations indicate that the carbon and hafnium isotope compositions of most late Phanerozoic kimberlites can be reproduced by adding between 2 and 15% of partially devolatilized sediments to the depleted kimberlite source. Using lower concentrations of Corganic (e.g., 100 μg/g) in subducted sediments due to higher degrees of carbon loss during metamorphism and deep subduction matches the carbon-isotope compositions of the lowest-εHf kimberlites in Fig. 3, thus covering the whole compositional spectrum of kimberlites. These findings are remarkable given the wide variability in age and geographic settings of kimberlites and the likely compositional variations of subducted sediments in space and time. Support for this model also comes from previous sulfur and nitrogen-isotope studies of kimberlites and mantle xenoliths metasomatized by kimberlite melts shortly before transport to the surface, which point to a sedimentary input in the source of some Cretaceous southern African kimberlites (–).

CO2 loss model

The exsolution of a CO2-rich fluid phase may drive the very fast ascent of kimberlite magmas () and could also modify the carbon and oxygen isotope composition of kimberlites (, ). We have modeled the change in δ13C and δ18O after variable degrees of CO2 exsolution by assuming CO2 contents of 20 ± 5 wt % in the primary melt, which can be considered a robust estimate (–). The δ13C value of the residual melt after CO2 loss was calculated using a Rayleigh distillation model (i.e., open system conditions) and using a high-temperature (~1200° to 1400°C) carbon-isotope fractionation factor αCO2-melt of 1.0022 between CO2 and carbonate-bearing melt (). The δ18O of the residual melt was computed using a similar approach, but assuming a αCO2-melt factor of 1.00276, which describes oxygen isotope fractionation between CO2 and silicate-carbonate (kimberlite) melts at 1200°C (). The variability in carbon-isotope compositions in the mantle source (δ13C = −5 ± 1‰) and CO2 concentrations in primary kimberlite melts (20 ± 5 wt %) was assessed by using Monte Carlo simulations (blocks of 300 simulations for each randomized parameter). These calculations indicate that the δ13C value of kimberlites with mantle-like initial 13C/12C can decrease to values lower than the mantle range (i.e., <−6‰) if more than 7 to 9 wt % of CO2 is exsolved at high temperature during ascent (Fig. 2B). Instead, changes in δ18O are negligible (less than −0.8‰; Fig. 2A). If degassing occurred (at least partly) at lower temperatures or in a closed system, carbon-isotope fractionation would be higher, and hence, the modeled trajectory of the residual melt composition in δ13C versus bulk CO2 (Fig. 2B and fig. S4) would be steeper and even less consistent with the kimberlite data. Although this model does not consider the role of exsolved H2O, the results would not change markedly because H2O exsolution does not affect carbon-isotope fractionation and CO2 is the dominant volatile species in kimberlite melts (, ).