SHV products are phenocryst-rich (30–45 vol%), with an assemblage of plagioclase, amphibole, orthopyroxene, titanomagnetite and quartz (<0.5%), and minor amounts of clinopyroxene occurring as microphenocrysts or as overgrowth rims on orthopyroxene, plus apatite and ilmenite (Humphreys et al., 2009; Edmonds et al., 2014). Petrological features such as mineral phases and enclave textures are similar throughout the eruptive phases, as described in Christopher et al. (2014). Whole rock compositions from all phases are largely andesitic and range from 57–64 wt% SiO₂ (Murphy et al., 2000; Zellmer et al., 2003; Plail et al., 2018), while groundmass glass compositions, published for Phases 1–3, are 70–80 wt% SiO₂ (Edmonds et al., 2002; Harford et al., 2003; Buckley et al., 2006; Humphreys et al., 2010).

Geochemical and geophysical studies point to both a two-tiered magma storage region (e.g., Aspinall et al., 1998; Barclay et al., 1998; Devine et al., 2003; Rutherford and Devine, 2003; Foroozan et al., 2010; Miller et al., 2010; Christopher et al., 2014; Edmonds et al., 2014) or a transcrustal mush system at Soufrière Hills (e.g., Edmonds et al., 2016; Gottsmann et al., 2020; Alshembari et al., 2024). Published measurements of dissolved H₂O and CO₂ contents within plagioclase and quartz-hosted melt inclusion glasses are 4.07–5.05 wt% H₂O and <60 ppm CO₂ for Phase 1 inclusions (Barclay et al., 1998), and ≤6.40 wt% H₂O and ≤546 ppm CO₂ for Phase 3 inclusions (Edmonds et al., 2014). These translate, via the solubility-pressure model VolatileCalc (Newman and Lowenstern, 2002), to pressure and thus depth estimates of ∼130 MPa, and 5–6 km for Phase 1 (Barclay et al., 1998) and ≤300 MPa for Phase 3 (Edmonds et al., 2014) which is equivalent to ∼7–11 km. Edmonds et al. (2014) attributes the higher CO₂ content of a few plagioclase- and orthopyroxene-hosted inclusions of 836 and 1032 ppm to CO₂ flushing, which occurs when CO₂- rich fluids are released from deeper in the magmatic system, in shallow conduit systems, or from carbonate sources, and interact with magmas stored in the upper or mid-crust (e.g., Rust et al., 2004; Marianelli et al., 2005; Spilliaert et al., 2006; Blundy et al., 2010; Caricchi et al., 2018).

However, these CO₂ measurements neglect melt inclusion bubble-hosted CO₂ and are thus likely severely underestimated (e.g., Hartley et al., 2014; Moore et al., 2015; Wieser et al., 2021). The estimated magma storage depths of 5–6 km are similar to earthquake hypocentral depths and seismic tomographic data (Aspinall et al., 1998; Miller et al., 2010). Mineral geochemistry also yields shallow storage depths (around 5–6 km), via Al-in-hornblende geobarometry (Rutherford and Devine, 2003) and clinopyroxene-melt equilibria (Christopher et al., 2014). However, iron oxide compositions (Devine et al., 2003) point to deeper storage regions >10 km, along with H₂O contents of enstatites (Edmonds et al., 2016), which indicate a magma storage region that is vertically elongated through the crust.

A two-tiered model is suggested from geodesy, where best-model fits to GPS data from Phase 1 identify a source at ∼6 km depth (Mattioli et al., 1998), seemingly switching to a deeper-fed region at 10.4 ± 2.1 km during the later phases (Mattioli et al., 2010). Foroozan et al. (2010) suggests 5 km and 17 km deep storage regions, also based on GPS data. More recently, geodetic modelling using 3D crustal mechanical and topographical data has proposed the presence of a vertically extended pressure source between ∼4 and 14 km depth (Gottsmann et al., 2020).

Bubbles in melt inclusions were analysed at the Laboratoire Magmas et Volcans (LMV), Clermont-Ferrand, France using a Renishaw inVia confocal Raman microspectrometer. This was equipped with a 532.1 ± 0.3 nm diode-pulsed solid state laser, a Rayleigh rejection edge filter (cut-off at about 50 cm⁻¹), and a CCD detector of 1040 × 256 pixels. For each analysis, a slit aperture of 20 µm (high confocality setting) and a grating of 2400 grooves/mm was used. A Leica DM 2500 M optical microscope with a motorised XYZ stage was used to focus samples, and ×50 or ×100 microscope objectives were used, dependent on bubble size. A spectral resolution better than 0.4 cm⁻¹, and spatial resolutions of few µm were achieved based on the applied conditions. To calibrate peak positions and check the linearity of the spectrometer, the 520.5 cm⁻¹ peak of Si and the two neon emission bands (568.982 and 576.442 nm) were used, as the neon bands bracket the peaks of CO₂, known as the Fermi diad (Δ, peaks at ∼1388 and 1285 cm⁻¹). In order to acquire CO₂ spectra used for quantification of CO₂ concentration, spectra were collected in a single window ranging from ∼725 to 1880 cm⁻¹, using the WiRE™ 4.4 software. Each measurement took 120 s (3 acquisitions of 40 s). Neon bands were measured before and after each analysis of CO₂ and the correction factor (^realΔ_Ne/^measuredΔ_Ne) for each measurement lies between 0.9987 and 1.0003. In order to quantify CO₂ concentration in the bubbles, fluid inclusions of pure CO₂ with known densities were used as standards and were analysed three times during each analytical session (Boudoire et al., 2023). Uncertainties associated with the reproducibility of the measurement determined on standards are <0.04 g/cc. The 60–1320 cm⁻¹ wavenumber range was subsequently examined for the identification of mineral phases known to occur in the bubble (sulfates, carbonates, etc.). The presence of other fluid species (i.e., liquid or vapour H₂O, H₂S, HS⁻) was also investigated by measuring the whole spectral range up to 4,000 cm⁻¹.

The concentration of H₂O and CO₂ in melt inclusions glasses were measured at the Natural Environment Research Council (NERC) Ion-Probe Facility at University of Edinburgh, UK using a Cameca IMS 7f-GEO paired with a 5 nA ¹⁶O⁻ beam.

Prior to analysis, each sample was pre-sputtered for 180 s ²⁴Mg²⁺, ²⁶Mg, ³⁰Si (counting times = 2s), ¹H (counting time = 1s) and ¹²C (counting time = 10 s) were analysed over 10 cycles with an electron multiplier. In order to separate the mass interferences of ²⁴Mg²⁺ and ¹²C, a mass resolving power of 2000 was applied. The curves of ¹H/³⁰Si vs. H₂O and (¹²C/³⁰Si)*SiO₂ vs. CO₂ for H₂O and CO₂ respectively, were used for calibration based on a set of known glass standards (Supplementary Figure S1; H₂O = 0.64–7.56 wt%; CO₂ = 0–10,380 ppm). Eight of fifteen orthopyroxene-hosted, and one of sixty-one plagioclase-hosted inclusions required calculation of H₂O by difference, due to high measured H₂O contents which exceeded that of the standards (7.89–8.76 wt%), leading to high total oxides of 102–104 wt%. H₂O was calculated to achieve totals of 100.44 wt%, the average total for the remaining orthopyroxene-hosted inclusions where the standard deviation is 0.54 wt%. Calculation of CO₂ required SiO₂ which was measured via EPMA. H has a matrix correction and does not require further correction. Background concentrations for H₂O and CO₂ were measured on nominally anhydrous minerals (plagioclase and orthopyroxene) before final concentration calculations. Pressure in the sample chamber was <6.80 × 10⁻⁸ mbar over the analytical session. Reproducibility (2σ) on known standards amounts to <10% for both H₂O and CO₂, with a detection limit of 0.003 wt% for H₂O, and 3 ppm CO₂.

Melt inclusion glass and host major elements were analysed at the University of Cambridge utilising a JEOL JXA-iHP200F HyperProbe with 15 kV accelerating voltage.

Major elements in the melt inclusion glasses were measured with a beam size of 5 µm. A 5 µm and a defocused beam were respectively applied for major element analysis of plagioclase and orthopyroxene host crystals. Beam current of 6 nA was applied across all analyses, and alkalis were measured first to mitigate for loss or migration.

Reproducibility (2σ) of major elements is based on repeat measurements of known rhyolitic glass standard AthoG and are <5% for SiO₂, TiO₂, Al₂O₃, FeO and CaO, <10% for K₂O, <13% for MgO and Na₂O.

Counting times for each element along with standards and associated diffraction crystals used for glass analysis are listed in Supplementary Table S1.

Whole rock samples from Phases 1-5 were analysed for major element composition at LMV, France using an Agilent 5800 ICP-OES instrument.

An induction furnace was used to melt 100 mg of each sample together with 300 mg of LiBO₂. The resulting product was dissolved in 1M HNO₃ until a final volume of 200 mL was achieved. “GH” and ‘BR’ from Centre de Recherches Pétrographiques et Géochimiques (CRPG), Nancy, France were used as standards for Si, Na and K, and Al, Ti, Fe, Mn, Mg, Ca and P respectively (Supplementary Table S2). The errors on reproducibility of the standards are <10% (2σ).

Post-entrapment modification of melt inclusion compositions is common, and occurs via diverse processes, resulting in compositions not representative of the parental melt. These processes include post-entrapment crystallisation that modifies both major and volatile elements (Anderson and Brown, 1993; Danyushevsky et al., 2000; Kent, 2008), the formation of bubbles that can be empty (Schipper et al., 2010; Steele-Macinnis et al., 2011), or contain vapour (Anderson and Brown, 1993; Moore et al., 2015; Moore and Bodnar, 2019), or aqueous species and solids (Schiavi et al., 2020). Other processes include diffusion of H⁺ into and out of inclusions (Gaetani et al., 2012), altering H₂O contents, and decrepitation leading to volatile loss (Neave et al., 2017). As these processes allow for misrepresentation of the major element and volatile composition of the melt, it is important to assess their extent in individual inclusions, and make corrections where possible, prior to further modelling.

The major element composition of seventy-seven SHV melt inclusions across phases 1, 2, 4, and 5 is plotted against K₂O as a representation of magma differentiation (Figure 5). There are trends in the overall dataset, and the dataset can also be divided into distinct groups, where the compositions differ based on Phase and host mineral. Overall, SiO₂ ranges from 69.63–78.15 wt% with total alkalis of 4.48–7.51 wt% across K₂O 1.61–2.84 wt%, classifying the inclusions as dacite and rhyolite based on total alkali vs. silica (Figure 2; Le Bas et al., 1986). Inclusions hosted in orthopyroxene are less evolved than those hosted in plagioclase, and inclusions are all more evolved than whole rock compositions. However, they overlap with groundmass glass compositions of 69.85–80.03 wt% for Phases 1–3 (Edmonds et al., 2002; Harford et al., 2003; Buckley et al., 2006; Humphreys et al., 2010). Increasing trends with differentiation are seen in SiO₂ and Na₂O vs. K₂O, while decreasing trends occur in CaO vs. K₂O. Trends remain relatively constant for TiO₂, Al₂O₃, FeO and MgO vs. K₂O.

In Phase 1, SiO₂ ranges from 70.06–74.04 wt%, and is one of the least evolved Phases. Major element data exist only for plagioclase hosted inclusions in Phase 2. This group generally overlaps with Phases 1 and 5 with SiO₂ of 70.36–75.37 wt%. Phase 4 stands out as the most evolved group, with SiO₂ ranging from 72.14–78.15 wt% in plagioclase hosted inclusions. Inclusions in Phase 5 are similar to those of Phase 1, with 69.63–74.89 wt% SiO₂. Overall, there is a trend of increasingly evolved melt compositions through Phases 1–4, before compositions revert to lower SiO₂ contents during Phase 5. Similar temporal evolution–with a distinction between the Phase 1-4 trend and Phase 5 – is also identified in CaO and Na₂O (Figure 5).

SHV melt inclusions are rich in volatiles, with H₂O in the glass of seventy-five orthopyroxene-hosted inclusions ranging from 5.38–7.74 wt%, and 2.47–6.40 wt% in plagioclase-hosted inclusions (Figure 7A) throughout the eruption. Melt inclusion H₂O values across all Phases are similar to those recorded for Phase 1 and 2 plagioclase- and quartz-hosted inclusions (1.20–6.86 wt%, reported by Barclay et al., 1998; Mann et al., 2013), and are also similar to values reported by Humphreys et al. (2009) and Edmonds et al. (2014) where H₂O in plagioclase- and orthopyroxene-hosted inclusions measured 0.03–6.40 wt% for Phase 3 samples. Values from this study are also roughly similar to H₂O measured in orthopyroxene cores (6-9 wt%) for Phase 3 data of Edmonds et al., 2016.

H₂O across the eruption exhibits a decreasing trend with K₂O. This trend is also apparent in individual Phase groups, where plagioclase-hosted inclusions indicate degassing associated with crystallisation, except in Phase 1, where most of the data fall between a very restricted K₂O range, but are consistent with the overall degassing trend. Separating the dataset into Phase groups illustrates not only the degassing trend in the Phases, but also the difference in H₂O content as the eruption progresses. Similar to major elements, there is a decrease in H₂O contents from Phase 1, at the beginning of the eruption, to Phase 4, with a return to high H₂O in Phase 5 (Figure 7A). Based on their relationship with K₂O, the conclusion is made that SHV inclusions experienced at most negligible H₂O loss, as inclusions do not fall out of trend at a given K₂O.

This distinction in Phase groups is however not seen with melt inclusion glass CO₂ concentrations, where CO₂ across the eruption does not trend with K₂O except for Phase 2, and are typically <1000 ppm, with 67 of 69 inclusions containing 13–937 ppm, while two inclusions which contain 1126 (SHV_P5_PLAG_002_MI4) and 1243 (SHV_P5_PLAG_042_MI2) ppm.

Bubbles at SHV contain CO₂ vapour along with solid phases crystallised on the bubble wall, identified via their main and other vibrations (Frezzotti et al., 2012; Figure 6). Solids are in the form of sulphates and are normally anhydrite and gypsum. No carbonates were observed.

Two groups can be identified from a total of 577 bubbles analysed by Raman spectroscopy. Group A consists of 91% of the bubble population, and are described as displaying a Fermi diad. Group B accounts for 9% of the bubble population, and their Raman spectra lack observable Fermi diads, indicating no CO₂ or CO₂ with very low densities. The 523 bubbles displaying Fermi diads have a CO₂ density range from 0.001 to 0.22 g cm⁻³, with an upper limit of 0.13 (Figure 8B) for those with glass analysed via SIMS. Overall, bubble CO₂ density exists below the critical density of CO₂ (0.468 g cm⁻³; Moldover, 1974), therefore, CO₂ exists solely as vapour and is not underestimated due to the presence of aqueous CO₂. In order to correct the total melt inclusion CO₂ to account for the bubble-hosted component, we take the mass of CO₂ in the bubble and add it back into the entire mass of glass in the MI and must thus account for their relative volumes. The contribution of the bubble-hosted CO₂ to the total inclusion CO₂ is based on the mass ratio between the bubble and the glass, and is therefore given by Equation 1: C O 2 b u b b l e = ρ C O 2 V b u b b l e ρ g l a s s V g l a s s x   10 6 (1) where [CO₂]_bubble is the amount that the concentration of CO2 in the entire inclusion will be corrected when CO₂ in the bubble is accounted for. CO₂ density, volume of the bubble(s) and volume of the total inclusion are calculated according to Section 4.2.1, and glass density is calculated for individual inclusions using DensityX (Iacovino and Till, 2019). Important sources of errors in calculating the amount of CO₂ that is contributed from the bubble are the estimation of bubble and inclusion volumes. For bubbles in the range 3–6 µm (majority orthopyroxene-hosted inclusions), volume errors are as large as 24%. However, the contribution of CO₂ from the bubble to total CO₂ in these inclusions are small in comparison to the concentration of CO₂ in the glass, and therefore the total error is small (∼11%). For bubbles over 6 μm, errors on volume calculations are 5%–10%. Additionally, large errors are introduced in calculation of the CO₂ density, especially at low densities (Supplementary Data Sheet), calculated using the densimeter equation of Lamadrid et al. (2017) in the absence of a calibration curve specific to the Raman spectrometer used. Overall, errors on the calculation of total CO₂ amount to ∼22%. Based on (i) CO₂ existing solely as vapour in the bubble(s), and (ii) carbonates being absent in the bubble, the contribution of CO₂ from the bubble of seventy-two inclusions for which there are SIMS and EPMA data ranges from 9 to 3767 ppm. However, for further calculations and models, we use bubbles that do not show signs of heterogeneous entrapment (bubbles <5% of total inclusion volume). The contribution of CO₂ from these bubbles to the total inclusion is 9–3145 ppm.

The summation of CO₂ in the glass and CO₂ in the bubble (5% threshold applied) gives a total CO₂ budget of melt inclusions of 131–3230 ppm, with an average of 1006 ± 684 ppm at SHV (Figure 7C), the highest for the Northern Group islands where maximum CO₂ is 720 ppm for St Kitts (Melekhova et al., 2017), and is most comparable with islands in the Central Group where CO₂ values extend up to 1507 and 4012 ppm for bubble-free melt inclusions from Guadeloupe and Dominica respectively (Balcone-Boissard et al., 2018; d’Augustin et al., 2020; Metcalfe et al., 2022).

While the bubble comprised ≤5% of total inclusion volume, the percentage of CO₂ sequestered to the bubble amount to 5%–76% for inclusions hosted in orthopyroxene, and 37%–99% for those hosted in plagioclase (Figure 8A), in agreement with published data on CO₂ loss to the bubble (e.g., Hartley et al., 2014; Wallace et al., 2015; Moore and Bodnar, 2019; Venugopal et al., 2020). It is noted that orthopyroxene-hosted melt inclusions trapped less evolved and less differentiated magma and therefore contain more H₂O than the plagioclase-hosted melt inclusions. Following this trend, it is expected that these inclusions should also contain higher CO₂ contents, but this is not the case, and we are unable to adequately explain this phenomenon.

H₂O and total CO₂ do not show simple open- or closed-system degassing trends overall or when separated by Phases (Figure 7B). However, CO₂ during Phase 1 appears 1082–1685 ppm lower than Phases 2, 4, and 5. The case can be made for CO₂ degassing in each Phase as large variations in CO₂ exist at relatively constant H₂O while being hosted in phenocrysts with similar enstatite and anorthite contents. Due to its low solubility, CO₂ begins to exsolve at high pressures and is the first volatile species to be degassed (Best, 2013).

The temperatures calculated represent magmatic temperatures at the time of melt inclusion trapping, and account for PEC. For orthopyroxene-hosted inclusions from Phase 1 and Phase 5, the temperature calculation was based on orthopyroxene-liquid equilibria using Equation 28a of Putirka, (2008). Temperatures from the two Phases are not highly varied, with a total range of 880°C–916°C, and a standard error of estimate (SEE) of 28°C. Phase 1 and 5 results are not significantly different, being 880°C–898°C for Phase 1 (three melt inclusions) and 884°C–916°C for eight Phase 5 inclusions. Plagioclase-liquid equilibria are applied to fifty-two plagioclase-hosted inclusions where H₂O was measured via SIMS, to reduce overestimation of temperature due to its strong dependence on H₂O. Applying Equation 24a of Putirka, (2008) yields a narrow temperature range across Phases, with temperatures ranging from 828°C to 893°C and SEE of 36°C. However, the plagioclase melt inclusion data do reveal a significant difference in temperature between Phases, where Phase 2 inclusions yield the hottest temperatures from 875°C to 893°C along with Phase 4 (849–891). Phases 1 and 5 are cooler, with temperatures of 828°C–881°C and 836°C–864°C respectively (Figure 9). While Phase 5 inclusions are more primitive than Phase 2 and 4, based on SiO₂, FeO and MgO contents (Figure 5), their temperatures are lower. This may indicate a lack of hotter, more mafic magma input, correlating with evidence from uranium-series disequilibria in mafic enclaves and their hosts that suggest the intrusion of mafic magma had halted by Pause 2 (McGee et al., 2019).

Temperatures calculated in this study are the first temperature estimates for Phases 2, 4 and 5 using melt inclusion data and are similar to those of Phase 1 (812°C–891°C) obtained via geothermometry and experimental phase equilibria (e.g., Barclay et al., 1998; Devine et al., 1998; Murphy et al., 2000).

Constraints on hornblende stability have been used to infer that the andesitic magma was originally cooler (≤840°C), before being reheated by injection of deeper basalt, which elevated the andesite temperature to ∼880°C (Barclay et al., 1998; Devine et al., 1998).

The crustal depth and thus pressure at which magma accumulates is a key influence on volcanic system behaviour, including exsolution of volatiles, crystallisation, and mineral growth (Huber et al., 2019). The mush model, whereby melt exists within a continuous crystalline framework (Marsh, 2004; Cashman et al., 2017) has been proposed for arc volcanoes including those of the Lesser Antilles, with storage regions being continuous (Dominica and Kick ‘em Jenny) or multi-leveled (Guadeloupe, Martinique and St Lucia; Metcalfe et al., 2023). Storage depths at SHV, Montserrat have been previously estimated using volatiles in melt inclusion glasses alone (e.g., Edmonds et al., 2014; Edmonds et al., 2016) and are clearly underestimated as bubble-hosted CO₂ was neglected. Recalculation is now possible with the new whole melt inclusion glass + bubble CO₂ results provided here, providing a better understanding of the plumbing system.

After accounting for PEC and bubble growth, the values of total CO₂, H₂O, melt inclusion composition, and temperature were input into the H₂O-CO₂ solubility model MagmaSat (Ghiorso and Gualda, 2015), hosted in VESIcal (Iacovino et al., 2021). MagmaSat permits calculation of saturation pressures for each melt inclusion based on specific composition-relevant pressure-solubility relations, and its usefulness is underpinned by the assumption that the melt trapped within inclusions is representative of the magma storage regions from which the crystals are derived (Cannatelli et al., 2016). Model results indicate that melt inclusion saturation pressures over the duration of eruptions at SHV are highly variable, with plagioclase-hosted inclusions yielding entrapment pressures spanning 1000–6800 bars (100–680 MPa). Pressures for orthopyroxene-hosted inclusions are more restricted in range, between 1900–3200 bars (190–320 MPa).

To correlate entrapment pressures to depths (MI_depth), Equation 2 is applied: M I d e p t h = p ρ g (2) where MI_depth (m), is based on pressure “p” (in Pascals) calculated above, Lesser Antilles crustal density “ρ” of 2660 kgm⁻³ (Christeson et al., 2008) and g = 9.81 m⁻².

Therefore, the pressures derived here indicate that SHV melt inclusions were trapped in orthopyroxenes and plagioclases at upper to mid crustal depths (5.7 ± 0.8 to 17 ± 2.5 km), with four plagioclase-hosted inclusions in the lower crust (20 ± 2.8 to 23 ± 1.3 km), and one deeper inclusion plotting at 26 ± 4.0 km. While plagioclase-hosted inclusions span the entire calculated range, orthopyroxene-hosted inclusions are restricted to 7.4 ± 0.1 to 12 ± 0.3 km, similar to depth estimates from the H₂O content of enstatites (mean of 10 km, over a total range of 4–16 km; Edmonds et al., 2016).

The depth of magma storage has been linked to 1) the depth of neutral buoyancy where the density of magma and country rock are similar (Hooft and Detrick, 1993), 2) rheological conditions of the crust (Mazzarini et al., 2010), and 3) to the location of pre-existing zones of crust weakness (Chaussard and Amelung, 2014). Alternatively, inferred magma storage depths may instead merely reflect the magmatic water concentrations (Rasmussen et al., 2022), where the maximum H₂O content of SHV magmas correlates to ∼13 km depth. However, melt inclusion data at SHV resolve depths greater than 13 km for 20% of inclusions, and are outside of the mafic-intermediate dataset range used in the study. This H₂O limit does not appear to apply to this volcanic system - perhaps due to the higher silica content of the SHV magma.

Our new melt inclusion H₂O and CO₂ data support a vertically extensive transcrustal magmatic system (Figure 9; Figure 10; e.g;, Cashman et al., 2017). The wide range of pressures and depths relate to inclusion trapping during crystal growth within a vertically extensive magma storage system that spans 5–17 km, with no distinct gaps in crystallisation depths. The transcrustal model can also explain the wide range of entrapment pressures inferred from H₂O and Al₂O₃ contents of Phase 3 orthopyroxenes (Edmonds et al., 2016), and 3D modelling of the ground deformation during Pause 5, where best solutions require melt injection into the base of a mush-based reservoir at ∼17 km depth (Gottsmann et al., 2020; Alshembari et al., 2024). Vertically extensive mush systems are the proposed magma storage types in arc environments, including the Lesser Antilles (Metcalfe et al., 2023), and are more importantly not uncommon to Montserrat, where melt inclusions from South Soufrière Hills reveal trapping pressures equivalent to 3–17 km depth using vapour-bubble-free melt inclusions hosted in olivine and clinopyroxene (Cassidy et al., 2015a; Cassidy et al., 2015b).

Phase 1 inclusions are evenly spread at 9.7 ± 0.7 to 16 ± 1.3 km, with one outlier at 7.3 ± 0.5 km. In Phase 2, depths span 7.2 ± 1.3 to 16 ± 2.7 km, with one outlier at 26 ± 4.0 km. Phases 4 and 5 also have similar depths to the first two phases at 5.7 ± 0.8 to 17 ± 2.5 and 9.9 ± 0.4 to 17 ± 1.3 km respectively, however, with outliers at 20 ± 1.2 to 23 ± 1.3 km.

Geochemical, geodetic and Fe-oxide data place a lower limit of ∼17 km for the SHV magma storage region in both dual and mush reservoir scenarios (Foroozan et al., 2010; Edmonds et al., 2016; Alshembari et al., 2024), consistent with the inferred entrapment pressure of 58 of the 63 melt inclusions characterised in this study. Thus, the overwhelming majority of the inclusions, which resolve depths <17 km, are taken as representative of a melt undergoing crystallization, degassing, and thus differentiation.