Melt inclusion constraints on mantle carbon heterogeneity within an individual mantle plume and at the global scale

The Earth’s mantle holds more carbon than its oceans, atmosphere and continents combined, yet the distribution of carbon within the mantle remains uncertain. Global variations in the carbon content of the depleted mantle have been inferred from carbon-trace element systematics of ultra-depleted mid-ocean ridge glasses and melt inclusions, though the origin of mantle carbon variability remains uncertain. A variety of observations suggest that the mantle plumes underlying ocean islands have high bulk carbon contents, but the characterisation of individual mantle plume components has not been possible. We supplement the existing melt inclusion data from Iceland with four new datasets, significantly enhancing the spatial and geochemical coverage. Within the combined dataset there is significant variation in melt inclusion CO2/Ba ratio, which is tightly correlated with trace element enrichment. The trends in CO2/Ba–Ba space displayed by our new data coincide with the same trends in data compiled from global ocean islands and mid-ocean ridges, forming a global array. The overall structure of the global CO2/Ba array is not a property of the source, instead controlled by CO2 degassing before melt inclusion entrapment, and CO2 loss by decrepitation following inclusion entrapment. The position of melt inclusion datasets within the global array allows us to filter out inclusions likely to have undergone decrepitation. Providing earlier CO2 degassing and magma mixing have not destroyed the mantle signal, the remaining inclusions may preserve mantle CO2/Ba values. Using the observed covariation of the filtered CO2/Ba ratio with radiogenic isotope tracers of mantle source, we find CO2 concentrations in the Icelandic depleted mantle of 102 −27 ppmw, typical of mantle underlying the mid-ocean ridge system, and high CO2 concentrations of 2.2 −1.5 wt% associated with the high He/He component of the Iceland plume. However, source and process are difficult to deconvolve for the eruptions most sensitive to enriched, or recycled, mantle components; the presently available data is consistent with recycled material being either enriched or depleted in carbon. Extending our approach globally, we find no systematic variation in CO2/Ba or CO2/Nb with enrichment (in the absence of a primordial mantle contribution). The lack of global co-variation suggests mechanisms recycling CO2 from the surface into the convecting mantle are decoupled from those recycling most trace elements. However, the variation of CO2/Nb with Ba/Nb suggests Ba may behave similarly to CO2 during recycling.


Introduction
The majority of Earth's carbon, that has not been sequestered to the core, resides in the mantle (Dasgupta and Hirschmann, 2010; the Earth's deep past, remain largely uncertain. A great deal of this uncertainty is associated with the efficiency of carbon recycling in subducting slabs; thermodynamic modelling suggests slab decarbonation should be limited in most subduction zones (Kerrick and Connolly, 2001), but the infiltration of hydrous fluids may mobilise carbonate from the subducting plate (Kelemen and Manning, 2015). Hirschmann (2018) placed empirical bounds on recycling efficiency by comparing the CO 2 /Ba ratios of the Earth's exosphere and upper mantle and found that at least 25%, and probably much more, of the carbon outgassed to the exosphere must have been returned to the mantle via subduction.
An alternative, but complementary, approach to constraining CO 2 recycling efficiency is to estimate the amount of carbon within recycled mantle domains where they are sampled at ocean islands. CO 2 -noble gas systematics of fluid inclusions in material erupted in Hawaii and La Réunion Island indicate their mantle sources are carbon rich (Trull et al., 1993;Boudoire et al., 2018), which is further supported by models that couple CO 2 emission rates with magma supply rates at Hawaii (Anderson and Poland, 2017) and high CO 2 contents of primitive melts (Tucker et al., 2019). However, both Hawaii and La Réunion are underlain by mantle plumes transporting primordial (high 3 He/ 4 He) mantle, which may contain high concentrations of carbon (Miller et al., 2019). In order to assess the efficiency of carbon recycling it is necessary to deconvolve the contributions of carbon from recycled and primordial mantle domains to volcanic outgassing.
Progress in identifying mantle carbon variability has been hampered by its tendency to exsolve from magmas into a CO 2 vapour phase. The solubility of CO 2 in basaltic melts is minimal at the low pressures of eruption (Stolper and Holloway, 1988), with the consequence that considerable CO 2 will have been lost by CO 2 vapour exsolution (degassing) in all but the most CO 2 poor melts.

Estimating mantle carbon
A number of approaches have been taken to circumvent the overprint of degassing on the mantle CO 2 signals. Rare ultra-depleted glasses erupted in oceanic fracture zones often contain sufficiently low concentrations of CO 2 that they are undersaturated in CO 2 vapour at the pressure of eruption. However, Matthews et al. (2017) suggested these magmas formed by mixing of trace-element and CO 2 enriched melts that are partially degassed with trace-element and CO 2 depleted melts; with the significant CO 2 vapour undersaturation of the depleted melts conferring CO 2 undersaturation onto the mixed magma. This process is difficult to definitively identify and could lead to underestimates of mantle CO 2 concentrations. Occasionally, extremely trace-element enriched mid-ocean ridge basalt (MORB) glasses retain gas rich vesicles, which may contain the full budget of CO 2 and other volatiles that exsolved from the magma (Cartigny et al., 2008;Jones et al., 2019). It is again unclear whether these enriched magmas lost CO 2 deeper within the crust, before beginning to retain vesicles. Furthermore, the MORB glass datasets only provide information about the mantle underlying the mid-ocean ridge system and are therefore biased towards characterising the depleted upper mantle.
Another approach uses melt inclusions, tiny droplets of magma trapped within crystals as they grow at depth within the crust. The elastic strength of the host crystals prevents decompression of the inclusion to the pressure of the surrounding melt, potentially allowing the melt to remain undersaturated in CO 2 vapour. Not only are melt inclusions ubiquitous in MORB, but they can also be readily found in ocean island basalt (OIB). Melt inclusions also offer an unrivalled view of the earliest stages of magma evolution, when the diversity of melt chemistry generated by near-fractional melting is still partially preserved (Sobolev and Shimizu, 1993). However, the CO 2 budget of melt inclusions remains susceptible to the same degassing and mixing processes as the ultra-depleted glasses. Postentrapment CO 2 loss via decrepitation of the host crystals is also likely to occur (Maclennan, 2017).
1.2. Ratios of CO 2 and trace-element concentrations One process affecting carbon that it may be possible to successfully control for is melting. Ratios of CO 2 to Ba or Nb concentrations have been used to achieve this (Saal et al., 2002;Le Voyer et al., 2017;Hauri et al., 2018;Michael and Graham, 2015;Shimizu et al., 2016;Tucker et al., 2019;Shimizu et al., 2019;Miller et al., 2019), since CO 2 , Ba, and Nb partition similarly during mantle melting (Saal et al., 2002;Rosenthal et al., 2015). Given this similarity in partitioning behaviour, it should be expected that pre-degassing magmatic values of CO 2 /Ba and CO 2 /Nb will not be fractionated from their mantle sources; the exception to this being that even small differences between the partition coefficients can cause large variations in the CO 2 /Ba and CO 2 /Nb of instantaneous fractional melts extracted from a previously depleted source (Rosenthal et al., 2015).
A number of previous studies (Saal et al., 2002;Le Voyer et al., 2017;Hauri et al., 2018) are predicated upon correlations between CO 2 and Ba or Nb indicating a lack of any fractionation between CO 2 and the incompatible trace elements, including the absence of any CO 2 degassing prior to, or after, inclusion entrapment. However, Matthews et al. (2017) demonstrated that such correlations are also a natural consequence of degassing and mixing of compositionally diverse near-fractional mantle melts, and such a process can explain the increasing variation in CO 2 /Ba with Ba concentration that many datasets exhibit.
Magma mixing is common in magmatic systems (Sobolev, 1996;Maclennan, 2008a;Shorttle, 2015;Shorttle et al., 2016), and the most CO 2 rich fractional melts are expected to saturate in CO 2 vapour at pressures within, or below, the crust (Dixon et al., 1995). Though the degassing and mixing process causes the average value of CO 2 /Ba within a dataset to be controlled primarily by the pressure of degassing, the maximum values seen in a dataset may still record the mantle value (or are at least a minimum estimate of the mantle value). Datasets that do not exhibit correlations likely have a more complex multi-stage degassing history, making them less reliable records of mantle CO 2 /Ba ratios. Shimizu et al. (2019) recently demonstrated that not all melt inclusion datasets require a partialdegassing and mixing process to explain their CO 2trace element systematics. For inclusions from the Siqueiros and Garrett transform faults, Shimizu et al. (2019) showed that the covariance of CO 2 with a broad suite of trace elements is compatible with the trapped melts having a single mantlederived CO 2 /Ba ratio for variable Ba concentration. The observed variability in CO 2 /Ba and CO 2 /Nb ratios arises from analytical imprecision. However, partially degassed datasets can also give rise to melt inclusions with a single CO 2 /Ba ratio at variable Ba concentration. If the trace-element and CO 2 enriched melts degas, but are then efficiently homogenised, a single enriched mixing endmember is produced, with a CO 2 /Ba ratio lower than the mantle value (Matthews et al., 2017, Appendix A).
When mixed with extremely trace-element and CO 2 depleted melts, a binary mixing array in CO 2 -Ba space will be generated with each melt having an identical CO 2 /Ba ratio. In this scenario, the mean CO 2 /Ba ratio of the datasets would then be only a minimum bound on the mantle value. The implications and likelihood of such a process for the Siqueiros and Garrett datasets are discussed in Section 5.7.
To understand what CO 2 /Ba and CO 2 /Nb ratios tell us about mantle carbon heterogeneity we need to consider how these elements are transferred to the melt. During mantle melting, the residue is likely to be completely depleted of CO 2 after only a few percent melting, either by the formation of carbonatitic or carbonated-silicate melts if a carbon phase is present (Dasgupta and Hirschmann, 2010), or by virtue of the incompatibility of CO 2 in silicate minerals (Saal et al., 2002;Rosenthal et al., 2015). Even refractory mantle lithologies are likely to have entirely lost their carbon to a melt phase long before silicate melting begins. Since such low-degree melts have a vanishingly small chance of escaping mixing during transport, we assume that melts in one location will not preserve the variable CO 2 /Ba and CO 2 /Nb ratios of small-scale heterogeneities, but instead will inherit an average value. Larger scale variations in the relative proportions of small-scale heterogeneities may allow the relative contributions of different mantle components to the CO 2 budget of the magma to be deconvolved.

Identifying mantle CO 2 heterogeneity
In this contribution we seek to assess the presence of carbon rich mantle reservoirs in the Earth. In order to deconvolve the signals of CO 2 degassing and mantle heterogeneity, a large number of melt inclusions from eruptions sampling a diversity of mantle sources need to be analysed for both their trace element and CO 2 concentrations. We approach this by considering mantle carbon heterogeneity on the scale of both a single well sampled mantle plume (Iceland), and a global scale with the more sparsely sampled mid-ocean ridges and ocean islands. Iceland offers an excellent opportunity for studying carbon heterogeneity within a single mantle plume due to the large number of previous melt inclusion studies (Hauri et al., 2018;Hartley et al., 2014;Neave et al., 2014;Bali et al., 2018;Schipper et al., 2016;Miller et al., 2019), which we supplement with four new datasets (Sections 2 and Section 3.1).
Comparisons between the CO 2 -trace element systematics present within the Iceland melt inclusion compilation and the global compilation enable us to make a new assessment of the effects of crustal storage and melt inclusion decrepitation (Section 4). Though this secondary crustal processing signal dominates the global CO 2 /Ba and CO 2 /Nb arrays, in Section 5 we demonstrate that a signal of mantle source CO 2 heterogeneity still remains in the Iceland compilation. By inversion of a three component mixing model (Section 5.4) we estimate CO 2 concentrations, and their uncertainties, in endmember mantle components within the Iceland plume. We find evidence that CO 2 concentrations higher than in the depleted mantle are present within the primordial mantle, and lower in recycled mantle (though poorly constrained).

Samples and Methods
Samples were collected from four fresh primitive eruptions in Iceland ( Figure 1) which represent diverse mantle sources, as indicated by their Sr-, Ndand He-isotope ratios. Háleyjabunga and Stapafell are in close proximity on the Reykjanes Peninsula in the Western Rift Zone, but preferentially sample more depleted and enriched mantle components respectively (Thirlwall et al., 2004). Stapafell erupted sub-glacially between 70 and 14 ka (Saemundsson et al., 2010) forming basal pillow basalts, from which samples were taken near 63 • 54.585'N, 22 • 31.409'W. Háleyjabunga was erupted as a subaerial lava shield at ∼13 ka (Saemundsson et al., 2010), from which olivine-phyric lava flow samples were taken from the Eastern side of the vent near 63 • 48.978'N, 22 • 39.099'W. Berserkjahraun is an eruption younger than 11 ka in the Snaefellsnes flank zone (Hjartarson and Saemundsson, 2014), and has extreme geochemical enrichment (Peate et al., 2010). Glassy olivine-and plagioclase-phyric scoria was collected from the crater at 64 • 95.915'N, 22 • 89.853'W. Heilagsdalsfjall is a lava flow erupted in the Northern Volcanic Zone with an age < 11.5 ka (Saemundsson et al., 2012). Glassy olivinephyric achneliths were collected from the crater at 65.49934 • N, 16.70224 • W.
Olivine and clinopyroxene crystals containing melt inclusions were extracted from each sample. The total CO 2 content of a melt inclusion is partitioned between the glass and vapour bubble (if present). Following Hartley et al. (2014), crystals were sequentially prepared for micro-Raman analysis of CO 2 vapour density in the bubbles, Secondary Ion Mass Spectrometry (SIMS) to determine CO 2 and trace element concentrations in the glass, and electron probe microanalysis (EPMA) to determine the major element chemistry of the inclusions and their crystal hosts. Any CO 2 detected in the inclusion-hosted bubbles was added back into the glass, assuming the inclusion was a homogeneous melt phase at the time of entrapment. Full information on sample preparation and analytical setup is given in the Supplementary Information. In Appendix A we explain why we do not apply a correction based on predicting the CO 2 vapour density in bubbles. As we are primarily interested in trace element ratios, and CO 2 concentrations at the time of eruption, we do not apply post-entrapment crystallisation corrections (though we provide sufficient information in Supplementary Tables that this may be performed).

Results
The four new datasets from Berserkjahraun ,Háleyjabunga, Heilagsdalsfjall, and Stapafell are described in Section . In Section 3.2 we consider our new datasets alongside others from Iceland to assess whether the mantle CO 2 /Ba ratio has survived the melting process. Finally, we discuss the processes most likely to be controlling the CO 2 -trace element systematics (Section 3.3).

CO 2 and trace elements
Melt inclusion CO 2 and Ba concentrations, with and without bubble corrections, are shown for each eruption in Figure 2. All four eruptions show considerable variation in their CO 2 /Ba ratios, but only Háleyjabunga and Stapafell display positive correlations between CO 2 and Ba. Only three Háleyjabunga melt inclusions analysed by SIMS contained bubbles in which CO 2 vapour was detected (a further three were measured by Raman only, as shown in Supplementary Figure 1d). Though five Stapafell melt inclusion hosted bubbles contained CO 2 vapour (Supplementary Figure 1e), only one was in an inclusion analysed by SIMS. No CO 2 vapour was detected in any of the bubbles hosted in Heilagsdalsfjall inclusions. A substantial fraction of Berserkjahraun inclusions contained bubbles in which CO 2 vapour was detected. In Háleyjabunga, Stapafell and Berserkjahraun, the highest CO 2 concentrations are in inclusions where corrections have been applied. For Háleyjabunga and Stapafell, the corrected inclusions do not have higher CO 2 /Ba than uncorrected inclusions. However, the corrections applied to the Berserkjahraun inclusions are extremely large and consequently significantly increase the CO 2 /Ba ratios recorded by the melt inclusion population. The large magnitude of the corrections arises from high CO 2 vapour densities up to 0.75 g cm −3 . Though this density is above the maximum density of vapour in coexistence with liquid CO 2 at room temperature (∼0.13 g cm −3 ), any heating by the Raman instrument's laser could move the bubble beyond the CO 2 triple point. The bubble volume fraction exerts some control on the corrected CO 2 concentration, but little dependence on CO 2 vapour density is observed (Supplementary Figure 2). A sub-population of three Háleyjabunga melt inclusions have much higher (bubbleuncorrected) CO 2 concentrations than the general population (distinguished on Figures 2 and 3), the origin of this feature is discussed further in Section 3.3.
Trace element concentrations in the four melt inclusion suites are shown in Figure 3. Melt inclusions from Háleyjabunga ( Figure 3a) show ex-treme variability in relative trace element enrichment, Maclennan (2008b) and Neave et al. (2018) argued this variability is most likely mantle derived. The Háleyjabunga matrix has a depleted trace element pattern (Skovgaard et al., 2001). Four Háleyjabunga inclusions have anomalously high Ba and Nb concentrations relative to the light rare earth elements; since their CO 2 concentrations are similar to the main population of inclusions they have not retained high CO 2 /Ba or CO 2 /Nb ratios. The Háleyjabunga inclusion with highest CO 2 /Ba ratio is part of the high-CO 2 population, but as we discuss in Section 3.3 this may not reflect a mantle value, and so we treat this inclusion as an anomaly here. Ignoring this anomalous inclusion, the highest CO 2 /Ba ratios (>200) are observed in inclusions with depleted trace element patterns, but are not anomalously depleted among the larger population of Háleyjabunga inclusions.
The Stapafell inclusions with highest CO 2 /Ba ( Figure 3b

Fractionation of Ba, Nb and CO 2 during melting
Primary magmas will only preserve mantle CO 2 /Ba and CO 2 /Nb ratios if Ba, Nb and CO 2 are not fractionated from each other during mantle melting. Crustal processes are not expected to fractionate Ba from Nb, but magma degassing will cause CO 2 to be fractionated from both. The melt inclusions least affected by degassing can be identified by taking the highest CO 2 /Ba ratios (Matthews et al., 2017).
However, neither this filtering process, nor observations of the distribution of CO 2 -trace element ratios, can conclusively rule out small degrees of fractionation of CO 2 from Ba and Nb during melting. Experimental work by Rosenthal et al. (2015) suggests that during mantle melting CO 2 should behave marginally more compatibly than Ba, and less compatibly than Nb. Therefore, if Nb and Ba have not been fractionated from each other during melting, then it is unlikely that CO 2 was fractionated from either. Fractionation of Ba and Nb will generate variations in the Ba/Nb ratio of melts with Ba and Nb concentration. Figure 4a shows this is not seen in any of the eruptions in this study, or the other Icelandic eruptions in the compilation. The scatter in Ba/Nb ratio arises largely from analytical imprecision, but could also reflect small scale source heterogeneity.
Another approach to assessing elemental fractionation during melting is to compare the variability of trace element concentrations. The more incompatibly an element behaves during melting, the greater its variability amongst instantaneous fractional melts of the mantle. If two elements in a dataset have the same mean-normalised variance, the melting process has not fractionated the elements from each other (Rudge et al., 2013). The relative variance of Ba and Nb for each Icelandic melt inclusion suite are shown together in Figure  4b, providing further evidence that Ba and Nb have not been fractionated from each other. Elements more compatible than Nb have been fractionated from each other, as demonstrated by their relative variances (Supplementary Figure 3).
Assuming that the behaviour of CO 2 during mantle melting can be modelled as an incompatible element, and that it has a partition coefficient between that of Ba and Nb (Rosenthal et al., 2015), the apparent lack of fractionation between Ba and Nb would suggest that primary magmatic CO 2 /Ba and CO 2 /Nb ratios will record mantle source ratios.
3.3. Processes controlling CO 2 -trace element systematics Háleyjabunga and Stapafell both display correlations between CO 2 and Ba concentrations (Figure 2). Additionally, Háleyjabunga shows the increased variance of CO 2 concentration with Ba concentration that is typical of datasets generated by partial degassing and mixing (Matthews et al., 2017). Stapafell exhibits somewhat more complexity, possibly due to a more prolonged history of mixing and fractional crystallisation. Heilagsdalsfjall and Berserkjahraun do not exhibit correlations between CO 2 and Ba concentrations, likely a result of more extensive degassing and mixing.
Three Háleyjabunga melt inclusions (without bubble corrections) form a sub-population with higher CO 2 concentrations than the other Háleyjabunga melt inclusions ( Figure 2). The higher CO 2 concentrations may reflect inclusion entrapment at an earlier stage of magma evolution, where the melts had higher CO 2 concentrations due to a greater storage pressure or by less mixing having taken place. The highest CO 2 /Ba ratio in the Háleyjabunga melt inclusions also comes from within this sub-population; if the population is derived simply from being trapped earlier, then this is the value of CO 2 /Ba most likely to reflect the Háleyjabunga mantle source. However, to fully assess this hypothesis by consideration of the sub-population's CO 2 -trace element systematics requires a much larger sample size. Alternatively, these inclusions may reflect CO 2 being added to depleted melts,   (Rudge et al., 2013), using the partition coefficients from Workman and Hart (2005).
possibly by dissolution of CO 2 vapour or by melting in the presence of graphite. Both processes can elevate CO 2 concentrations and CO 2 /Ba ratios of depleted melts, and are discussed further in Appendix B. Since we cannot be certain these inclusions represent mantle CO 2 /Ba ratios, we choose to discount them from the following discussion.

Systematics of the Global Melt Inclusion Array
In Figure 5 we compare our new CO 2 and Ba data from Háleyjabunga, Stapafell, Heilagsdalsfjall and Berserkjahraun, with other Icelandic melt inclusion suites and those from ocean-island and midocean ridge settings. We exclude inclusions from arc volcanoes since their enrichment in H 2 O complicates our understanding of CO 2 solubility. In both the Iceland compilation ( Figure 5a) and the global compilation (Figure 5b), there is a striking negative correlation between CO 2 /Ba and Ba concentration. This correlation is bounded above and below by lines of constant CO 2 concentration, though bubble-corrected inclusions (unfilled symbols) break through this upper bound. For the most depleted inclusions the low bound corresponds to approximately 100 ppmw CO 2 , and the high bound to 1000-1400 ppmw CO 2 . The same systematic is observed in CO 2 /Nb-Nb space (Supplementary Figure 4).

4.1.
A global array reflecting mantle CO 2 heterogeneity? The negative correlation between CO 2 /Ba and Ba concentration could be caused by source heterogeneity. The trace-element enriched mantle reservoirs that contribute disproportionately to the Baenriched melts may not have corresponding CO 2 enrichments. However, as the negative correlation has the same position and gradient in both the Icelandic and global datasets, this would imply depleted and enriched mantle components occur with the same distinct CO 2 /Ba and CO 2 /Nb ratios everywhere. Radiogenic isotope constraints demonstrate that there is substantial variation in the chemistry of enriched and depleted components globally, even between highly incompatible elements that are difficult to fractionate during partial melting (e.g. U and Rb) (Stracke et al., 2005). It is therefore perhaps unlikely that these heterogeneous components would be characterised by single CO 2 /Ba and CO 2 /Nb ratios. Furthermore, the presence of CO 2 -rich shrinkage bubbles in many of these inclusions demonstrates that many of the melts must have contained more CO 2 in the past, which would have placed them originally above the global array. This is particularly apparent for melt inclusions from Berserkjahraun, Miðfell, Laki and Surtsey (Figure 5a).

A global array controlled by degassing?
A simpler explanation than source heterogeneity is the operation of a process limiting the CO 2 concentration within the melt inclusion. The most obvious process is degassing. In this case, the lower limit of CO 2 on the array corresponds to CO 2 solubility for magma storage in shallow crustal magma chambers. The small number of inclusions with lower CO 2 concentrations than this probably owe their undersaturation to degassing during eruption, or mixing with extremely CO 2 -depleted partial mantle melts (Matthews et al., 2017). The upper bounds on CO 2 concentration in depleted melt inclusions of 1000-1400 ppmw correspond to saturation pressures in the range of 1-3 kbar, depending on the solubility model chosen (Figure 6, Supplementary Figures 5 and 6). This pressure range coincides with that expected for olivine decrepitation (Wanamaker et al., 1990;Maclennan, 2017), where the crystal undergoes brittle failure when it can no longer support the pressure difference between the inclusion and its surroundings. Experiments by Wanamaker et al. (1990) predict that overpressures of up to 2.2 kbar may be supported. This is very close to the maximum entrapment pressure of the moderately-depleted inclusions when calculated with the CO 2 solubility model of Iacono-Marziano et al. (2012), demonstrated in Figure 6. Figure 5b demonstrates the presence of two different upper limits; one for more depleted eruptions and one for more enriched eruptions. If decrepitation is the process responsible for the upper limit of the array, the limits for both enriched and depleted eruptions ought to correspond to the CO 2 solubility at 2.2 kbar. However, the dependence of CO 2 solubility on melt composition allows significantly different CO 2 concentrations to equilibrate at the same pressure in different melts. Melt polymerisation tends to reduce CO 2 solubility, whilst the presence of cations with an affinity for forming carbonate complexes increases CO 2 solubility (Shishkina et al., 2014).
The enriched inclusions are significantly richer in total-alkalis (Na 2 O + K 2 O) for similar SiO 2 con-tents (Supplementary Figure 7), but this does not translate into significant differences in the π * compositional parameter used by Shishkina et al. (2014) in their CO 2 solubility model ( Supplementary Figure 8). Though the concentration of alkali elements in the enriched inclusions acts to increase the CO 2 solubility, their lower MgO and CaO concentrations tends to decrease CO 2 solubility (Dixon, 1997). As the decrepitation limit should not vary strongly with eruption enrichment, it is likely that the higher CO 2 bound for the enriched melts is solubility related, even if not captured by π * .
A compositional-solubility effect creating the different bounds to the global array is supported by recent experimental work by Allison et al. (2019). By performing new experiments at intermediate pressure and alkaline compositions, Allison et al. (2019) show that existing CO 2 solubility models provide inadequate predictions of CO 2 solubility in more alkaline melts. Though the least alkaline experimental melt composition used by Allison et al. (2019) is more alkaline than the majority of eruptions from which melt inclusions are displayed in Figure 5, their empirical expression for CO 2 solubility likely incorporates other compositional effects of importance for the more enriched eruptions considered here. Figure 6 demonstrates that applying the empirical expression to the more enriched eruptions brings their saturation pressures into line with the depleted eruptions, when the CO 2 solubility model by Iacono-Marziano et al. (2012) is applied to the inclusions from depleted eruptions. However, a discrepancy exists if the saturation pressures for depleted eruptions are calculated using the models by Shishkina et al. (2014) and Eguchi and Dasgupta (2018a) (Supplementary Figure 5).

Identifying mantle CO 2 heterogeneity
The gross structure of both the Icelandic and Global melt inclusion arrays shown in Figure 5 is primarily controlled by low-pressure processes: storage in shallow magma chambers and olivine decrepitation (Section 4). However, many eruptions sit almost entirely within the bounds imposed by shallow storage and decrepitation, or have partitioned a significant fraction of their CO 2 budget into bubbles and thereby overcoming the decrepitation threshold. Many of these datasets may still preserve the CO 2 /Ba ratios of their melts at the time of entrapment.

This Study
Miðfell (Miller et al., 2019)   However, all the melt inclusion datasets considered here show significant internal variability in the CO 2 /Ba ratios they preserve. Matthews et al. (2017) argued that the diversity of CO 2 /Ba ratios within eruptions indicates that they have experienced partial CO 2 loss by degassing, but can regain a positive (albeit scattered) correlation between CO 2 and Ba by magma mixing. We consider the implications of concurrent degassing and mixing for identifying mantle CO 2 /Ba in Section 5.1. Another plausible explanation for diversity in CO 2 /Ba ratios is the mapping of source CO 2 /Ba heterogeneity into its melts. In Section 5.3 we argue such small heterogeneities are likely to be homogenised for individual eruptions during melt generation and transport. Applying the methodology developed in Sections 5.1 and 5.3 to the melt inclusion compilation from Iceland, supplemented with our new data, we identify and place quantitative constraints on CO 2 heterogeneity within the Icelandic mantle plume (Section 5.4). We then explore in Section 5.7 whether this approach yields information on the global variability of CO 2 in the upper mantle, as sampled by mid-ocean ridge basalts (MORBs).

Has concurrent degassing and mixing removed
the mantle CO 2 /Ba signal? A consequence of concurrent degassing and magma mixing is that the average value of CO 2 /Ba for an eruption is more sensitive to the degassing pressure than to the CO 2 /Ba ratio inherited by the primary melts from their mantle source. Since the degassing, mixing and decrepitation processes all act to lower the CO 2 /Ba ratio recorded by melt inclusions, it follows that the maximum CO 2 /Ba value observed in melt inclusions from a single eruption is the value least affected by CO 2 loss and therefore most likely to preserve the mantle CO 2 /Ba ratio.
Any given melt inclusion dataset samples a population of melts with CO 2 /Ba ratios ranging from very low up towards the mantle value. The likelihood of the dataset to sample the true maximum CO 2 /Ba ratio (the mantle value) depends on the relative proportions of these melts, described statistically by a probability distribution. An eruption's CO 2 /Ba distribution may have a long lowprobability tail towards the mantle value, depending on the extent of degassing prior to mixing and the efficiency of mixing. Without knowledge of the distribution of melts entering the crust, and the history of degassing and mixing they subsequently un-dergo, it is impossible to constrain the shape of this distribution and therefore assess the likelihood of a given melt inclusion suite preserving the mantle CO 2 /Ba ratio in one or more of the inclusions.
Though we cannot tell whether an individual dataset preserves the mantle CO 2 /Ba ratio when taken in isolation, comparison with a large number of melt inclusions from a diversity of eruptions can provide a context in which to consider the maximum CO 2 /Ba values preserved by an individual eruption. If the maximum CO 2 /Ba ratios of eruptions vary systematically with their mantle sources, and systematic biases from crustal processing can be ruled out, then greater confidence may be placed in the recorded CO 2 /Ba ratios reflecting their mantle source.
To avoid interpreting analytical artefacts, Rosenthal et al. (2015) and Hirschmann (2018) avoid using the highest CO 2 /Ba values, suggesting these analyses may have incorporated carbon derived from cracks in the sample. Standard practice when collecting melt inclusion datasets is to avoid analyses in the vicinity of cracks, and the 12 C count rates would vary erratically in such a situation. In carefully collected datasets such analytical artefacts will be eliminated, and it is more likely that the highest CO 2 /Ba ratios instead reflect sampling of low probability-density tail to the CO 2 /Ba distribution. Furthermore, the inclusions with the highest CO 2 /Ba and CO 2 /Nb have CO 2 concentrations that are within the range of the rest of the dataset, and not anomalously high as might be expected if it were an analytical artefact.
In Supplementary Figures 9 and 10 we demonstrate the difference between taking the dataset maximum CO 2 /Ba and CO 2 /Nb ratios, as opposed to the upper limit of the main CO 2 /Ba and CO 2 /Nb data density. The dataset maximum CO 2 /Ba ratios do not coincide with extreme values of Ba/Nb (Supplementary Figure 11) or anomalous trace element chemistry ( Supplementary Figure 12), suggesting the high values are also not derived from analytical artefacts in trace element concentration measurements. The lack of coincidence between high CO 2 /Ba and extreme Ba/Nb provides further evidence that the highest CO 2 /Ba ratios do not arise from fractionation between Ba, Nb and CO 2 during melting.

Preservation of small-scale mantle CO 2 heterogeneity in melts
In applying the method described above, we are implicitly assuming that the CO 2 /Ba ratio of a single melt inclusion is representative of the heterogeneous mantle that contributes melt to the whole eruption.
However, mantle-derived Pbisotope heterogeneity is present within melt inclusions from single Icelandic eruptions, including single hand-specimens from Háleyjabunga and Stapafell (Maclennan, 2008b). This observation clearly demonstrates that heterogeneity in Pbisotope ratios survives the melting and melt transport processes operating beneath Iceland.
However, CO 2 , Ba and Nb are much more incompatible than Pb and are likely to be almost entirely removed from the solid residue in the first increments of melting. Though channelised mantle flow can transport significant mantle derived chemical heterogeneity in melts, the deepest melts are likely to be well mixed (Spiegelman and Kelemen, 2003;Kelemen et al., 1997), and may subsequently be diluted in varying amounts by higher degree melts. We therefore consider it likely that single eruptions will have mixed out any heterogeneity in CO 2 /Ba and CO 2 /Nb ratios that was present within their mantle source regions. If this is correct, we must look for carbon heterogeneity on an eruption to eruption basis rather than within single eruptions.
By taking the maximum CO 2 /Ba ratio of a melt inclusion dataset we are assuming this is representative of the average mantle sampled by the eruption, as argued above. When we make the comparison between this CO 2 /Ba value, and the bulk-rock radiogenic isotope ratios, we must assume that the isotope ratios are representative of the same mantle average. Since Pb, Sr, and Nd are more compatible than CO 2 and Ba, melts containing these elements are generated at higher melt fractions, and source heterogeneity can survive melting and transport (Spiegelman and Kelemen, 2003;Kelemen et al., 1997). If the whole rock lava from which the radiogenic isotope measurements are made does not reflect complete mixing of these diverse melts, the Pb-, Sr-, and Nd-isotope ratios may be decoupled from the maximum CO 2 /Ba and CO 2 /Nb seen in melt inclusions.
Additionally, higher extents of melting are required for a mantle component to contribute to the Nd, Sr, and Pb budget of a magma, than for it to contribute to the CO 2 , Ba, and Nb budget. Where the extent of melting is controlling the relative contribution of depleted and enriched mantle to the Pb-, Nd-, and Sr-isotope ratios, rather than source composition, the radiogenic isotopes and the maximum CO 2 /Ba and CO 2 /Nb in melt inclusions may reflect averages of different proportions of the enriched and depleted components in the source.
For the influence of these processes to be properly assessed a much larger number of plausibly undegassed melt inclusion datasets is required. We, therefore, do not consider the possibility of CO 2 /Ba and radiogenic isotope decoupling further. A more robust comparison could be made with radiogenic isotopes of an element at least as incompatible as Nb. One such system is 3 He/ 4 He which we use in complement to Pb-, Nd-and Sr-isotopes in the following analysis.

Heterogeneity within the Icelandic plume
Radiogenic isotope observations require multiple enriched and depleted source components within the Icelandic mantle (Thirlwall et al., 2004;Peate et al., 2010;Shorttle et al., 2013). We will use our new data to evaluate whether these mantle components also differ in their CO 2 content. The Pb-, Nd-, and Sr-isotope systems are sensitive to ancient melt extraction and recycling events, whilst the 3 He/ 4 He system is diagnostic of primordial (or undegassed) mantle reservoirs.
Though there is a large quantity of melt inclusion data from Iceland, from a number of spatially and temporally restricted eruptions with diverse geochemistry, it is not yet sufficient to fully resolve the CO 2 concentrations in the diversity of mantle components beneath Iceland. We therefore proceed on the basis of finding the minimum number of distinct components required to explain the coupled CO 2 /Ba and radiogenic isotope observations, whilst acknowledging the end-members we invert for are some average of the full diversity of heterogeneities present in the Icelandic mantle. Figures 7a and 7b demonstrate there is no systematic covariation of CO 2 /Ba within 206 Pb/ 204 Pb-208 Pb/ 204 Pb space, nor 87 Sr/ 86 Sr-143 Nd/ 144 Nd space. Though not shown here, this is also true of any combination of Pb-isotopes with Sr-and Nd-isotopes. Whilst the lithophile isotope systems excel at separating contributions from enriched and depleted mantle reservoirs, the primordial mantle lies in the middle of lithophile isotope space (Hart et al., 1992). Helium isotopes, in contrast, reveal contributions from primordial man-tle components clearly, as primitive mantle has extremely high 3 He/ 4 He ratios (e.g., Class and Goldstein, 2005).
When 3 He/ 4 He ratios are plotted with 143 Nd/ 144 Nd ratios (Figure 7c) a covariation of CO 2 /Ba with position across the space does emerge: mantle material with high (primordial) 3 He/ 4 He is associated with high magmatic values of CO 2 /Ba, whilst melts derived from depleted mantle (with R/R a ∼ 8) have moderate CO 2 /Ba and those from enriched mantle (with R/R a ≤ 8) have low CO 2 /Ba. The same trends are observed if CO 2 /Nb is used in place of CO 2 /Ba (Supplementary Figure 13), and a similar co-variation is observed if the upper-limits of the CO 2 /Ba and CO 2 /Nb ratios in the main data population are used in place of the dataset maxima ( Supplementary Figures 14 and 15).
The association of high CO 2 /Ba with primordial undegassed mantle material was first made by Miller et al. (2019), based on the elevated CO 2 /Ba ratios seen in the Miðfell eruption in Iceland. By making a comparison to melt inclusions from Borgarhraun (Hauri et al., 2018), Miller et al. (2019) argue that the elevation in CO 2 /Ba in Miðfell is not caused by crustal processing, since both eruptions have similar trace element concentrations. Though it could be argued that Borgarhraun has seen more extensive CO 2 loss by magmatic degassing prior to inclusion entrapment, this seems unlikely given the high (6-8 kbar) pressures of crystallisation (Winpenny and Neave and Putirka, 2017) and the Borgarhraun melt inclusions' distance from the decrepitation threshold ( Figure 5a). Our new data from Háleyjabunga lends further support to hypothesis of Miller et al. (2019) that primordial mantle domains are associated with high CO 2 /Ba. However, the same argument cannot be used when making a comparison to the more enriched eruptions, such as Berserkjahraun and Stapafell, as they have significantly higher trace element concentrations in their melt inclusions. If these melts started off with the same or higher CO 2 /Ba ratio than Borgarhraun, the Stapafell and Berserkjahraun melts should have considerably higher CO 2 concentrations. Such high CO 2 concentrations make these melts more susceptible to degassing and the inclusions more susceptible to CO 2 gas loss during decrepitation; crustal processing will therefore tend to introduce a systematic bias towards low CO 2 /Ba ratios in enriched magmas. Regardless of  Supplementary Table 3. The symbols used are the same as in Figure 5, but each symbol is colored for the maximum CO 2 /Ba ratio observed in melt inclusions from that eruption, as described in the text. Grey circles show the compilation of Iceland whole-rock isotope analyses, data sources are given in the Supplementary Information.
whether the low CO 2 /Ba in the enriched Berserkjahraun magmas is a result of low CO 2 /Ba in the enriched source or not, we require at least three end-members to deconvolve the primordial mantle signal.
5.4.1. Inverting for mantle endmember CO 2 concentrations Here we consider mixing between three mantle end-members: 1. Depleted mantle (DM)-containing low trace element concentrations with radiogenic Nd-and He-isotope ratios (R/R a ∼ 8), typical of the MORB-source mantle 2. Enriched mantle (EM)-containing higher trace element concentrations with unradiogenic Nd-isotopes and radiogenic He-isotope ratios (R/R a ≤ 8), representing recycled material 3. Primordial mantle (PM)-containing trace element concentrations and isotope ratios typical of bulk silicate Earth.
The numerical values of their geochemical properties are given in Table 1.
In principle, a mixing surface which closely matches the covariation of CO 2 /Ba across 143 Nd/ 144 Nd-3 He/ 4 He space can be found. Sohn (2013) developed a mathematical formulation for n-dimensional mixing hyperbolae, when mixing between three end-members this reduces to a surface in three-dimensions. Sohn (2013)  (2) The subscripts 'DM,' 'PM,' and 'EM' represent endmember values for the depleted, primordial and enriched components respectively. The parameters K 1 , K 2 , K 3 and K 4 represent ratios of source trace element concentration ratios: i.e., source trace element concentrations will trade off against each other to yield identical surfaces. By combining the concentration ratios into four parameters, the number of free parameters in the inversion is reduced (Sohn, 2013).
In order for our prior assumptions about mantle source chemistry to be applied consistently, we take a Bayesian approach to estimating the CO 2 /Ba ratio and CO 2 concentration of the mantle sources. The prior distributions for end-member properties are described in Table 1 and its caption. For the K n parameters log-uniform distributions from 10 −6 to 10 6 were used. We implement the Bayesian calculation using the importance nested sampling Monte-Carlo inversion routine 'Multinest' (Feroz and Hobson, 2008;Feroz et al., 2009Feroz et al., , 2013 with the pymultinest python wrapper (Buchner et al., 2014). For simplicity we ignore the uncertainty in the 143 Nd/ 144 Nd and 3 He/ 4 He ratios of the samples, calculating the log-likelihood of each mixing model based only on the misfit between the CO 2 /Ba of the eruption and modelled mixing surface. The total log-likelihood is calculated by summing contributions from each eruption: ln(L) = n 1 ln (L n (µ n , σ n , x n )) where µ n is the observed CO 2 /Ba ratio, and σ n is its one standard deviation uncertainty. The value of CO 2 /Ba ratio on the modelled mixing surface (x n ) is calculated by numerically solving Equation 1, given the observed 143 Nd/ 144 Nd and 3 He/ 4 He ratios of the samples. Since the Beserkjahraun, Stapafell, and Surtsey melt inclusions have maximum CO 2 /Ba ratios that fall close to the decrepitation threshold ( Figure 5), we treat these observations as minimum bounds on the CO 2 /Ba mixing surface. Models which place the mixing surface at lower values of CO 2 /Ba ratio are penalised by an exponential contribution to the log-likelihood: Haleyjabunga, Miðfell and Borgarhraun preserve inclusions with high CO 2 /Ba ratios considerably below the decrepitation threshold; therefore, they provide the most direct constraint on the mixing surface. Their contribution to the log-likelihood function reflects their squared normalised distance from the surface: This expression is equivalent to the expression for log-likelihood assuming the constraints follow a normal distribution, with the constant terms removed (Lee, 2012). Inverting the Sohn (2013) mixing equations provides estimates for only CO 2 /Ba ratios of the three end-members, and not their CO 2 concentrations. To estimate CO 2 concentration we employ a Monte-Carlo calculation pairing values drawn from the posterior CO 2 /Ba distributions of each endmember with a randomly selected value from a prior distribution of Ba concentrations of each mantle source (Table 1).

Inversion results
The results of the inversion are shown in Figure 8 and Table 2, they represent the first CO 2 concentration estimates to be made for different mantle components in one location. Figure 9 shows that our estimate for the depleted mantle CO 2 concentration (102 +21 −14 ppmw) compares favourably with the concentrations inferred from MORB melt inclusion suites (Rosenthal et al., 2015;Le Voyer et al., 2017), undersaturated glasses (Michael and Graham, 2015) and with estimates derived from CO 2 -He systematics (Marty, 2012;Tucker et al., 2018). Our new estimate for primordial mantle (2.15 +0.94 −1.49 wt%) is much higher than previous estimates of the CO 2 content of plume mantle, or bulk mantle as estimated by Marty (2012), whilst our estimate of CO 2 concentration in the recycled mantle (68 +1170 −64 ppmw) is significantly lower than the plume and bulk mantle estimates. This difference most likely arises from previous estimates averaging high-and low-CO 2 mantle reservoirs.

The Depleted Mantle
When the DM CO 2 /Ba prior is left open, the inversion routine finds two likelihood maxima; one solution places the high mantle CO 2 concentrations into the DM-endmember rather than the PMendmember. Since previous work has shown the depleted mantle has moderately low CO 2 concentrations (e.g., Saal et al., 2002;Hirschmann, 2018; Figure 8: Results of the Bayesian inversion for mantle endmember CO 2 concentrations, assuming Haleyjabunga, Miðfell, and Borgarhraun (large symbols in panel a) preserve mantle CO 2 /Ba ratios, whilst Berserkjahraun, Stapafell, and Surtsey (small symbols in panel a) provide only minimum bounds. The colour shading in panels a and b represents the CO 2 /Ba ratio of the fitted surface, with the scale saturating below 50 and above 1200. The eruptions are coloured for their CO 2 /Ba observations on the same scale. The black contours on panel a also show CO 2 /Ba ratio. The dashed red lines on panels a and b show the extreme mixing boundaries. Panels c-e show the posterior distributions for end-member CO 2 /Ba. Panel f shows the trade-off between PM CO 2 /Ba and K 2 as a 2D histogram with shading indicating data density. Panels g-j show the posterior distributions for K 1 -K 4 . The vertical lines in panels g and i show the values for the solid mantle calculated using Nd and Ba estimates from Workman and Hart (2005) and Stracke et al. (2003). Panels k-m show the posterior CO 2 concentrations calculated from the posterior CO 2 /Ba distributions. Panel n shows the trade-off between mantle CO 2 concentration and K 2 .  Kurz and Jenkins (1981) 5.0-7.4

Parameter
See caption d

50-95
See caption e Ba (ppmw) 0.58±0.32 Workman and Hart (2005) 2.55 ± 0.25 Stracke et al. (2003) f 6.9 ± 1.0 Palme and O'Neill (2003) Table 1: The priors set on the parameters in the inversions reported in the main text. The sources from which the priors are based are shown in the table. Where ranges are quoted for isotope ratios a uniform (or log-uniform) distribution between those bounds is used, for trace element concentrations it is a log-uniform distribution. Otherwise a normal distribution with the quoted mean and standard deviation are used. a the prior was chosen in order to obtain the solution where DM has moderate CO 2 concentrations, as inferred from previous studies. b The lower bound is lower than the most extreme basalts (e.g. Stracke et al., 2005), and the upper bound is just lower than the lowest Icelandic basalt. c The low Nd-isotope ratio bound for enriched mantle end-member is taken to be lower than the least radiogenic ocean island basalt (Stracke et al., 2005), and the high bound is taken to be lower than the least radiogenic observation from Iceland (Figure 1). d Lower bound is the R/Ra inferred for the EM2 mantle component by Jackson et al. (2007), and the upper bound is the ingrown R/Ra of the enriched component after 1 Ga. e Lower bound is the highest R/Ra observed in basalts derived from the Iceland plume (Stuart et al., 2003). The upper bound is the R/Ra predicted for entirely undegassed BSE-like mantle (Class and Goldstein, 2005). f Following Shorttle and Maclennan (2011), the trace element concentrations are taken as a third recycled crust (Stracke et al., 2003), and two thirds depleted mantle (Workman and Hart, 2005 Table 2: Summary of the inversion results described in the main text. Median values and the 5% and 95% confidence intervals are given for the CO 2 /Ba, CO 2 /Nb and CO 2 concentration of each of the three mantle endmembers (DM, EM and PM). Hauri et al., 2018), and the primordial component has high CO 2 concentrations (Miller et al., 2019), we choose the DM CO 2 /Ba prior such that only the likelihood-maxima corresponding to this low DM CO 2 solution is found (Table 1).

The Enriched Mantle
A consequence of treating the most enriched eruptions as minimum bounds in the inversion is the EM endmember becomes poorly constrained, and may have a CO 2 /Ba ratio both lower or far in excess of the depleted mantle. Though the median EM CO 2 /Ba ratio (∼10) is significantly lower than for the DM endmember (∼130), the posterior EM CO 2 /Ba distribution depends strongly on the lower bound on the CO 2 /Ba prior. The posterior distribution of K 3 (Figure 8i) is shifted to higher values than that expected from the Nd/Ba ratios of the DM and EM endmembers, shown by the verti-cal lines (Workman and Hart, 2005;Stracke et al., 2003). This offset in K 3 most likely reflects fractionation during melting. The calculated CO 2 /Ba ratio of the DM end-member is comparable to previous estimates for the upper mantle (Rosenthal et al., 2015;Hirschmann, 2018).

The Primordial Mantle
The calculated CO 2 /Ba ratio of the PM endmember is extremely high. However, CO 2 /Ba strongly trades-off against K 2 (the ratio of He/Ba ratios of PM and DM), this reflects the control of the K parameters on the shape of the surface and the distance of the PM end-member from the constraints. No significant trade-off is seen between PM CO 2 /Ba ratio and the PM 3 He/ 4 He and 143 Nd/ 144 Nd ratios within the range of their priors.
Two alternative inversions were also performed. One where we do not assume Berserkjahraun, Stapafell and Surtsey are degassed and instead assume they also preserve mantle CO 2 /Ba ratios (Supplementary Figure 16). In the second inversion we take the most extreme values of CO 2 /Ba ratio observed in Haleyjabunga (Section 3.3) as the mantle source value (Supplementary Figure 17). The resulting surfaces are contorted and do not replicate the observations satisfactorily.

Difficulties in inferring the recycled mantle CO 2 budget
We have demonstrated that reasonable constraints can be placed on the CO 2 budget of depleted and primordial mantle reservoirs. Difficulties remain, however, for identifying the CO 2 content of recycled mantle material. The first of these difficulties is made clear by the quantitative analysis above: eruptions which sample predominantly recycled mantle components are also extremely trace element enriched. If the CO 2 /Ba of the recycled mantle material is similar or greater than the depleted mantle, the high trace element concentrations in enriched basalts results in a greater likelihood of CO 2 loss due to degassing and decrepitation.
A further difficulty is the possibility for decoupling between the extremely incompatible elements CO 2 , Ba and Nb, and the moderately incompatible elements Pb, Nd and Sr. As discussed in Section 5.3, this is unlikely to be an issue for the signal of high CO 2 /Ba with high 3 He/ 4 He, as He also behaves extremely incompatibly during mantle melting (Heber et al., 2007).

Implications of carbon rich mantle reservoirs
The extremely high CO 2 /Ba ratio inferred for the primordial mantle in our inversion (Section 5.4) stems from the observations of high CO 2 /Ba in melt inclusions from two eruptions. One way of rationalising this observation is if the primordial material is extremely depleted in lithophile trace elements, including Ba. Whilst it has been suggested that the high 3 He/ 4 He mantle may be more depleted than bulk silicate Earth (Hart et al., 1992), if the high 3 He/ 4 He mantle were to have similar CO 2 concentrations to the depleted mantle it must be at least an order of magnitude more depleted than the depleted mantle to explain the high observed CO 2 /Ba ratio. Though our estimate for the primordial mantle CO 2 concentration (2.15 +0.94 −1.49 wt%) is higher than most previous mantle CO 2 estimates, it is considerably lower than the carbon content of CI and CM chondrites and similar to the carbon content of the enstatite chondrites, the likely building blocks of the Earth (Wasson and Kallemeyn, 1988).
Our estimate for the primordial mantle CO 2 concentration is primarily constrained by the inclusions from Háleyjabunga and Miðfell. These inclusions have slightly lower CO 2 concentrations (∼500 ppmw), than inclusions from Borgarhraun (∼1000 ppmw), the eruption constraining the depleted mantle CO 2 concentration. Magmas produced from a mantle source as rich in carbon as inferred for the primordial component might be expected to also be extremely CO 2 rich. However, in our model the Háleyjabunga and Miðfell inclusions have only a small contribution from the primordial component ( Figure 8) causing their CO 2 concentrations to become substantially diluted, and allowing them to retain their mantle CO 2 /Ba signal. The low trace element concentrations in Háleyjabunga and Miðfell are further evidence of substantial dilution by extremely trace element depleted melts.
This extremely high primordial mantle CO 2 budget is considerably greater than can be stored in the peridotite mineral assemblage (Rosenthal et al., 2015), though storage in silicate minerals is implicitly assumed when we treat CO 2 as a trace element during melting. For typical mantle oxygen fugacities (f O 2 ) carbon is stabilised in a reduced phase (either graphite or diamond) below the melting region (Dasgupta and Hirschmann, 2010). As mantle decompresses, mineral equilibria shift and cause the reduced carbon to be oxidised to carbonate (CO 2− 3 ). At typical mantle potential temperatures 19 this carbonate-bearing assemblage is above the carbonated peridotite solidus, meaning the transition is accompanied by 'redox' melting (Dasgupta and Hirschmann, 2010). Dilution of this melt by higher degree silicate melts is unavoidable; the carbon liberated from the residue by redox melting will have the apparent behaviour of an extremely incompatible element and therefore remain coupled to Ba and Nb.
If carbon abundances are sufficiently high to influence redox equilibria, or the mantle has anomalously low f O 2 , reduced carbon may persist into the melting region. In Appendix B we demonstrate that the presence of reduced carbon during melting will result in melts with extremely high CO 2 /Ba ratios, and distinctive CO 2 -trace element systematics.

Global upper mantle CO 2 heterogeneity
Despite the unparalleled quantity of melt inclusion data available for Iceland, we are unable to place robust constraints on the carbon content of recycled mantle components. This limitation is due, in part, to having few eruptions constraining the mixing surface close to the EM-component. An alternative approach for assessing the effect of recycled material on the mantle CO 2 budget is to assess global variations in CO 2 /Ba and CO 2 /Nb ratios along the mid-ocean ridge system. The fairly uniform R/R a ∼ 8 for MORBs means the influence of primordial (high 3 He/ 4 He mantle) is excluded from the dataset. We also include eruptions from Iceland that have R/R a ∼ 8 (Borgarhraun and Berserkjahraun), in addition to Heilagsdalsfjall as nearby eruptions all have R/R a ∼ 8 (Harðardóttir et al., 2018), though He-isotope measurements have not been made on Heilagsdalsfjall material. Using the global dataset allows us to populate the DM-EM array with more depleted (but still variably enriched) eruptions, i.e., those which are less susceptible to CO 2 loss by magmatic degassing. Figure 10 shows the compilation (filtered for eruptions with R/R a ∼ 8) of CO 2 /Ba and CO 2 /Nb ratios plotted against ratios used as indices of enrichment: 143 Nd/ 144 Nd (as used in Section 5.4), 206 Pb/ 204 Pb and Ba/Nb. Though not shown here, similar observations can be made with 207 Pb/ 204 Pb, 208 Pb/ 204 Pb and 87 Sr/ 86 Sr. Some caution must be exercised with comparing the radiogenic isotopes globally, as source age is convolved with enrichment in the isotope ratio values. As discussed in Section 1, we advocate taking the highest values of CO 2 /Ba and CO 2 /Nb from melt inclusion datasets as the best proxy for the mantle, in preference to the dataset averages often used (Saal et al., 2002;Le Voyer et al., 2017;Hauri et al., 2018). In Figure 10 we highlight the difference between the two approaches for the datasets which preserve correlations between CO 2 and the trace elements (filled and unfilled blue circles). We distinguish melt inclusion datasets which do not retain a CO 2 -trace element correlation (orange circles) since they likely experienced greater extents and a more complex history of CO 2 degassing and magma mixing. CO 2 -undersaturated (ultra-depleted) MORB glasses are shown as grey circles. Figures 10a, 10b and 10c show there is no systematic variation in CO 2 /Ba with eruption enrichment. Likewise, Figures 10d and 10e show there is no covariation of CO 2 /Nb with radiogenic isotope tracers of eruption enrichment, however a correlation is observed with Ba/Nb (Figure 10f). Hirschmann (2018) also described this feature in his compilation of global CO 2 /Ba and CO 2 /Nb ratios, though he found a stronger correlation between CO 2 /Nb and Ba/Nb than shown here. This discrepancy arises since Hirschmann (2018) includes datasets associated with plumes (which we do not include as they may have contributions from high 3 He/ 4 He mantle) and presents the undersaturated glass measurements as a set of averages.
The highest CO 2 /Ba ratios observed in datasets with CO 2 -trace element correlations coincide with the highest CO 2 /Ba ratios seen in the undersaturated glasses, with the exception of the Garrett inclusions where the mean CO 2 /Ba ratio observed in the melt inclusions is closest to the undersaturated glasses. The highest CO 2 /Ba ratio measured in undersaturated glasses from the Siqueiros fracture zones (146±31, Michael and Graham, 2015) is intermediate between the mean and highest Siqueiros melt inclusion CO 2 /Ba ratios (though it is within uncertainty of the highest inclusion value). This might suggest that the melts that were preserved as melt inclusions may have experienced some CO 2 loss, prior to homogenisation of the trace element and CO 2 enriched melts, and before mixing with depleted melts to form a binary mixing array with constant CO 2 /Ba ratio (Matthews et al., 2017, Appendix A). Equally, the undersaturated glasses may reflect a subtly different mantle source. However, the coincidence of the mean CO 2 /Ba ratio of the Garrett melt inclusions with the highest CO 2 /Ba ratios observed in undersaturated glasses lends sup-  Figure 10: Compilation of data from the global spreading ridge system, filtered for R/Ra ∼ 8. Melt inclusion datasets are separated into those exhibiting correlations between CO 2 and trace elements (blue circles), and those without a correlation (orange squares). For comparison we also show the CO 2 /Ba and CO 2 /Nb ratios that were reported for the datasets exhibiting correlations; these values were derived by taking dataset averages rather than by taking the maximum values as advocated here. Glasses that are undersaturated in CO 2 vapour at the pressure of eruption are shown as grey circles. The 143 Nd/ 144 Nd, 206 Pb/ 204 Pb and Ba/Nb ratios for the melt inclusion datasets are taken from published analyses of whole-rock or matrix glass from the same eruptions, for the which the sources are listed in Supplementary Table 3 Wanless et al. (2014). The undersaturated glass CO 2 /Ba and CO 2 /Nb ratios were reported by Shimizu et al. (2016); Michael and Graham (2015). port to the argument of (Shimizu et al., 2019); that the Garrett melt inclusions represent a rare case of the inclusions preserving largely undegassed melts.
Most of the datasets without a CO 2 -trace element correlation have CO 2 /Ba ratios lower than datasets with a correlation, this could reflect a greater extent of CO 2 degassing and mixing seen by the datasets that lack a correlation. The undersaturated glasses have CO 2 /Ba ratios that overlap with the range displayed by the melt inclusion datasets. The lowest CO 2 /Ba ratios (<50) seen in the undersaturated glasses are significantly lower than the CO 2 /Ba ratios seen in melt inclusion datasets at similar Ba/Nb. Matthews et al. (2017) demonstrated that undersaturation at eruption cannot be taken as evidence that CO 2 has not been lost; the undersaturated nature of the glasses may come from mixing of melts undersaturated in CO 2 vapour with melts that saturated in CO 2 vapour during magma storage. We therefore think it most likely that the undersaturated glasses with lowest CO 2 /Ba lost substantial CO 2 vapour during magma storage.
The strength of the correlation between CO 2 /Nb and Ba/Nb shown in Figure 10f arises from the low CO 2 /Nb ratios associated with the samples with low Ba/Nb. The low Ba/Nb samples are from eruptions in fracture zones, including the melt inclusion datasets from the Garrett and Siqueiros fracture zones (Shimizu et al., 2019). Hirschmann (2018) suggests that the correlation may arise from Nb being fractionated from Ba and CO 2 during mantle melting. Whilst we have shown that Nb and Ba are not fractionated from each other during the melting processes that generated the Iceland melt inclusion suites (Section 3.2); the same is not true of Siqueiros. Supplementary Figure 18 demonstrates that the Ba/Nb ratio of the Siqueiros inclusions increases with Ba concentration, and the normalised variances of Ba and Nb differ, both observations that suggest Ba and Nb have been fractionated from each other during melting. However, CO 2 and Ba have comparable normalised variance, suggesting that they are sufficiently incompatible that they haven't been fractionated from one another and that the CO 2 /Ba ratio should therefore have been largely unaffected by the melting process. Michael and Graham (2015) suggested that ultradepleted glasses are generated where melt aggregation is incomplete, with a particularly low contribution from the deepest melts. This partial aggregation provides a mechanism for fractionating Nb from Ba and CO 2 that is likely to only occur in fracture zones, and therefore may explain the origin of the low CO 2 /Nb and Ba/Nb population in Figure 10. If fracture-zone eruptions (low Ba/Nb) are discounted from Figure 10, a correlation remains between Ba/Nb and CO 2 /Nb, though it is weaker. Given that Ba and Nb have not been fractionated from each other during melting in the petrogenesis of these non-fracture zone eruptions, this correlation is most likely a source feature and may reflect a decoupling of Nb and Ba concentrations during recycling (Hirschmann, 2018).
As discussed in Section 5.3, it is possible that the variation in 143 Nd/ 144 Nd observed between Icelandic eruptions is decoupled from CO 2 /Ba and CO 2 /Nb by melt transport, and so our proxy for source components could be invalid. A more reliable record of the local mantle average Nd-, Sr-and Pb-isotope ratios is ridge-segment average values (Shorttle, 2015). Supplementary Figure 19 nonetheless demonstrates that CO 2 /Ba and CO 2 /Nb do not co-vary with ridge-segment average Nd-and Pb-isotope ratios.

Whole mantle CO 2 mass balance
In Section 5.4.2 we demonstrated that reasonable constraints may be placed on the CO 2 concentration in a depleted mantle component (DM) and a primordial component (PM) in Iceland, but there is little constraint on the CO 2 content in the recycled mantle component (EM). In Section 5.7 we showed that there is little variation in the CO 2 /Ba ratio with degree of enrichment in the upper mantle. If the variable enrichment of the upper mantle is primarily controlled by varying amounts of recycled material (e.g. Chauvel et al., 2008) this might be taken to suggest there is little difference in the CO 2 /Ba ratio of the depleted mantle and recycled mantle.
Marty (2012) estimated a whole mantle averaged CO 2 concentration of 2130±1390 ppmw. By making the approximation that the whole mantle can be described as a mixture of the three mantle components, DM, EM, and PM, we inverted for using the Iceland data, we can assess the range of EM CO 2 concentrations that are consistent with the bulk mantle CO 2 estimate made by Marty (2012). Since the relative proportions of the three components are poorly known we allow them each to vary from 0-100%. Figure 11 shows the CO 2 concentration of the EM component required to satisfy this mass balance, as Figure 11: The concentration of CO 2 required in the EM component to satisfy mass balance assuming a bulk mantle CO 2 concentration of 2130 ppmw (Marty, 2012), plotted here in a ternary diagram of mass fraction. The calculation assumes the DM mantle component has 102 ppmw CO 2 , and the PM component has 2.2 wt%. The dashed line indicates the locus of solutions where the EM and DM components have the same CO 2 /Ba ratio, assuming a Ba concentration of 0.563 ppmw in the DM component, and 2.55 ppmw in the EM component (Table 1). a function of the relative proportions of DM, EM and PM in the whole mantle. In the calculation the CO 2 content of the DM and PM components is taken as the median value of their respective posterior distributions from Section 5.4.2. If the depleted mantle and enriched mantle have the same CO 2 /Ba ratio, 8-10% of the mantle must be primordial to satisfy this mass balance.
The large low shear velocity provinces (LLSVPs) near the base of the mantle have been suggested to host primordial geochemical signatures (Castillo, 1988), though these account for only 2% of the mantle's volume (Hernlund and Houser, 2008). Whilst the LLSVPs may be too small to host the required mass of carbon, mass-balance calculations using lithophile elements suggest up to 50% of the mantle might be primordial (Turcotte et al., 2001). Having 10% of the mantle bearing primordial high CO 2 concentrations is, therefore, perhaps not unreasonable. Further support comes from 2D mantle convection models; Ballmer et al. (2017) demonstrated that 10-15% of the mantle could remain isolated from convection if those domains are anomalously rich in bridgmanite.

Conclusions
We have leveraged a large new melt inclusion dataset of trace element and CO 2 concentrations in geochemically diverse Icelandic eruptions, alongside existing suites of melt inclusions, to place new constraints on the interplay of source and process on CO 2 /Ba and CO 2 /Nb ratios. The key results are: 1. Though there is a global covariation of CO 2 /Ba with enrichment in melt inclusions, this is a result of olivine decrepitation limiting the CO 2 concentration in melt inclusions. Decrepitation may be avoided if a significant fraction of a melt inclusion's CO 2 budget is partitioned into a vapour bubble. When the CO 2 contained in vapour bubbles is added back into their coexisting glass, substantial heterogeneity in CO 2 /Ba between Icelandic eruptions remains. 2. The maximum CO 2 /Ba and CO 2 /Nb values within Icelandic melt inclusion datasets co-vary with radiogenic isotope proxies for the contributions of depleted, primordial and enriched mantle components to the melts. The highest CO 2 /Ba and CO 2 /Nb are associated with primordial mantle derived melts, and the lowest with enriched mantle derived melts. Whilst it is likely the low CO 2 /Ba and CO 2 /Nb associated with the enriched mantle is due to the greater tendency for enriched melts to be affected by degassing, we argue the elevated CO 2 /Ba and CO 2 /Nb ratios for primordial mantle derived melts relative to depleted mantle is a consequence of elevated CO 2 in the solid primordial mantle. 3. Using the observed covariation of CO 2 /Ba and CO 2 /Nb with mantle source we invert for the CO 2 concentration in the Icelandic mantle sources, and find concentrations of 102 +53 −27 ppmw for the depleted mantle, 2.2 +3.8 −1.5 wt% for the primordial mantle, and 67 +22000 −67 ppmw for the enriched/recycled mantle. Though the CO 2 concentration is likely to be greater in the enriched mantle than the depleted mantle, this could correspond to similar CO 2 /Ba and CO 2 /Nb ratios in both sources. 4. The global mid-ocean ridge CO 2 /Ba and CO 2 /Nb observations show no systematic variation of CO 2 /Ba with enrichment, but do show a variation in CO 2 /Nb with Ba/Nb ratio; these observations may suggest similar behaviour of Ba and CO 2 during recycling, but different behaviour of Nb. 5. Our new estimates of mantle source CO 2 con-centrations are consistent with estimates of the total quantity of CO 2 in the mantle. The standard deviation normalised to the mean for each trace element is shown as a function of its partition coefficient (D) as reported by Workman and Hart (2005). CO 2 * indicates bubble-corrected CO 2 concentrations. Panel b shows the Pearson correlation coefficient for 1/El vs CO 2 /El (where El is a trace element), as used by Matthews et al. (2017). The vertical dashed line shows the partition coefficient of CO 2 during melting as determined by Rosenthal et al. (2015).

Acknowledgements
where C n and C ( n − 1) are the CO 2 concentrations in the residual porosity after, and before the nth melting increment, respectively. This can be applied recursively to yield: where C 0 is the CO2 concentration in the melt at the point of graphite exhaustion. Setting dF n+1 = dF n = dF , Examples of these calculations for variable fO 2 , are shown in Figure B.2. Graphite saturation during melting decreases the CO 2 content of the most enriched melts, and increases the CO 2 concentration in the more depleted melts. Since the degassing step substantially reduces the CO 2 concentration of the most enriched melts in most runs of the model, the most pronounced effect of graphite saturation is on the depleted melts. Depleted melts generated from a very reduced melting region have extremely high CO 2 /Ba ratios ( Figure B.3). Dissolution of CO 2 , exsolved from enriched melts, into undersaturated depleted melts will have the same qualitative effect on the CO 2 -trace element systematics. Though graphite-saturated melting is unlikely to be important unless the mantle is significantly more reduced than indicated by Fe 3+ /ΣFe ratios of Icelandic basalts , or MORB (Cottrell and Kelley, 2011;Berry et al., 2018), it is feasible that extremely high mantle CO 2 concentrations may exert an influence on redox equilibria. Whilst the CO 2 -trace element systematics of Háleyjabunga and Miðfell rule out equilibrium with graphite for most of the melting region, it is not implausible that the most extreme values of CO 2 /Ba could be generated in this way. Significantly more melt inclusion analyses are required to fully characterise the high-tail of the CO 2 /Ba distributions. Figure B.2: Panel a shows the melt fraction as a function of pressure, calculated as described in the main text. The solid blue lines in panels b-f show the calculated CO 2 concentration of melts generated at variable oxygen fugacity, buffered relative to the FMQ buffer. Below the grey line the mantle is graphite saturated. The red dashed line shows the concentration of CO 2 in the melts in a mantle sufficiently oxidised that graphite is never saturated during melting, i.e. the result used in the main text of this paper. Figure B.3: Results of the mixing and degassing model for melts produced during graphite saturated melting (Error! Reference source not found.) and degassing at 2000 bar. The horizontal red and green lines show the source CO 2 /Ba ratio, and the ratio inferred from orthogonal-distance regression, respectively. The shading shows the degree of CO 2 undersaturation in bars, at the pressure of degassing