Oxygen isotope evidence for progressively assimilating trans-crustal magma plumbing systems in Iceland

Magmas underneath active volcanoes are stored over a large range of depths in so-called trans-crustal magmatic systems (Cashman et al., 2017). Mantle-derived mafic magmas can assimilate the overlying crust during ascent and overprint original chemical signatures. Magma interaction with crustal rocks is best followed by observing changes in isotopic compositions. In Iceland, oxygen isotope ratios (δ¹⁸O) of crustal rocks deviate significantly from mantle values due to high-temperature interaction of low-δ¹⁸O fluids with the Icelandic crust (Gautason and Muehlenbachs, 1998). As a result, oxygen isotopes have been widely used to study the role of crustal rocks in Iceland's basalt petrogenesis, either in the form of crustal material subducted into the mantle (source contamination) and/or as contaminants throughout the magmatic column (crustal contamination) (Eiler et al., 2000a; Kokfelt et al., 2006; Thirlwall et al., 2006; Bindeman et al., 2008; Hartley et al., 2013).

Fresh mid-oceanic ridge basalt (MORB) glasses typically have δ¹⁸O values in the range of +5.5‰ ± 0.2‰ (Eiler et al., 2000b), whereas Icelandic basaltic glasses and melt inclusions (MIs) from Iceland's rift zones have δ¹⁸O as low as +2.5‰ (Breddam, 2002; Burnard and Harrison, 2005; Peate et al., 2010; Hartley et al., 2013; Halldórsson et al., 2016). Although the origin of this shift toward ¹⁸O-depleted values is a matter of debate, there is general consensus that assimilation of low-δ¹⁸O, hydrothermally altered crust (Eiler et al., 2000a; Hartley and Thordarson, 2013) and δ¹⁸O heterogeneities in the mantle source (e.g., Thirlwall et al., 2006) control δ¹⁸O variations of Icelandic glasses. However, our understanding of where and to what extent melts are affected by crustal contamination in trans-crustal magmatic systems across Iceland is limited, primarily because it is challenging to quantify the δ¹⁸O values of components that are truly mantle-derived.

Our objectives are to pinpoint the depths in the Icelandic crust at which contamination affects δ¹⁸O values of Icelandic basalts and to quantify the extent of crustal contamination as melts migrate through the Icelandic crust. We address these objectives with δ¹⁸O and trace element (TE) analyses of a well-characterized subglacial and Holocene basalt sample suite of melt inclusions (MIs) and groundmass glasses (Caracciolo et al., 2020, 2021) from the Bárðarbunga volcanic system (Figs. S1–S3 in the Supplemental Material¹). Located in the Eastern Rift Zone, the Bárðarbunga volcanic system is one of the most active systems in Iceland (Larsen et al., 2015), and it is situated above the inferred location of the Iceland mantle plume (Harðardóttir et al., 2018). The Bárðarbunga volcanic system is an ideal candidate for evaluating the effects of crustal contamination, because the crust reaches > 40 km thickness (Jenkins et al., 2018) and the plumbing architecture is likely controlled by multilevel stacked reservoirs in which melts are processed over a range of depths (Hansen and Grönvold, 2000; Maclennan, 2019; Caracciolo et al., 2020, 2021). Our new results suggest that the Bárðarbunga volcanic system is a progressively assimilating, multi-level magmatic system, and that this process is likely to occur in other parts of the active rifts in Iceland.

Oxygen isotope and TE analyses were performed by secondary ion mass spectrometry (SIMS) on MIs (n = 133) and groundmass glasses (n = 29) at the NordSIMS facility at the Swedish Museum of Natural History (Stockholm, Sweden). Also, oxygen isotope analyses (n = 16) were performed via laser fluorination (LF) at the University of Texas at Austin, USA, on groundmass glasses from the same localities (see the Supplemental Material for analytical methods). SiO₂-corrected SIMS δ¹⁸O values of MIs vary between +3.2‰ and +6.4‰, whereas groundmass glasses have δ¹⁸O values between +2.6‰ and +5.5‰, which on average are lower than those of MIs (Fig. 1; Fig. S7). Most MIs (78%) record δ¹⁸O > +4‰, while 66% of groundmass glasses have δ¹⁸O values >+4‰. In contrast, LF data of groundmass glasses, which are generally in good agreement with SIMS data (Fig. S5), record a narrower range of δ¹⁸O values, between +3.7‰ and +4.2‰ (Fig. 1A). Oxygen isotope ratios of MIs and groundmass glasses correlate with melt MgO content (Fig. 1A). Primitive MI compositions (MgO > 8 wt%) record the largest spread in δ¹⁸O values (+3.4‰ to +6.4‰) (Fig. 1A) and in TE ratios (Fig. S6), and the variability of TE ratios becomes narrower as MgO decreases (Fig. S6). The most primitive MIs preserve the most incompatible TE-enriched (Zr/Nb < 8, La/Sm > 2.2) and depleted (Zr/Nb > 15, La/Sm < 1.3) signatures, whereas the most evolved MIs and glasses record lower δ¹⁸O and intermediate TE ratios.

Evidence suggests that the Icelandic mantle is isotopically, chemically, and lithologically heterogeneous, and many studies have demonstrated that the mantle underneath Iceland contains a geochemically enriched, ¹⁸O-depleted component (Skovgaard et al., 2001; Macpherson et al., 2005; Kokfelt et al., 2006). The relationship between Sr-Nd-Pb isotope signatures and low δ¹⁸O values (down to +4.3‰) found in lavas from Reykjanes Peninsula indicates the presence of a geochemically enriched, low-δ¹⁸O mantle domain beneath the Reykjanes Peninsula (Thirlwall et al., 2006). A comparable low-δ¹⁸O component has also been documented in samples from north, south, and central Iceland (Breddam, 2002; Maclennan et al., 2003; Burnard and Harrison, 2005; Macpherson et al., 2005). The geochemical features of this enriched component likely reflect a mantle source that contains recycled subducted oceanic lithosphere (Breddam, 2002; Gurenko and Chaussidon, 2002; Macpherson et al., 2005; Thirlwall et al., 2006; Peate et al., 2010).

Collectively, MIs and groundmass glasses from the Bárðarbunga volcanic system exhibit a large variation of δ¹⁸O values (Fig. 1). Notably, melts with δ¹⁸O similar to MORB are only found in some depleted to moderately enriched primitive MIs (Zr/Nb > 10), whereas primitive enriched MIs (Zr/Nb = 7–8) have lower δ¹⁸O (Fig. 1B; Fig. S7). Assuming that depleted mantle (DM) and enriched mantle (EM) components are present underneath the Bárðarbunga volcanic system as elsewhere in Iceland (Thirlwall et al., 2004, 2006; Macpherson et al., 2005) and that variations in TE ratios, such as Zr/Nb (e.g., Fitton et al., 1997), reflect source heterogeneity, we tested to see whether our data set could be reproduced by binary mixing between the DM and EM domains (Fig. 1B). Our modeling shows that by taking into account high degrees of partial melting of DM (15%) and small degrees of partial melting for EM (5%) (Stracke and Bourdon, 2009), we are only able to reproduce a small subset of the MIs, with most of the data falling to lower values than the modeled envelope (gray field in Fig. 1B).

Decreasing δ¹⁸O values with decreasing MgO content in our set of MIs and groundmass glasses (Fig. 1A) are consistent with those of previous studies (Hemond et al., 1988; Nicholson et al., 1991; Hartley et al., 2013) and likely indicates that crustal assimilation processes play a fundamental role in controlling melt δ¹⁸O values by driving them toward increasingly lower value. Indeed, the relationship between TE ratios, δ¹⁸O, and MgO contents suggests that as melt evolution proceeds, melts acquire a low-δ¹⁸O character (Fig. 1) and the TE compositional variability collapses to a narrower range (Fig. S6) as a result of concurrent mixing and crystallization (Maclennan, 2008) coupled with assimilation of hydrothermally altered, ¹⁸O-depleted Icelandic crust.

We tested the idea of three distinct endmember components (EM, DM, and the crust) by modeling the assimilation of low-δ¹⁸O basaltic crust into mantle-derived melts by binary mixing processes. Binary mixing was modeled between the different pairs of endmembers and for different TE ratios (Fig. 2). The modeling was carried out assuming a crust with δ¹⁸O of 0‰, in agreement with δ¹⁸O values measured in drill core samples of the altered upper Icelandic crust (Hattori and Muehlenbachs, 1982), and TE ratios Zr/Nb = 11.8 and La/Sm = 1.5 (see Table 1, and the Supplemental Material).

Overall, the model shows that the distribution of magma compositions from Iceland's rift zones (Marshall et al., 2022) is consistent with the existence of DM and EM components that undergo partial melting. In particular, the δ¹⁸O values of mixtures lying along the mantle array (gray field in Fig. 1B) are then progressively shifted toward even lower δ¹⁸O values upon assimilation. Therefore, any melt composition associated with the Bárðarbunga volcanic system lies within the binary mixing lines according to a three-step process (Fig. 2):

1. Enriched and depleted mantle domains undergo partial melting, producing enriched and depleted primary melts.

2. Primary enriched and depleted melts mix in different proportions.

3. The mixed melts ascend throughout the crust. Their initial mantle-like δ¹⁸O value is lowered as they progressively assimilate low-δ¹⁸O crustal material while rising toward higher levels in the crust.

This three-step process can explain the full range of δ¹⁸O values, TE data in the Bárðarbunga volcanic system, and most melts erupted across Iceland's neovolcanic rift zones (Fig. 2; Fig. S7).

Following the mixing equation outlined in Sohn (2013) (see the Supplemental Material), we quantitatively derived the extent of oxygen isotope exchange required to explain the observed δ¹⁸O values and TE contents (Fig. 2) in each of the Bárðarbunga volcanic system MIs and groundmass glasses. The extent of oxygen isotope exchange is adopted here as a proxy for the amount of crustal material assimilated. The majority of glasses and MIs in the Bárðarbunga volcanic system require between 10% and 30% oxygen isotope exchange to explain their low-δ¹⁸O value (Fig. 3A; Fig. S9A), which is in agreement with thermodynamic limits calculated for primitive basaltic magmas (Heinonen et al., 2022). Up to 55% oxygen isotope exchange is required to explain the lowest δ¹⁸O values recorded by groundmass glasses. However, the calculated extent of exchange is strongly dependent on the chosen δ¹⁸O of the assimilant, which is difficult to constrain and likely to be heterogeneous across the crust. For example, if melts were to assimilate crust with δ¹⁸O = −2‰, we can reproduce the lowest δ¹⁸O values (around +2.6‰) with 35–40% oxygen isotope exchange (Fig. S9C).

Having constrained the extent of oxygen isotope exchange, we next seek to establish where this process occurs within the Bárðarbunga volcanic system. The equilibration pressure of glasses and MIs can be estimated by applying the Olivine–Plagioclase–Augite–Melt (OPAM) barometer (Yang et al., 1996; Hartley et al., 2018). OPAM equilibration pressures for the Bárðarbunga samples are 1.0–6.3 kbar (3.5–22.5 km), and ~60% of the samples are in the 2–4 kbar range (7.1–14.3 km) (Caracciolo et al., 2020). The distribution of equilibration pressures (Fig. 3A) is consistent with a multi-tiered, trans-crustal magmatic system in which melts fractionate and mix within an interconnected stacked sill network (Maclennan, 2019; Caracciolo et al., 2020, 2021). Equilibration pressures correlate with the extent of oxygen isotope exchange (Fig. 3A). The most contaminated melts, which record the lowest δ¹⁸O values, show the lowest equilibration pressures in the 1–2 kbar range (3.5–7.1 km). For example, melts stored and equilibrated at ~4 kbar (14.3 ± 4.7 km) have experienced 15% crustal contamination on average, but melts equilibrated at ~1.5 kbar (5.4 ± 4.7 km) have experienced 35% contamination on average (Fig. 3A). Although this correlation has a low R² (0.19), it is highly statistically significant (p < 0.0001, n = 86). The low R² value is largely due to the uncertainties associated with the OPAM barometer and the estimate of the end-member compositions. Therefore, we argue that as melts are transferred upwards throughout the Bárðarbunga volcanic system, their δ¹⁸O values lower as they become more crustally-contaminated (Fig. 3B). Our data suggest that most of the oxygen isotope exchange that affects the δ¹⁸O of melts occurs in the upper and uppermost middle crust, between 3.5 km and 10 km. Although the scarcity of data > ~4.5 kbar (~15 km) does not allow us to provide reliable estimates of the extent of oxygen isotope exchange in the deep crust, our data suggest that lower crustal material has little effect on δ¹⁸O values due to the lack of interaction of lower crustal rocks with hydrothermal and/or meteoric waters. Indeed, the lower crust lacks altered, low-δ¹⁸O rocks, and it is likely built by intrusions (Greenfield and White, 2015).

This study shows that melts supplied to the Bárðarbunga volcanic system, and possibly to other volcanic systems along Iceland's rift zones, experienced different degrees of crustal assimilation that influenced their δ¹⁸O values. We demonstrate that as melts ascend through the crust and approach shallow levels at around 3.5–10 km depth, they acquire lower δ¹⁸O values because they assimilate more crustal material. We therefore envision the plumbing system beneath the Bárðarbunga volcanic system as a progressively assimilating, multi-tiered system in which magmas are processed through a large range of depths (3.5–22.5 km) within stacked-sill reservoirs. Such complex multi-tiered volcanic systems that undergo continuous assimilation are likely common in all long-lived magmatic systems in thick crust, such as arc volcanoes and continental rifts.

A. Caracciolo was supported by the University of Iceland Research Fund (HI17060092) and by the Nordic Volcanological Center. S.A. Halldórsson and E.W. Marshall acknowledge support from the Icelandic Research Fund (grants 196139–051 and #195638–051). The involvement of S.A. Halldórsson was partly in relation to Horizon 2020 project EUROVOLC, which is funded by the European Commission (grant 731070). M. Kahl acknowledges funding by the German Research Foundation (grant KA 3532/2-1) The NordSIMS ion microprobe facility acknowledges support by the Swedish Research Council (grant 2017-00671), the Swedish Museum of Natural History, and the University of Iceland; this is NordSIMS publication 705. We thank K. Lindén for laboratory assistance in Stockholm, R. Sohn for help with mathematical equations, and J. Cullen for laser fluorination analyses. We also thank U. Schaltegger for editorial handing and J. Troch and an anonymous reviewer for helpful reviews.