Depth of Magma Storage Under Iceland Controlled by Magma Fluxes

The compositions of volcanic materials are sensitive to physical conditions in the underlying magmatic system. When basaltic melts are saturated in olivine‐plagioclase‐augite prior to eruption, their compositions can be used to estimate the pressure at which they last equilibrated. We developed PyOPAM, an open‐source tool that runs in Python, and use this refreshed liquid‐barometer to investigate the relationship between final depths of magma storage and magma flux. We first tested PyOPAM using 312 experimental glasses compiled from literature and found that the 1σ uncertainty is 1.13 kbar (±3 km). PyOPAM was then applied to a data set of 13,400 analyses from Iceland, where suspected controls on magma flux are well constrained. Of these, 3807 analyses return robust pressure estimates, constraining final pre‐eruptive magma storage depths for 23 of the 30 Icelandic volcanic systems. Our results indicate that magma storage pressures on Iceland are linked to melt‐flux from the mantle. This finding is consistent with previous models linking storage depths and spreading rates on the global mid‐ocean ridge system. In addition, we provide clear evidence that the magma flux, rather than spreading rate alone, is the key control on the distribution of melt at spreading centers. Increased melt flux is associated with shoaling of pre‐eruptive storage depths, indicating that mantle melt fluxes dictate the long‐term stabilization of extensive magmatic storage regions at depths shallower than 10 km. Quantitative relationships between mantle melt flux and storage depths can be used to test computational models of transcrustal magmatic systems.

We aim to investigate and quantify the relationship between magma storage depths and suspected controls to advance closer to a consensus on this matter. We apply new implementation of a liquid barometer to glass and whole-rock compositions from Iceland. The intersection of the Mid-Atlantic Ridge (MAR) with a plume at Iceland results in the generation of oceanic crust with variable thickness (15-44 km, Jenkins et al., 2018) in rift zones with varying spreading rates (0-20 mm yr −1 ). Iceland, therefore, provides an excellent opportunity to disentangle the effects of spreading rate and melt supply on magma storage depths. Abundant petrological data and geophysical constraints are available for Iceland, providing information about the likely physical controls on magma storage depths (e.g., spreading rates, crustal thickness, magma-flux) which can be compared with our petrological estimates.
For mantle-plume-fed ocean island volcanoes, magma flux has been proposed as a control on magma storage distribution and compositional diversity (Geist et al., 2014). Frequently active, high-melt flux ocean island volcanoes, such as those in the western Galápagos archipelago, are able to sustain shallow magma systems, eruptive magma from storage zones as shallow as 1 km beneath volcanic centers (Bagnardi & Hooper, 2018;Bagnardi et al., 2013;Geist et al., 1998;Stock et al., 2018). While low-melt-flux ocean island volcanoes, both in the Galápagos archipelago and elsewhere, are characterized by high-pressure, upper mantle magma storage (Clague, 1987;Gleeson et al., 2020;Kahl et al., 2021).
This idea of mantle productivity and resulting magma supply controlling crustal magma storage depth and compositional variability has also been suggested for non-plume-related intraplate volcanism. For example, large, compositionally diverse volcanic complexes within monogenetic volcanic fields are suggested to reside above loci of high mantle productivity (McLeod & White, 2018). Increased magma supply results in ascent, storage, and then eruption of more evolved magma compositions than seen elsewhere in an otherwise basaltic monogenetic field (e.g., Jeju island, South Korea (Brenna et al., 2012); Dunedin Volcano, New Zealand (Baxter et al., 2022;McLeod & White, 2018)). Similar relationships have also been proposed for arc settings, where high magma flux controls the growth of magma systems (Paterson et al., 2011), and keeps the lower crust warmer and primed to facilitate rapid magma ascent to upper crustal levels, rather than capturing the magma (Cox et al., 2020), resulting in staged, transcrustal magmatic systems spanning lower, middle, and upper crustal levels (De Silva & Kay, 2018).

Geological Setting: Iceland
Iceland is one of the most active and productive subaerial volcanic regions on Earth. Many Icelandic eruptions have affected Europe, whilst the largest ones likely had global impacts (Thordarson & Larsen, 2007). Astride the MAR, Iceland is composed of multi-layered, laterally variable crust 15-44 km thick (Darbyshire et al., 2000;Jenkins et al., 2018). The upper discontinuity at ∼20 km underlies all of Iceland, while a second, deeper crustal discontinuity at 25-44 km depth is only present in specific regions and is thickest beneath the Vatnajökull ice cap. This multi-layered, anomalously thick oceanic crust is attributed to plume-ridge interactions (Jenkins et al., 2018). The mantle plume presently located roughly beneath Öraefajökull (Figure 1; Shorttle et al., 2010), increases mantle potential temperature and melt supply (White & McKenzie, 1995). Iceland's geological geometry is produced by interactions between the MOR and the plume, with the Reykjanes, West, North and East Volcanic Zones (RVZ, WVZ, NVZ, and EVZ, respectively) representing the spreading axis for the MAR (Figure 1; Árnadóttir et al., 2008;Einarsson, 2008;Gudmundsson, 2000;Saemundsson & Saemundsson, 1974).
Spreading rates have been linked to magma storage depths in numerous studies on MOR (Chen & Morgan, 1996;Morgan & Chen, 1993a, 1993b. Iceland has a spatially and temporally diverse range of spreading rates. Spreading and volcanism in northern Iceland has been predominantly confined to the NVZ rift zone for the last 6-7 million years, while activity in the south is more complex due to the eastward progression of the plate boundary over the last 15 million years (Sigmarsson et al., 2020). The continuous volcanic zone from the Reykjanes Peninsula to the north of Langjökull in the WVZ has been rifting over the last 6-7 million years (Saemundsson & Saemundsson, 1974), with the Hengill volcanic system marking the triple junction between the obliquely rifting RVZ, the WVZ, and the transform zone of the South Icelandic Seismic Zone (SISZ). Plate spreading is divided between WVZ-RVZ and the parallel EVZ, separated by the Hreppar region (Sigmarsson et al., 2020; Figure 1). Activity formerly along the WVZ-RVZ has shifted to the increasingly active EVZ (Árnadóttir et al., 2008;Thordarson & Larsen, 2007). The southward-propagating tip of the EVZ is marked by the Vestmannaeyjar system (Jakobsson, 1979;Mattsson & Oskarsson, 2005). The Torfajökull volcanic system marks the transition between the East Volcanic Rift Zone (denoted as EVZ-I), the South Iceland Flank Zone (EVZ-II) and the SISZ (Geirsson et al., 2012;Sigmarsson et al., 2020). The EVZ contains the four most productive volcanic systems in Iceland (Grímsvötn, Bárðarbunga-Veiðivötn, Katla, and Hekla) and has produced 80% of verified eruptions throughout the holocene (Thordarson & Hoskuldsson, 2008;Thordarson & Larsen, 2007). Grímsvötn and Bárðarbunga-Veiðivötn are within EVZ-I, while Hekla and Katla are within EVZ-II. Two other flank zone volcanic belts on Iceland are the Öraefajökull Volcanic Belt (OVB), a suspected embryonic rift (Thordarson & Hoskuldsson, 2008), and the Snaefellsnes Volcanic Belt (SNFS).
There are various considerations to be made when selecting equilibrium mineral-melt pairs for barometery (Mollo et al., 2013;Neave & Putirka, 2017;Neave et al., 2019;Putirka, 2008). Different minerals and different crystal populations may reflect different crystallization depths in a plumbing system, while earlier growth phases of minerals can be chemically decoupled from equilibrium melts through a number of processes. Melt inclusions may also reflect processes occurring at a variety of depths within a volcanic system. Complications can arise with this approach due to post-entrapment processes (Hartley et al., 2014;Moore et al., 2015;Maclennan, 2017;Rose-Koga et al., 2021). For investigating final depth of storage prior to eruption, the ideal approach is to use a melt barometer which is based on the constrains and pressures at which erupted liquids last equilibrated. Here we use a barometer to establish the pressure at which a basaltic carrier liquid was last multiply saturated in olivine, plagioclase and clinopyroxene.
Equations developed to parameterize this boundary include calculating mineral-melt phase equilibria based on mass balance, stoichiometry, and single-distribution coefficients (Weaver & Langmuir, 1990), empirical estimation from glass major element compositions based on models describing mineral-basaltic melt equilibria for olivine, plagioclase, and augite at 1 atm (Danyushevsky et al., 1996), determining the OPAM boundary in P-T space by parameterizing divariant equilibrium melt in a simple CaO-MgO-Al 2 O 3 -SiO 2 (CMAS) system, CMAS + Na 2 O (Walter & Presnall, 1994) from 7 to 35 kbar, CMAS + FeO (Shi, 1992), and CMAS + FeO + Na 2 O (Shi, 1993) at 1 atm. Grove et al. (1992) presented a method to locate two saturation boundaries as a function of liquid composition and pressure of crystallization from 0.001 to 10 kbar, quantitatively predicting phase boundaries for a given melt composition. Yang et al. (1996) (Y96), considered the influence of composition and pressure on OPA-saturated melts and augite composition within the CMAS system. The effects of additional components (molar proportion of Na, K, Fe, Ti, and Si) and pressure on the molar proportion of Ca, Al, and Mg in OPA-saturated melts were incorporated into three equations. Multiple linear regression determined the coefficients for these components then used for empirical equations describing liquid composition along the OPA cotectic as a function of P and T (see Text S1 in Supporting Information S1).
The approach of Y96 was the basis of many more studies using OPA-saturated melts to investigate pressure of crystallization. Michael and Cornell (1998) rearranged Equation 2 from Y96 (Text S1 in Supporting Information S1), the expression for molar proportion of Ca in relation to additional components and pressure, to calculate the pressure of given glass samples. They opted to ignore the equations describing molar proportions of Mg and Al, aware that the rearrangement of Y96 Equations 1 (Text S1 in Supporting Information S1) and 3 (Text S1 in Supporting Information S1) were highly sensitive to inter-laboratory errors. It should be noted that rearrangement of equations that describe fits of a linear regression is mathematically incorrect.
A more recent implementation of OPAM barometry from Kelley and Barton (2008) (KB08), calculates a series of liquid compositions along the OPA cotectic at 0.1 kbar increments with the Y96 parameterization. Normative mineral compositions were then calculated for these liquids, with pressure regressed against components for calculated and original compositions. From the six calculations produced, a mean was calculated as the equilibration pressure. While this method uses all three equations of Y96, the additional steps are unnecessary and increase the uncertainty of the final outcome. Voigt et al. (2017) presented a new parameterization of experimental data to estimate the OPAM equilibrium pressure, including Cr and Fe 2+ and Fe 3+ as separate individual variables in addition to those used in the Y96 parameterization.
Many of these methods, especially graphical, rely on major elements, subject to errors introduced by inter-laboratory biases (Danyushevsky et al., 1996;Presnall & Hoover, 1984;Yang et al., 1996). Many equations developed were for limited pressure ranges (Danyuskevsky et al., 1996 (7-35 kbar), Walter & Presnall, 1994, Shi, 1992, 1993 (variations of CMAS system at 1 atm)). Grove et al. (1992) and Y96 used experimental samples created at a range of pressures (0.001-10 kbar) to parameterize individual variables-the equations with which KB08 built their spreadsheet. Although no systematic errors were identified, these methods still relied on determining presence of OPA within glass and were otherwise unable to exclude non-OPA-saturated compositions. Hartley et al. (2018) presented an implementation that addressed earlier flaws, including a statistical test to screen for glass that was not OPA-saturated, and constrained uncertainties associated with the Y96 parametrization. The statistical screen was particularly valuable for applying the OPAM barometer to olivine-hosted melt inclusions, which could not be visually assessed to have formed in the presence of OPA at depth. Prior to filtering, equilibration pressures were systematically overestimated by 1.0 kbar and subject to random error of 2.5 kbar (1σ), for 157 experimental compositions we compiled from the literature. Systematic overestimation and random error were reduced to 0.32 and 1.32 kbar (1σ), respectively, with the implementation of this statistical filtering. When the Y96 parameterization was compared to the (Voigt et al., 2017) parameterization, there was less than 0.1 kbar difference between the equilibrium pressures calculated by either approach in the case of Icelandic compositions.

Barometer Implementation
For our application of OPAM (Baxter & Maclennan, 2022b), we followed the implementation of Hartley et al. (2018), which used a statistical test for screening melt compositions and, if they were three-phase-saturated, assessing the robustness of calculated pressures. The first step of this implementation calculates the χ 2 misfit between measured X Mg , X Ca , and X Al content and those predicted by Y96. The misfit χ 2 is calculated using the following standard expression: where X i are the measured cation fractions of Mg, Ca, and Al, Y i are those predicted from the Y96, and σ is the analytical uncertainty on glass measurements by EMPA, estimated as ±5% for major element oxide components (Neave et al., 2015). Equations 1-3 from Y96 were then used to estimate the Y i as a function of pressure and other measured sample compositional properties (X Si , X Fe , X Ti , X K , X Na ) and functions 1-3 from Table 2.15, Yang et al. (1996) were used to return the pressure that minimized the χ 2 misfit. The quality of fit between the observed and predicted melt compositions is evaluated using the minimum χ 2 . The resulting significance estimate, termed by Hartley et al. (2018) as the probability of fit (P F ), allows for the melt compositions to be assessed for three-phase saturation, and, thus, the calculated OPAM equilibration pressures to be assessed for robustness. If melt compositions deviated from three-phase saturation and, therefore, are inappropriate for OPAM barometry, the low probability of fit returned allows for these samples to be discarded from further investigation with the barometer.

Test Data Set
We validated our pyOPAM code with a test data set of 312 experimental glasses compiled from the Library of Experimental Phase Relations (LEPR; Hirschmann et al., 2008) and literature (Baxter & Maclennan, 2022c).
This includes experiments carried out over a range of pressures (0.001-14 kbar) on tholeiitic to alkaline basalt compositions (1.01-9.5 wt% alkalis; Figure 2; Table 1). Glasses were only included in our test data set when observed alongside phases of olivine, plagioclase, and clinopyroxene (± one of spinel, magnetite, or ilmenite).

Barometer Validation
Of the 312 OPA-saturated glasses, 62 contained ≤4 wt% MgO, indicating that these compositions were inappropriate for testing the barometer as applied here; these samples were removed before the remainder of the experimental glasses were processed. The remaining 250 glasses were henceforth considered unfiltered OPAM results. The pressure estimates from our implementation for the unfiltered glasses produced a mean absolute error (MAE) of 2.07 kbar (Table 1; Figure 2a). Absolute difference between predicted and known experimental pressures is ΔPr (we use Pr to disambiguate it from P for probability). ΔPr for 110 experimental glass compositions that pass the probability of fit (P F ) threshold set to 0.7, ranges within 3.9 kbar (one exception at P F = 0.84 of ΔPr = 4.7; Figure 2b). These return a MAE of 1.28 kbar (Table 1). There is little difference in the MAE returned when P F ≥ 0.8 and P F ≥ 0.9, returning 1.26 and 1.23 kbar, respectively, but 28% more experimental glasses failed the filter at P F ≥ 0.9 (Table 1). We retain the P F value set by Hartley et al. (2018), which was as such to eliminate false positives from passing into filtered datasets. Henceforth, our filtered results are analyses which passed P F ≥ 0.8. This filter is retained for our subsequent treatment of the Icelandic data set. In comparison, unfiltered compositions processed with the (Kelley & Barton, 2008) (KB08) implementation have a more scattered ΔPr (Figure 2e) than our unfiltered results, producing an MAE of 3.77 kbar, and our filtered results when processed with the KB08 implementation produced an MAE of 2.86 (Table 1).
Filtered samples only include tholeiitic to alkali basalts, with total alkalis ranging 1.0 to 5.0 wt % ( Figure 2d). From the samples containing either spinel, magnetite or ilmenite (in addition to OPA), only samples with spinel passed the filtering. There is no apparent relationship between ΔPr and alkalinity (Figures 2a and 2b), although pressures returned from alkaline samples are more likely to be underestimated ( Figure 2a). No sample with ≥6 wt% total alkalis exceeded a probability of 0.8 ( Figure 2.1b), perhaps indicating the limits of applicability of the Y96 parameterization. Y96 developed their parameterization for MORB tholeiitic basalts, and  noted that OPAM was unsuitable for transitional and alkali basalt compositions. While our findings indicate that it is possible for some alkaline compositions to produce reasonable OPAM estimates, we highlight the possibility for systematic errors at total alkali contents of ≥5 wt%.

Scrutiny of Barometer and Effectiveness of Our Statistical Filter
Ideal samples for OPAM barometry should contain glass and olivine, plagioclase, and clinopyroxene crystals, demonstrably in equilibrium with the glass. Assessment of equilibrium will involve petrographic checks not only for the presence of the OPAM phases but also for crystal morphologies that are indicative of crystal-melt equilibrium. Furthermore, the use of expected crystal-melt exchange coefficients (e.g., Blundy et al., 2020;Saper et al., 2022 for Mg-Fe in olivine, Mollo et al. (2013) for clinopyroxene, Neave and Namur (2022) for plagioclase) can be used to verify multiple saturation and mitigate against the effects of entrainment of a disequilibrium crystal cargo. Microanalysis of crystal-free, clean matrix glasses provide the best estimates of carrier liquid compositions. Aphyric or glassy samples may still represent the liquid, even if analyzed and reported as whole-rock X-ray fluorescence (XRF) analyses. If all crystals formed during in situ solidification and are have not accumulated in the sample by motion relative to the carrier liquid, then whole-rock compositions can represent liquid compositions. A risk with whole-rock analyses is that any antecrysts or xenocrysts present can be incorporated into the composition. The majority of bulk-rock compositions are, thus, inappropriate for OPAM barometry as they reflect a combination of the liquid and the accumulated crystal compositions (Passmore et al., 2012). For our study using literature data, information on phases present with glasses was not always reported. Additionally, for some sites around Iceland, there are little to no glass analyses available.
To test barometer performance, we investigated how crystal accumulation changes the OPAM pressure estimates, and the effectiveness of the probability filter. With experimental glasses formed at known pressures, we calculated the effect of incrementally adding a proportion of crystals to the glass composition. Using sample Y0204 formed at 3 kbar from Neave et al. (2019), we demonstrate the effect of adding 1%-40% crystals of olivine, plagioclase, and clinopyroxene, and a combination thereof, into the glass composition (Figures 3a and 3b). Olivine results in pressure overestimation. No glass + olivine combinations modeled passed the filter set to P F ≥ 0.8 ( Figure 3b). Glass + clinopyroxene crystals or glass + plagioclase crystals modeled combinations did pass P F ≥ 0.8, as long as crystal input did not exceed 2% (Figure 3b). Additional plagioclase or clinopyroxene results in underestimation of OPAM depth. Few model combinations return robust OPAM depth estimates, with glass (98%) + plagioclase (1%) + olivine (1%) and glass (97%) + olivine (1%) + plagioclase (1%) + clinopyroxene (1%) being the only two combinations which pass our P F filter for this sample. These models returned pressures within 0.6 kbar of experimental pressure, within the range of error generated by glass-only samples.
Crystal input percentage and combinations were tested for a range of pressures (Figures 3c and 3d). Samples modeled to contain ≤5% clinopyroxene or ≤3% plagioclase or ≤1% olivine added to glass compositions produced robust estimates. Successful combinations modeled rarely exceeded 98:1:1 mixtures of glass + olivine + plagioclase or clinopyroxene, or 97:1:1:1 mixtures of glass + olivine + plagioclase + clinopyroxene. The only exception was a 96:2:2 mixture of glass + olivine + plagioclase. These models returned pressures within 0.8 kbar of experimental pressure, within the range of error generated by glass-only samples. The probability filter excludes all other combinations and concentrations. The MAE for modeled mixtures that passed the probability filter when set to 0.80 is 1.19 kbar lower than the error generated for all glass samples in our test data set. This demonstrates that whole-rock samples that do not reflect a true magmatic liquid are filtered out by this implementation of the OPAM barometer. Even if the samples consist of a melt and a cotectic assemblage, these will be filtered out by the probability of fit set to ≥0.8 if they contain ≤4% of any crystal combination or deviate from three-phase saturation. If a sample is the result of two liquids mixing, and neither petrographic nor microanalysis techniques are able to identify and exclude the presence of disequilibrium assemblages and textures, pyOPAM is still able to identify and exclude a resulting non-OPA saturated melt (see Figure S2, Text S2 in Supporting Information S1).

Icelandic Data Set
We compiled ∼13,400 analyses of Icelandic samples from literature and ran these through our pyOPAM script (Baxter & Maclennan, 2022a). To demonstrate the effectiveness of the statistical filter for OPAM pressure estimates, we included a range of compositions in our Icelandic data set wider than the calibrated range of the barometer. This compilation includes both glass probe analyses and whole-rock analyses of fully crystalline, aphyric, and glassy samples (e.g., pillow lava rinds, hyaloclastites, tephra etc.). Melt inclusions are excluded from this study because our aim is to estimate the depth of magma storage immediately before eruption, and melt inclusions often sample melt that equilibrated deeper within magmatic plumbing systems. Sample localities, ages, origins, and eruption series are included in our Icelandic data set (Baxter & Maclennan, 2022a) if reported in the literature. If these details were omitted and could not be recovered, samples were removed from the data set. If lava/ hyaloclastite samples were not associated with a volcanic system of origin in the literature, they were attributed to the volcanic system within whose footprint the sample had been collected. Areas of a volcanic system footprint were taken from Thordarson and Larsen (2007). If samples could not be attributed to a specific volcanic system, due to being outside of any region of a volcanic system, or if there was uncertainty reported over the origin of tephra samples analyzed from cores, volcanic system information was omitted. If the volcanic zone of origin was still known, this information was included. Samples that could not to be attributed to a volcanic system and volcanic zone were removed from our data set. If the location name of a sample collection site was given without longitude-latitude, approximate coordinates were assigned based on the coordinates of the location named in the literature. Samples for which the origin was known (i.e., erupted from a specific central volcano) but were not collected from that volcanic footprint (i.e., tephra samples from offshore core/ice core) were allocated coordinates  (Williams et al., 2020). (e) Known experimental pressure for created glasses versus calculate pressure. Large, blue-gray points indicate pressures calculated using method given in spreadsheet from Kelley and Barton (2008)

Testing for Three-Phase Saturation
To demonstrate the effectiveness of the statistical filter for OPAM pressure estimates, we included a range of compositions in our Icelandic data set wider than the calibrated range of the barometer (Figure 2c). The first filter in our implementation removes samples with MgO ≤4 wt%, after which ∼11,588 samples remained. With the filter set at P F ≥ 0.8, 3807 samples passed, all of which were tholeiitic to alkali basalts (Figure 2d). Approximately 35% of successful samples are whole-rock analyses. In our Icelandic data set, while some studies report mineral phases observed affirming three-phase saturation of olivine, plagioclase, and clinopyroxene (Budd et al., 2016;Caracciolo et al., 2020;Eason & Sinton, 2009;Gurenko & Sobolev, 2006;, some literature reports lava composition, but not phase assemblage (Hardarson et al., 1997;Sinton et al., 2005). Furthermore, tephrochronology studies which report glass compositions do not report the presence of mineral phases at all (Bourne et al., 2015;Harning et al., 2018;Streeter & Dugmore, 2014). When mineral assemblages are not reported, we take an alternative approach in confirming compositions in our data set are three-phase-saturated.
We compare the spread of filtered samples on major element plots ( Figures S3-S5 in Supporting Information S1). Relationships between FeOt, CaO, and Al 2 O 3 with MgO wt% are controlled by the removal of olivine, plagioclase, and clinopyroxene. However, there are some outliers which deviate from this array that are unlikely to be three-phase-saturated. This can be corrected when additional filters using the minimum FeOt wt%, CaO wt% anf Al 2 O 3 wt% of the basalts used in the original (Yang et al., 1996) parameterization. We remove samples with <7.9 FeOt wt%, <8.6 CaO wt%, and <11.15 Al 2 O 3 wt%, which eliminates the remaining outliers.

Comparison Between Glass and Whole-Rock OPAM Pressure Estimates
From the 3807 analyses which passed our filtering criteria, 1354 were whole-rock samples. The RVZ, WVZ, and HVZ are mostly represented by whole-rock analyses likely due to fewer recent eruptions (Thordarson & Larsen, 2007). Some volcanic systems are only represented by whole-rock analyses. Although glass analyses are preferable, limited availability for systems with relatively few recent eruptions reduces the available glass samples. OPAM pressure estimates generated by whole-rock analyses are compared with glass-only analyses for the Þeistareykir, Bárðarbunga-Veiðivötn, Grímsvötn, Hekla, and Katla volcanic systems, and specific selected eruptions from these systems that are well represented by glass and whole-rock analyses prior to filtering (Figure 4).
For the samples of the Borgarhraun eruption within the Þeistareykir system, only glasses passed the pyOPAM filter. For the rest of the Þeistareykir system, only whole-rock analyses passed pyOPAM filters, even though both the Borgarhraun eruption and wider Þeistareykir system have numerous glass and whole-rock analyses in the unfiltered data set ( Table 2). The Vatnafjöll and Þjórsá fissure events are represented by only whole-rock compositions, partially an artifact of very few glass analyses in this data set. Þjórsá whole-rock samples produce surprisingly shallow OPAM pressures. This could be due to samples being collected from lava selvages, quenched melt that experienced some in situ equilibration and crystallization during emplacement. This highlights that field context for samples is crucial for assessing if barometry recovers storage conditions or later, post-emplacement processes. Remaining systems and individual eruptions are represented by glass and whole-rock analyses after filtering ( Figure 4; Table 2). OPAM pressures estimates from whole-rock and glass compositions were similar, with interquartile ranges overlapping for many examples (Figure 4), except for the Hekla volcanic system and Veiðivötn fissure eruptions (Figure 4; Table 2). Some examples (e.g., Grímsvötn and Laki) displayed whole-rock estimates spanning a wider range than estimates from glass, or vice versa for other examples (e.g., Eldgjá Fires), however medians from glass and whole-rock generate OPAM pressures were within the range of uncertainty (±1.26 kbar). Hekla is a striking exception, where whole-rock generates shallow OPAM pressure estimates, while glass generates a significantly deeper estimate, ∼0.35 kbar and ∼9.37 kbar, respectively.

Using pyOPAM
As with any petrological tool and modeling, there are ideal conditions for when using the OPAM barometer.
Meeting these conditions will reduce the uncertainty of the final result produced by pyOPAM. Samples for OPAM barometry should ideally be basalt glasses ( Figure 2d) containing olivine, plagioclase, and clinopyroxene crystals, demonstrably in equilibrium. This can be demonstrated through crystal morphology textures or using expected crystal-melt exchange coefficients.
We demonstrate with our large data set of diverse glass and whole-rock compositions, that pyOPAM removes non-basaltic, non-OPA saturated samples (Figures 2c and 2d; Sections 2.3 and 2.4.1). While this barometer was initially constructure for tholeiites from MORBs (Yang et al., 1996), we demonstrate it can be used for transitional alkali-basalt compositions. Samples with an alkali content greater than 5 wt% fail internal tests (Table 1.) and are removed (Figures 2c and 2d). When used on whole-rock compositions, it is best to avoid samples containing crystals out of equilibrium (see Section 2.3.1). However, when modeled composition that include cumulated crystals or other additions (antecrysts, xenocrysts) fail our P F tests are removed (Section 2.3.1; Figure 3; Figures S3-S5 in Supporting Information S1). Our modeling demonstrates that when two OPA-saturated liquids are mixed, the mixture may still pass pyOPAM filters and give an OPA-saturation depth somewhere between the two original liquids (Text S2 and Figure S2 in Supporting Information S1). However, if an OPA-saturated liquid and non-OPA-saturated composition are mixed then the resulting product typically fails the P F threshold we set.

Icelandic Magma Storage Pressures
From the ∼13,400 compositions run through pyOPAM, 3807 samples passed all filters, providing pressure estimates for 23 volcanic systems within seven volcanic zones (   (Table 3). These volcanic systems are three of the four most active volcanic systems in Iceland (Thordarson & Larsen, 2007), while the fourth, Hekla, frequently erupts intermediate to silicic compositions, not represented in this study. These four volcanic systems within the EVZ account for 80% of verified historical eruptions in Iceland (Thordarson & Larsen, 2007). This higher productivity and eruption frequency results in a bias of these systems being represented by numerous samples, while less productive volcanic systems are represented by far fewer samples (Table 3).

Volcanic Rift Zones
OPAM pressure estimates span similar ranges across all volcanic zones but have strikingly variable median OPAM pressures ( Figure S7 and Table S1 in Supporting Information S1). The EVZ-I with 1.51 kbar (IQR of 0.83-2.23 kbar) is similar to the SNVZ with 1.84 kbar (IQR of 1.05-2.28; Figure S7 and Table S1 in Supporting Information S1). In contrast to the SNVZ, the NNVZ returns a much deeper median pressure estimate of 5.27 kbar (IQR of 4.22-6.38 kbar). The WVZ has the deepest median OPAM estimate out of all volcanic rift zones, of 5.51 kbar (IQR of 2.79-8.40 kbar). The RVZ, adjacent to the south end of the WVZ, has a median of 3.89 kbar (IQR of 2.16-6.55 kbar), while the HVZ, adjacent to the north end of the WVZ, has a median pressure of 5.37 kbar (IQR of 4.66-6.77 kbar). The EVZ flank zone (EVZ-II) has a lower OPAM median than any volcanic rift zone (1.10 kbar and an IQR of 0.30-2.26 kbar), while the other flank zones are not well represented in this data set ( Figure S7, Table S1, and Text S4 in Supporting Information S1). Note. For filtered samples the minimum, maximum, interquartile range and median are presented.

Volcanic Systems
There is greater variability of OPAM pressure estimate medians and distributions between individual volcanic systems within their larger volcanic zone ( Figure 5, Table 3). The median OPAM pressure estimate for the Northern Volcanic Zone (NNVZ and SNVZ considered together) decreases southwards ( Figure 5, Table 3).
The volcanic system Fremrinámur displays a bimodal distribution of OPAM pressures, potentially marking the transition between the NNVZ (containing Þeistareykir and Krafla based on isotopic compositions (Shorttle et al., 2010)) and the SNVZ (includes Askja and Kverkfjöll; Figure 5). The volcanic systems within the SNVZ have shallow median OPAM pressure estimates, similar to those of EVZ-I. OPAM pressures are similar for the Askja system (1.93 kbar) and the adjacent Bárðarbunga-Veiðivötn (1.93 kbar) system ( Figure 5; Table 3). Kverkfjöll has the shallowest OPAM median for all of Iceland, of 0.19 kbar ( Figure 5; Table 3).
Overall, median OPAM pressures for the EVZ decrease from 1.93 kbar for the Bárðarbunga-Veiðivötn volcanic system at the northern end of EVZ-I to −0.20 kbar for Eyjafjallajökull toward the southern end of the EVZ-II flank zone ( Figure 5; Table 3). Two systems that deviate from this are off-axis Hekla-Vatnafjöll (median of 8.79 kbar) and Vestmannaeyjar (mean of 7.58 kbar), which resides at the most southern tip of the southward propagating EVZ (Figure 1). The alkaline Vestmannaeyjar system has only two samples that successfully passed the pyOPAM filter, and while these estimates fall between other barometry estimates applied to Vestmannaeyjar of 5.6-4.33 kbar (Schipper et al., 2016) and 10 kbar (Furman et al., 1991), no further discussion of this limited data set is merited here.
Volcanic systems within the WVZ have the greatest median OPAM pressure estimates in Iceland, from 9.18 kbar for Hveravellir, north of Langjökull, to 4.98 kbar for Hengill (Figures 1 and 5). The off-axis Grímsnes and Geysir systems lie geographically close to the WVZ, and are represented by few samples (15 and 7, respectively). The limited Geysir results are not be considered representative for the system. Aside from these two off-axis systems, overall median OPAM pressures decrease southwards along the WVZ toward Hengill. The Hengill system marks the transition between the WVZ and adjacent RVZ (Figure 1), and has a bimodal distribution of OPAM pressure estimates; one cluster is marked by the lower quartile of 2.42 kbar, and the other is marked by the upper quartile of 7.36 kbar ( Figure 5; Table 3). Median OPAM pressure estimates decrease for RVZ volcanic systems, from 6.42 kbar at Brennisteinfjöll to 2.61 kbar for the Svartsengi-Reykjanes system at the southwestern-most end of the Reykjanes Peninsula. Eastward of the northern WVZ, the two volcanic systems comprising the HVZ, Hofsjökull, and Tungnafellsjökull ( Figure 4) have slightly shallower median OPAM values of 5.37 and 5.58, respectively ( Figure 5; Table 3).

Magma Storage Variations Across Icelandic Volcanic Rift Zones
The pre-eruptive storage depth of basalt varies across Iceland. The depth at which basaltic magma equilibrates under rift-zone volcanoes is linked with crustal thickness and spreading rate, indicators of magma flux variation. For converting pressure estimates to depth estimates, we assume that the Icelandic crust has two distinct layers (Jenkins et al., 2018). The upper 20 km of crust is estimated to have an average density of 3000 kg/m 3 , while the lower crust is estimated to have a density of 3050 km/m 3 (Darbyshire et al., 2000). Median depth of basaltic magma storage decreases southward in the NVZ, correlating with crust thickening from ∼18 km below the Þeistareykir volcanic system, to ∼40 km below Vatnajökull (Jenkins et al., 2018), where the southern end of Askja and Kverkfjöll approach and overlap with the northern edges of Bárðarbunga and Grímsvötn (Figures 1 and 6). The NVZ spreading axis runs through the Krafla, Fremrinámur, and Askja volcanic systems, with an estimated spreading rate of 17.4 ± 0.3 mm/yr (Drouin et al., 2017;Sigmarsson et al., 2020). Median OPAM pressures increase southwards as crustal thickness decreases from ∼40 km below southern Kverkfjöll to ∼30 km below the Grímsvötn system (Jenkins et al., 2018) (Figures 5 and 6).
The relationship between spreading rates (a contributor to magma flux) and magma storage depths is evident in southern Iceland, where plate spreading is distributed over the WVZ-RVZ and EVZ rift zones. Most spreading is accommodated by the EVZ-I (Árnadóttir et al., 2008;Scheiber-Enslin et al., 2011;Sigmundsson et al., 1995Sigmundsson et al., , 2020. The rift axis is located approximately along Vatnaöldur-Veiðivötn fissure systems (Geirsson et al., 2012), spreading ∼19.0 ± 2.0 mm/yr at the northern end, decreasing to 11.0 ± 0.8 mm/yr around Torfajökull, where the SISZ intersects the EVZ. This marks the transition from rifting EVZ-I and the EVZ-II flank zone.  Figure 6). Background colors represent spreading rate changes from north to south (see text for references). Volcanic systems with split colors (Hengill and Fremrinámur), represent volcanic systems which appear to be a transition point from one volcanic zone to another. Bandwidth for creating violins was left as the default selected by the Silverman method (Silverman, 1984). Default selection was left based on the recommendation, made by Rudge (2008), to stick with default unless an adequate reason exists for deviation. Silverman's rule works well for near-Gaussian densities, whereas the Sheather and Jones method (red lines) is more flexible. Most systems show unimodal distributions and there was very little change between the two methods.  Table S1 in Supporting Information S1). These rift zone volcanic systems with shallow final storage share similar characteristics. Both the EVZ-I and SNVZ are subject to the fastest spreading rates observed in Iceland, and partially reside above thick crust generated by melting above the mantle plume conduit (Jenkins et al., 2018). This combination indicates that the overall magma flux in these volcanic systems is the highest observed across Iceland's rifting volcanic zones.
A relationship between spreading rate and final storage pressures is exhibited along the WVZ to the RVZ, where spreading rate increases from the WVZ southwards toward the Reykjanes Peninsula, as median storage pressures decrease from MOHO/lower crustal levels to mid-crustal depths (Figures 5 and 6). Spreading accommodated on the WVZ-RVZ increases southwards from ≤2.6 ± 0.9 mm/yr north of Langjökull to 7.0 ± 0.4 mm/yr along the obliquely rifting RVZ (Hjartardóttir et al., 2016;LaFemina et al., 2005;Sigmundsson et al., 1995;Travis, 2012). Crustal thickness decreases southwards toward Prestahnjúkur (Figures 5 and 6). While both northern WVZ and eastern RVZ are underlain by ∼30 km thick crust (Jenkins et al., 2018), the southwards-increasing spreading rate coincides with a decrease in median OPAM pressure estimates from the northern WVZ to the SW along the Reykjanes Peninsula ( Figures 5 and 6) Between the northern WVZ and southern NVZ, median storage pressures in the HVZ volcanic system reside between the medians for their neighboring volcanic systems (Figures 5 and 6) The HVZ experiences negligible spreading (∼1 mm/yr, Einarsson, 2008), but is underlain by crust that decreases from ∼35 km in the east to ∼25 km thick in the west ( Figure 6; Jenkins et al., 2018). The HVZ is represented by relatively few samples that pass the OPAM filter (n = 27), so results are preliminary; however, the HVZ median OPAM pressure of 5.22 kbar is greater than estimates for the SNVZ, and shallower than the northern WVZ ( Figures 5 and 6), likely reflecting decreasing magma flux from the east to the west.
Variations in the last magma storage pressures between volcanic systems across different Icelandic rift zones can be explained as a function of magma supply. There is a well-recognized connection between the depth of magma storage, spreading rates, crustal thickness, and magma supply in MOR settings, where magma flux is recognized as a key component influencing both morphology and the thickness of crust created at MORs (Carbotte et al., 1998;Macdonald, 2001;Morgan & Chen, 1993a, 1993b. At fast spreading ridges where magma supply is high, crustal thickness is less variable and a warm thermal structure enables melt lenses to form and stabilize at shallower depths (Carbotte et al., 1998;Chen & Morgan, 1996;Christeson et al., 2019;Morgan & Chen, 1993a, 1993b. These shallow melt storage regions often persist along the axis of spreading ridges. At intermediate to slow spreading ridges, the crustal thickness is variable and thermal structure is cool. Magma supply is temporally and spatially inconsistent and crystallization pressures deepen and are more variable, while ephemeral melt lenses are observed in the lower crust or upper mantle, if at all (Carbotte et al., 1998;Morgan & Chen, 1993a, 1993bWanless & Behn, 2017). Geophysical observations confirm that for slow spreading centers, where magma supply is sporadic (Bach & Früh-Green, 2010), storage zones are deep, if detected at all (Purdy et al., 1992;Searle et al., 2010). In contrast, for fast spreading ridges with stable magma supplies (Christeson et al., 2019), shallow magma storage is prevalent (Lin & Morgan, 1992).
Although variations in magma flux and spreading rate can account for much of the variability in magma storage depth along Iceland rift zones, it cannot account for all of it. The decrease in final depths of storage along the NVZ southwards coincides with increasing crustal thickness toward Vatnajökull, indicative of increasing magma flux inspite of near-constant spreading rate along the volcanic zone. The influence of the Icelandic mantle plume on melt generation must also be considered. Elsewhere, it has been suggested that a first-order control on magma distribution and diversity for mantle-plume-fed volcanism is the magma supply from plumes themselves, with higher flux rates leading to shallow crustal storage, and low-flux volcanoes being fed by deep, upper mantle magma storage zones (Bagnardi et al., 2013;Geist et al., 1998Geist et al., , 2014Gleeson et al., 2020). The increasing crustal thickness toward Vatnajökull indicates increased input from the melt generated by plume-driven upwelling (Maclennan, Mckenzie, & Gronvold, 2001) as opposed to the passive, plate-driven upwelling of hot mantle that generates 15-20 km thick crust where the rift zones intersect the coasts at Þeistareykir and the Reykjanes Peninsula. Melt production rates decrease as distance from the Icelandic plume center increases (Flovenz & Saemundsson, 1993;Martin & Sigmarsson, 2007;Sigmarsson et al., 2008;Sigurdsson, 1970). This results in the decreasing crustal thickness from the southern end of the NVZ at Vatnajökull, across the HVZ and to the northern end of the WVZ (both areas with negligible spreading rates), which in turn reflects decreasing magma flux (Figures 5 and 6). Magma storage depths increase with decreasing magma flux for volcanic zones east to west, with increasing distance from the plume center ( Figure 6). In the plume-associated crustal accretion in Iceland, rift-zone magma fluxes vary not only as a function of spreading rate but also as a function of crustal thickness, and therefore Iceland provides an ideal setting in which to test models of the link between magma flux and storage depths.

Quantifying Magma Flux Across Iceland
We quantify magma flux across Iceland by multiplying the crustal thickness and spreading rate for areas in the rift zones then dividing by the area to provide a flux in km 3 per km 2 per year (see Text S5 in Supporting Information S1, and Baxter & Maclennan, 2022a for details). This flux is compared with final storage pressures ( Figure 6). Sites located along rift zones with increased magma flux have shallower magma storage depths ( Figure 6). Slower spreading results in slower mantle upwelling beneath the ridge, lower temperatures and lower degrees of melting (Bottinga & Allegre, 1978;Reid & Jackson, 1981). Lower degrees of melting, in turn, result in thinner crust. Thinner crust, and thus lower spreading rate, influences crustal thermal structure, influencing magma storage depths (Morgan & Chen, 1993a, 1993b. In two instances, the magma flux estimates were recalculated. For one section, which includes the site of the Eldgjá Fires eruption, it has been suggested that the rifting rate in this area is 12 mm/yr (Einarsson, 2008), rather than the negligible rifting of flank zone volcanic systems. When recalculated, magma-flux estimates fall closer to the trend presented for rift zone sites. Another outlier site of interest lies between the Tungnafellsjökull and Bárðarbunga-Veiðivötn systems. This site is off-axis, so when spreading rate is reset to be midway between the on-axis spreading rate of 19 mm/yr for the BárðarbungaVeiðivötn system (LaFemina et al., 2005) and the slow spreading of the HVZ (∼1 mm/yr; Árnadóttir et al., 2008;Einarsson, 2008), the resulting recalculated magma flux for this site falls along the trend presented by rift zone sites (Figure 6c).

Flank Zones
There are limitations to the approach we have taken to investigating the relationship of magma flux and final magma storage depths. Only basaltic magma can be investigated with this barometer, excluding information recorded in other compositions. Iceland is dominated by basaltic eruptions; historical proportions of compositions for mafic: intermediate: silicic are 91:6:3%, respectively (Thordarson & Larsen, 2007). Any relationship between intermediate and silicic magma storage depths cannot, therefore, be revealed by our approach.
Volcanic systems within flank zones also have variable magma storage depths, and variable productivities. The young OVB and peripheral SNFS are responsible for 1.0% and 0.5% of historical eruptions, respectively, and 0.3% and 2.6% postglacial events (Thordarson & Larsen, 2007). The EVZ-II flank zone, however, contains two of the four most active and productive volcanic systems in Iceland (Hekla and Katla) and is a large contributor to the overall productivity of the EVZ, responsible for ≥80% of historical eruptions in Iceland. Katla and Hekla account for 73.2% and 17.7% of successful OPAM pressure estimates from EVZ-II, respectively (Table 3), yet they have vastly different median OPAM estimates ( Figure 5). Hekla is responsible for 95% of historical intermediate eruptions and is only exceeded by Öraefajökull and Askja in producting silicic eruptions (Thordarson & Larsen, 2007), compositions not suitable for study with the OPAM barometer. A large contributor to this compositional variability between Hekla and Katla and different storage depths of their basaltic magmas is the thermal structure of the crust. The geothermal gradient in Iceland decreases with increased distance from rift zones and plume center (Flovenz & Saemundsson, 1993;Martin & Sigmarsson, 2007;Sigmarsson et al., 2008;Sigurdsson, 1970). As a result, off-axis volcanic systems formed within relatively cold crust capture basaltic melt, which fractionates to move evolved compositions (Martin & Sigmarsson, 2007. Hekla, an off-axis volcanic system (Geist et al., 2021), is underlain by cooler crust, favoring the capture and fractionation of basaltic magma, resulting in more evolved compositions being finally erupted. The basaltic samples from Hekla indicate deep storage at ∼25.8 km depth, which coincides with other estimates of deep storage magma chamber depth (Geirsson et al., 2012). These represent magma batches not captured in the colder crust ( Figure 5). Magma otherwise stalled at shallower levels would evolve to compositions unsuitable for the OPAM barometer. In contrast, Katla and Eyjafjallajökull return shallow median OPAM pressures (Table 3). While little to no spreading is currently taking place along the EVZ-II, propagation of the rift continues, with the southern most tip of the rift is marked by Surtsey, Vestmannaeyjar system (Jakobsson, 1979;Schipper et al., 2015). Significant lithospheric stretching, and consequent thinning, precedes development of active rift systems (White & McKenzie, 1989) has likely begun beneath Katla and Eyjafjallajökull, resulting in a warmer thermal structure of the crust below, allowing stable shallow storage of basaltic magma.

Conclusions
PyOPAM, our new tool for estimating the pressures at which basaltic magmas last equilibrated, is best applied to samples that have been petrographically examined to ensure they contain OPA phases in equilibrium with melt. However, such assessments are not always possible. We trialed pyOPAM with a large data set of diverse compositions and demonstrate that it is a useful tool for estimating the final depths of storage for basaltic magmas, even if all ideal testing of crystal-equilibrium cannot be completed.
We applied the PyOPAM barometer to an Icelandic data set of ∼13,400 analyses compiled from the literature. Of these, 3460 samples successfully passed the filter and met the criteria for being appropriate for OPAM barometry.
Using our data set of basaltic melt storage depths across Iceland, we conclude that the contemporary melt supply rate from the mantle is the primary control over final depths of storage for basalt magma through its effect on the thermal structure of the crust. A secondary control on is linked to inherited thermal structure in flank zones or systems at the fringes of rift zones.
Low magma flux results in greater final depths of magma storage while high magma flux allows shallow basaltic magma storage to be stabilized. This is well-demonstrated in Icelandic rift zones. Median storage depths increase from the Western Volcanic Zone, where deep storage at the ultra-slow spreading northern end of the volcanic zone transitions to shallower depths by the south end of the WVZ and the Reykjanes Peninsula, where magma flux in estimated to increases. As magma flux increased southwards from the northern section of the North Volcanic Zone to the East Volcanic rift Zone, stable, shallow basaltic magma storage coincides with the location of the thickest crust in Iceland.
This approach can be applied to volcanic systems outside of Iceland. Where crustal thermal structure and magma flux can be constrained, the likely final depth of basaltic magma immediately prior to eruption can be estimated. Volcanic systems with low magma flux extract basaltic magma for eruption from deep storage zones. For volcanic systems with high magma flux, and with a warm crustal thermal profile, basaltic magma is stored at shallow depths. Even with a similar magma flux, colder crust results in deeper basaltic magma storage, as shallower basalt stalls and fractionates to more evolved compositions with the potential for producing more explosive eruptions.

Data Availability Statement
Datasets of compiled analyses from LEPR for the test dataset inspecting the uncertainty and error of pyOPAM (Baxter & Maclennan, 2022c)