Melt inclusion bubbles provide new insights into crystallisation depths and CO2 systematics at Soufrière Hills Volcano, Montserrat

Howe, TiVonne A.; Christopher, Thomas E.; Moune, Séverine; Tuffen, Hugh; Schiavi, Federica

doi:10.3389/feart.2024.1509409

ORIGINAL RESEARCH article

Front. Earth Sci., 14 January 2025

Sec. Volcanology

Volume 12 - 2024 | https://doi.org/10.3389/feart.2024.1509409

Melt inclusion bubbles provide new insights into crystallisation depths and CO₂ systematics at Soufrière Hills Volcano, Montserrat

TiVonne A. Howe¹*

Thomas E. Christopher^2,3

Séverine Moune^4,5

Hugh Tuffen¹

Federica Schiavi⁴

¹Lancaster Environment Centre, Lancaster University, Lancaster, United Kingdom
²Montserrat Volcano Observatory, Flemmings, Montserrat
³Seismic Research Centre, The University of the West Indies, St Augustine, Trinidad and Tobago
⁴Laboratoire Magmas et Volcans, Observatoire de Physique du Globe de Clermont-Ferrand, Université Clermont Auvergne, Clermont-Ferrand, France
⁵Institut de Physique du Globe de Paris, Université Paris Cité, Paris, France

Improved understanding of the magmatic system of Soufrière Hills Volcano, Montserrat (SHV) is needed to inform future hazard management strategy, and remaining uncertainties include the depth of magma storage and the source of ongoing gas emissions. Eruptive activity between 1995 and 2010 has been proposed to be sourced from either a dual chamber or transcrustal mush-based magmatic system, with volatile solubility models using H₂O and CO₂ from melt inclusion (MI) glass estimating depths of 5–6 km. To date, published SHV MI volatile data have neglected the vapour bubbles now known to sequester the bulk of MI magmatic carbon. Total CO₂ concentrations in SHV magma are therefore underestimated, together with volatile-derived entrapment pressures and inferred magma storage depths. Here, we present a new dataset of volatile (H₂O and total CO₂) and major element concentrations in plagioclase- and orthopyroxene-hosted SHV MI, that span almost all of the eruptive activity (Phases 1, 2, 4, and 5), and include the first measurement of bubble-hosted CO₂ for SHV and indeed the Lesser Antilles Arc. Analyses were conducted using Raman spectroscopy, ion microprobe, and electron probe analysis. Dacitic–rhyolitic MI occur within andesitic whole rock compositions. Volatiles in MI glass are similar to published studies (H₂O 2.47–7.26 wt%; CO₂ 13–1243 ppm). However, bubble-hosted CO₂ contributes 9–3,145 ppm, to total inclusion CO₂ with 5%–99% (median 90%) of CO₂ sequestered within bubbles, and total CO₂ concentrations (131–3,230 ppm) are significantly higher than previously published values. Inferred entrapment depths from our dataset range from 5.7 to 17 km – far greater than previous estimates – and support a vertically elongated magmatic system where crystallisation spanned both upper- and mid-crustal depths. Our CO₂ measurements enable new estimation of CO₂ sources and fluxes. As a total of 4.5 Mt of CO₂ was held in SHV magma during the aforementioned phases, the maximum amount of CO₂ that can be emitted from a batch of SHV magma is ∼1500–1750 tonnes/day. Measured CO₂ fluxes are significantly higher, indicating additional input of CO₂ into the system from greater depths. Our study shows that including bubble-hosted CO₂ redefines understanding of the SHV plumbing system.

1 Introduction

Mitigation of risk at active island arc volcanoes is a significant challenge (Joseph et al., 2022) that requires detailed understanding of the nature of the magmatic plumbing system and the origin and significance of its gas emissions. Soufrière Hills Volcano is arguably one of the most monitored and studied arc volcanoes, with detailed datasets and models constraining the geochemical, petrological, seismological and geodetic components of the volcanic system (e.g., Aspinall et al., 1998; Barclay et al., 1998; Devine et al., 2003; Ryan et al., 2010; Edmonds et al., 2014; Plail et al., 2018). The current eruption began in 1995, and has produced both intermittent explosions and dome growth between 1995 and 2010 (Wadge et al., 2014). While magma extrusion has not been observed since 2010, ongoing unrest includes 1) deformation of the volcano caused by melt injection into a crustal reservoir, proposed to be dual sourced at ∼5–6 km and ∼17 km depth (Neuberg et al., 2022) or mush based and vertically extensive, spanning the upper- and mid-crustal regions with a base at ∼17 km (Alshembari et al., 2024), 2) seismicity mainly in the form of volcano-tectonic earthquakes that are attributed to pressurisation and fracturing related to magma injection and migration (Smith, 2013), and 3) emission of SO₂ at an average rate of 374 ± 140 tonnes/day (from February 2010 to December 2014) with occasional higher fluxes (∼10×) accompanying volcano tectonic earthquake swarms and associated with migration of magma (Christopher et al., 2015).

During the ongoing hiatus of surface activity, some significant advances in technologies and methods are yet to be applied to the Soufrière Hills system. In particular, the application of Raman spectroscopy as a tool for measuring volatiles in vapour bubbles (e.g., Hartley et al., 2014; Moore et al., 2015; Moore and Bodnar, 2019) now permits more accurate measurement of total CO₂ in magmatic melt inclusions. Melt inclusions are small parcels of melt trapped during crystal growth and can represent the pre-eruptive magma, giving insight into the evolving geochemical and physical environment of subsurface magma within the plumbing system (e.g., Kent, 2008; Moore et al., 2015). The post-entrapment generation of bubble(s) allows sequestration of a proportion of the magmatic volatile species in the vapour phase–in particular, CO₂. Methodological development now allows the vapour bubble CO₂ to be quantified via Raman spectroscopy, revealing that in some cases more than 90% of the total CO₂ content of the inclusion resides in the bubble. This implies that previous measurements of melt inclusion CO₂, made solely on dissolved CO₂ within melt inclusion glass, could have severely underestimated magmatic CO₂ contents (Moore et al., 2015; Wieser et al., 2021). As a consequence, constraints on magma storage depths from CO₂-dependent volatile solubility models have been underestimated (e.g., Hartley et al., 2014; Moore et al., 2015; Hanyu et al., 2020). In this study, we use ion microprobe (SIMS) and Raman spectroscopy measurements to provide the first full CO₂ contents of melt inclusions from Soufrière Hills Volcano, using melt inclusions from four of the five phases of the 1995–2010 eruptive activity. This leads to refined and more realistic constraints on the storage depths and volatile systematics of this important eruption.

2 Geological setting

The Lesser Antilles Arc trends roughly N-S, and extends ∼850 km (Figure 1). The arc is the product of subduction of the Atlantic Plate at a 67° (ENE) vector beneath the Caribbean Plate, at an overall rate of ∼2 cm/year (DeMets et al., 2000), erupting 5 km³Ma^-1km^-1 of magma over a 300-year period from 1680 to 1980 (Wadge, 1984). The arc consists of two lines of volcanism, separated by age, marked by a central boundary at Martinique (Figure 1 inset). To the west exists the volcanic islands, where volcanism occurred during the Neogene to present (Macdonald et al., 2000), and to the east, islands where volcanism prevailed during the Upper Jurassic to the Lower Oligocene, with their volcanic basements now covered by sedimentary rocks (Bouysse et al., 1990). The arc can also be divided based on characteristics such as magma chemistry, seismicity, and overall structure (Balcone-Boissard et al., 2023; Metcalfe et al., 2023). The islands in the north, from Saba to Montserrat, can be tholeiitic (St Kitts, Redonda) and calc-alkaline (Saba, Montserrat), producing andesite (Rea, 1974; Brown et al., 1977; Rea and Baker, 1980; Baker, 1984; Macdonald et al., 2000). In the central and southern islands, andesites, basaltic andesites and basalts are most prevalent (Brown et al., 1977; Macdonald et al., 2000). Beneath the arc, the crust can be divided into four layers based on crustal structure, with the Mohorovičić discontinuity (Moho) varying between 25–37 km depth (Melekhova et al., 2019).

Figure 1

Figure 1. Map of Montserrat showing location of the four stratovolcanoes. Soufrière Hills is located in the south, denoted by the red star. Inset map illustrates the Lesser Antilles Arc and the western and eastern arcs.

2.1 Soufrière Hills Volcano, Montserrat

Montserrat is the sixth island from the north of the active volcanic chain and forms part of the northerly group of islands (Figure 1; Macdonald et al., 2000). It consists of four stratovolcanoes, three of which are dormant (Silver Hills, Centre Hills and South Soufrière Hills), last erupting 0.96 ± 0.25 million years ago (Rea, 1974). The active Soufrière Hills Volcano (SHV) is a volcanic complex located in the southern part of Montserrat, and its most recent eruption commenced on 18 July 1995, with the extrusion of crystal-rich andesitic magma (Robertson et al., 2000; Sparks and Young, 2002). This involved predominantly effusive activity, leading to the growth of lava dome complexes punctuated by dome-collapse events. Other activity involved explosive Vulcanian events of up to VEI 2–3 (Robertson et al., 2000; Cassidy et al., 2018), and less intense ash venting (Cole et al., 2014).

A total of 988 × 10⁶ km³ (Wadge et al., 2014) of material was produced over five phases of dome building activity that are interspersed with “pauses,” marked by a cessation in lava extrusion (Table 1). SHV has produced cycles of activity on both the sub-daily and sub-annual scale. Sub-daily cycles in Phase 1 consisted of explosions coinciding with ground deformation and seismicity (Voight et al., 1999), while in Phase 2, seismicity coincided with peaks in SO₂ flux (Young et al., 2003). This cyclic behavior ceased during Phase 3, and occurred again in Phases 4 and 5 in the form of seismic cycles where swarms of hybrid earthquakes merged to form continuous tremor (Cole et al., 2014; Odbert et al., 2014).

Table 1

Table 1. Characteristics of eruptive phases.

The 15-year eruptive activity is extensively detailed in a number of studies including Robertson et al. (2000); Kokelaar (2002); Sparks and Young (2002); Harford et al. (2003) for Phases 1-2 and Wadge et al. (2014) for Phases 15, and is summarised in Table 1 below.

2.1.1 Constraints on petrology and the plumbing system

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).

3 Methodology

3.1 Sample details and preparation

Samples from the five eruption phases were sourced from the Montserrat Volcano Observatory’s rock catalogue and are listed in Table 2. They were crushed and separated into different sized fractions, and plagioclase and orthopyroxene phenocrysts were hand separated from the 500–1000 µm fraction under a binocular microscope. For this study, plagioclase and orthopyroxene were selected due to their abundance and optical properties allowing for easy preparation and analysis of melt inclusions.

Table 2

Table 2. List of samples used in this study along with dates of production, brief descriptions and the analytical techniques applied.

Individual phenocrysts were mounted on glass slides and gently polished to a flat surface using 2400 silicon carbide lapping paper and inspected under a microscope. Samples containing glassy melt inclusions were further polished with 3 and 1 µm alumina paper until the inclusions were within ∼20–30 µm of the surface in preparation for Raman spectroscopy. Glassy melt inclusions occurred in Phases 1, 2, 4, and 5, with Phase 3 inclusions being crystallised, likely reflecting a longer cooling history, and are therefore not considered for the analytical techniques applied in this study.

Inclusions analysed by Raman spectroscopy were 10–166 µm along their longest axis and were cuboidal (plagioclase-hosted) or ellipsoidal (orthopyroxene-hosted). All plagioclase-hosted inclusions contained at least one vapour bubble ranging from 3 to 52 μm, and orthopyroxene-hosted inclusions were either bubble free or contained at least one bubble, where bubble size ranged from 2 to 23 µm.

Following Raman spectroscopy, inclusions >25 µm and therefore large enough to be analysed by SIMS were polished with 3 and 1 µm aluminium oxide paper, (to avoid carbon contamination posed by diamond paste) until the inclusion glass was exposed at the surface. These samples were then washed in acetone to dissolve any adhesive resin, mounted in indium, and gold coated for SIMS analysis.

Following SIMS, samples were lightly polished with 0.25 µm diamond paste to remove the gold coat and carbon coated for EPMA.

Whole rock samples from all five phases were also crushed and separated into the 100 µm fraction in preparation for ICP-OES.

3.2 Analytical techniques

3.2.1 Raman spectroscopy

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⁻¹.

3.2.2 Secondary ion man spectrometry (SIMS)

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₂.

3.2.3 Electron probe microanalysis (EPMA)

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.

3.2.4 Inductively coupled plasma optical emission spectrometry (ICP-OES)

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σ).

4 Results

4.1 Whole rock and host compositions

Whole rock compositions across the five phases are andesitic, spanning 57.84–59.97 wt% SiO₂ (Figure 2; Table 3) and contain the mineral assemblage outlined in Section 2.1.1. Samples from this study are comparable to published data across the five phases (Murphy et al., 2000; Zellmer et al., 2003; Plail et al., 2018), and are less evolved than groundmass glass compositions (70–80 wt% SiO₂) published for Phases 1–3 (Edmonds et al., 2001; Edmonds et al., 2002; Harford et al., 2003; Buckley et al., 2006; Humphreys et al., 2010).

Figure 2

Figure 2. Total alkali vs. silica (TAS) plot of whole rock samples and melt inclusions from SHV. Whole rock are predominantly andesite composition except for two samples. Melt inclusions are dacitic to rhyolitic. Published data for whole rock are from Murphy et al., 2000, Zellmer et al., 2003; Plail et al., 2018. Published data for melt inclusions are from Barclay et al., 1998, Devine et al., 1998, Edmonds et al., 2001; Humphreys et al., 2010.

Table 3

Table 3. Whole rock compositions of samples from Phases 1-5 in wt%. Total iron is given as Fe₂O₃ and H₂O is loss on ignition.

Fifty-eight of the melt inclusions across the four phases being studied are hosted in An_48-58 plagioclase, with three being hosted at An_61-62 and two at An_68-71. Orthopyroxene data are only available for Phases 1 and 5 due to a restricted number of melt inclusions that are sufficiently large, glassy and crystal-free. Orthopyroxene phenocrysts from both phases occupy restricted compositional range of En_57-59.

All plagioclase phenocrysts across the studied Phases are out of equilibrium with the whole rock, at a total K_D range of 0.16–0.54 (applicable K_D range for equilibrium = 0.05–0.15; section 4.2.1). However, equilibrium is achieved between 11 of 16 Phase 1 inclusions and the average groundmass glass composition for Phase 1 (Edmonds et al., 2002; Harford et al., 2003; Buckley et al., 2006). None of the fourteen Phase 2 inclusions are in equilibrium with the average groundmass glass composition for Phase 2 (Edmonds et al., 2002; Buckley et al., 2006; Humphreys et al., 2010). Equilibrium is unable to be calculated between inclusions and the average groundmass compositions for Phases 4 and 5 as there are no published values of Phase 4 and 5 groundmass glasses.

The five analysed orthopyroxene melt inclusions from Phase 1 are out of equilibrium with the average whole rock (K_D = 0.49–0.55; section 4.2.1), along with the groundmass glass (Edmonds et al., 2001; Harford et al., 2003; Buckley et al., 2006) at a K_D range of 0.18–0.19. Equilibrium is also not achieved between the nine Phase 5 inclusions and the average whole rock composition (K_D = 0.49–0.55), and equilibrium between mineral and groundmass glass for Phase 5 cannot be calculated due to lack of measured groundmass glass compositions for Phase 5.

4.2 Melt inclusions

4.2.1 Post-entrapment modifications

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.

4.2.1.1 Post entrapment crystallisation (PEC)

Assessment of PEC for plagioclase-liquid pairs was carried out by two methods. Firstly, we considered the anorthite-albite exchange between the inclusion and its host which varies with temperature, where $K_{D} = 0.10 \pm 0.05$ for inclusions trapped <1050°C and $K_{D} = 0.28 \pm 0.11$ for those trapped at ≥1050°C (Putirka, 2008). As magmatic temperature at SHV are ∼850°C based on petrological and experimental studies (Sparks and Young, 2002), the lower temperature K_D value of 0.10 ± 0.05 was used to assess equilibrium. According to this equilibrium test, 52% of inclusions were in equilibrium with their hosts, with a total K_D range of 0.02–0.07. However, due to the hydrous nature of the inclusions, K_D may not be an accurate indicator of equilibrium for SHV inclusions, as equilibrium is affected by H₂O degassing (Humphreys et al., 2016), and a second method was employed. As an alternative test for equilibrium, the Al₂O₃-SiO₂, Al₂O₃-MgO, Al₂O₃-K₂O and MgO-K₂O systematics of the melt inclusions in relation to the established liquid line of descent (LLD) for SHV whole rock and groundmass glasses were assessed, as PEC leads to a decrease in Al₂O₃ and an increase in MgO (e.g., Figure 3; Nielsen, 2011). As no deviation from any of the tested LLDs occurred, this is interpreted as an indication of no PEC occurring in the plagioclase-hosted inclusions, and therefore no correction being required.

Figure 3

Figure 3. Various liquid lines of descent for whole rock and groundmass glass compositions, in relation to melt inclusion compositions. Red lines indicate the amount (∼15%) and direction that compositions of plagioclase-hosted inclusions would follow as a result of post-entrapment crystallization. Therefore, SHV melt inclusions do not show any indication of post-entrapment crystallisation, as data do not follow the red lines. (A) Al₂O₃ vs. K₂O. (B) MgO vs. K₂O. (C) Al₂O₃ vs. MgO. (D) Al₂O₃ vs. SiO₂.

Orthopyroxene-liquid pairs were tested for equilibrium according to the K_D threshold of ${K_{D}}_{p x - l i q}^{F e - M g} = 0.29 \pm 0.06$ (Putirka, 2008). All inclusions were out of equilibrium with their hosts at a K_D range of 0.05–0.22, and PEC ranges from 0.96% to 3.33%. Due to the low amount of PEC, the compositions of orthopyroxene-hosted melt inclusions do not require correcting, as this process has been shown to have negligible effects on both major and volatile elements up to 11% PEC (Moretti et al., 2018).

Melt inclusion compositions are listed in Table 4.

Table 4

Table 4. Compositions for seventy-eight melt inclusions across Phases 1 - 5 for Soufrière Hills Volcano. Major element oxides in wt% are measured via EPMA. H₂O and CO₂ in the glass are measured by SIMS. CO₂ in the bubble is measured via Raman spectroscopy. PEC refers to post-entrapment crystallisation as assessed in Section 4.2.1.

4.2.1.2 Bubble growth

After entrapment of melt, bubbles can be formed in response to the pressure-volume-temperature relationship between host mineral and melt. Based on this relationship, bubbles can grow via 1) post-entrapment crystallization, 2) diffusive H⁺ loss and 3) differential thermal contraction (Roedder, 1979; Anderson and Brown, 1993; Lowenstern, 1995; Kent, 2008; Hartley et al., 2014; Wallace et al., 2015; Aster et al., 2016; Hanyu et al., 2020; Wieser et al., 2021). The process with the largest effect occurs due to differing thermal expansivities between host mineral and melt during cooling from high trapping temperatures to the glass transition temperature, after which bubble growth is suppressed (Moore et al., 2015). However, bubbles are also known to nucleate in the melt prior to entrapment, and can grow by diffusion of volatiles from the inclusion glass into the bubble, coalescence, decompression during ascent or Ostwald ripening (Cashman and Mangan, 1994; Best, 2013). It is therefore necessary to identify bubbles that have grown homogenously post-entrapment, in order to prevent overestimation of the bubble volatile content caused by volatile-bearing bubbles being trapped at the time of melt inclusion formation. Homogenous bubble growth has been reported to be 5%–12% for a range of volcanic systems (e.g., Hartley et al., 2014; Aster et al., 2016; Hanyu et al., 2020), and we adopt the lower end of 5% bubble volume to distinguish bubbles that nucleated and grew post-entrapment, in comparison to those trapped with the melt.

Previous studies combining CO₂ analysis in vapour bubbles and their host melt inclusions have shown that up to 90% of CO₂ can be sequestered to bubbles in the form of vapour (Hartley et al., 2014; Moore et al., 2015; Venugopal et al., 2020). However, CO₂ can exist in its liquid form and also as carbonates that can store up to 50% of CO₂ in the bubble. This is also true for sulphur-bearing minerals precipitated on bubble walls, which can store up to 60% of sulphur originally trapped in an inclusion (Esposito et al., 2016; Schiavi et al., 2020), and H₂O of which up to 16% can be sequestered (Esposito et al., 2016). Overall, melt inclusion bubbles have the ability to not only store large amount of CO₂ and S, but also H₂O and major and minor elements that constitute carbonates, sulphates, sulphides, halides and other minerals such as Na, Ca, Mg, Fe or Cu (Schiavi et al., 2020; Venugopal et al., 2020).

The concentration of CO₂ sequestered to bubbles post-entrapment is calculated by mass-balance equations (Moore and Bodnar, 2019), taking into consideration: 1) the volume fraction of the melt inclusion that is constituted by the bubble, and 2) the density of CO₂ measured by Raman spectroscopy (Supplementary Data Sheet). Total inclusion and bubble volumes were estimated from photomicrographs, assuming a spherical shape for bubbles, an ellipsoid shape for orthopyroxene hosted inclusions and cuboidal shapes for inclusions hosted in plagioclase, based on their 2D appearances from a polished surface (Figure 4). The two perpendicular axes were measured using a Leica DM4500 P LED microscope on the Leica Application Suite software, and the third unseen axis was calculated using the arithmetic mean of the measured axes. This method is associated with an average 5% error, but a 1σ error of −48% to 37% (Tucker et al., 2019). Uncertainty in bubble sizes was ±2 µm,w which yielded bubble volume uncertainties of 6%–24%. (Mean 11%; Supplementary Data Sheet). Density calculations were undertaken by firstly processing the Raman spectra of individual bubbles using the WiRE™ 4.4 spectral analysis software. After baselines were applied to each spectrum with >500 counts using a polynomial curve, the Fermi diad was truncated at 1200 and 1500 cm⁻¹. Each peak was fitted with a mixed Gaussian-Lorentzian curve, and the Fermi diad split was calculated as the difference between the centres (in wave number) of the two peaks. In the absence of a CO₂ densimetry curve specifically calibrated for the instrument used, the experimental equation of Lamadrid et al. (2017) was adopted to calculate CO₂ density. Fermi diad peaks with counts <500 or asymmetrical peaks to which curves could not be readily fitted were not used to quantify CO₂.

Figure 4

Figure 4. Examples of melt inclusions adopting different shapes. (A) SHV_P1_PLAG_010 is a plagioclase hosted inclusion taking the 2D shape of a rectangle, assumed to be cuboidal in 3D. (B) SHV_P5_OPX_006_MI1 is an orthopyroxene hosted inclusion taking the 2D shape of an oval, assumed to be ellipsoidal in 3D.

4.2.1.3 H₂O loss

Loss of H₂O in melt inclusions was assessed from H₂O vs. K₂O, which shows no significant variation of H₂O at a given K₂O (Figure 7), and thereby indicates no significant H₂O loss in SHV melt inclusions.

4.2.2 Major element compositions

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.

Figure 5

Figure 5. Harker diagram showing major elements plotted against K₂O for assessment of variation with magma differentiation. The composition of inclusions differs with phase, where there is an increase in evolution through Phases 1 to 4, followed by less evolved inclusions at Phase 5. Grey zones are data reported in the literature for Phases 1-3 from Edmonds et al. (2002), Buckley et al. (2006), Humphreys et al. (2010), Mann et al. (2013). (A) SiO₂ vs K₂O, (B) MgO vs K₂O, (C) TiO₂ vs K₂O, (D) CaO vs K₂O, (E) Al₂O₃ vs K₂O, (F) Na₂O vs K₂O, (G) FeO 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).

4.3 Volatile element compositions

4.3.1 Melt inclusion glass

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.

4.3.2 Melt inclusion bubble

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.

Figure 6

Figure 6. Raman spectra of gaseous and solid phases found in Soufrière Hills bubbles. (A) Raman spectra of bubble from SHV_P2_OPX_001_MI1 displaying Fermi diad indicating the presence of CO₂. (B) Raman spectra of the bubble of plagioclase-hosted inclusion ‘SHV_P1_PLAG_030_MI1’ showing Fermi diad. (C) Raman spectra of bubble in SHV_P5_OPX_002_MI1 showing presence of sulphate (anhydrite) crystals on bubble wall with peaks at 430, 500, 611, 629, 676, 1018, and 1131 cm⁻¹.

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} = \frac{ρ_{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.

4.3.3 Total CO₂

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).

Figure 7

Figure 7. (A) H₂O vs. K₂O in the glass shows degassing within Phases and an overall decrease of H2O with time, reverting to high H2O during the last Phase of eruption. (B) Total CO₂ vs. H₂O does not follow simple open or closed degassing trends. (C) Concentration of CO₂ in the glass, bubble and total CO₂ of the MI. Grey areas are values reported in the literature. H₂O and CO₂ data from Barclay et al. (1998); Mann et al. (2013) & Edmonds et al. (2014). Error bars show 10% errors on H₂O, and are calculated individually for CO₂ based on the contribution from both the bubble and the glass. Where error bars are not shown, the error is smaller than the symbol size.

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.

Figure 8

Figure 8. (A) Percentage of CO₂ contributed by the bubble vs. bubble to inclusion ratio. In orthopyroxene-hosted inclusions, 5%–76% of CO₂ in the total melt inclusion is contributed by the bubble. In plagioclase-hosted inclusions, this value is 37%–99%. Overall, there is a median of 90% contribution from the bubble to the total inclusion CO₂ (B) There is no correlation with CO₂ density as a function of bubble volume percentage. The red line on both plots denote the 5% bubble threshold of homogenous vs. heterogeneous bubble growth.

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).

5 Discussion

5.1 Magma storage conditions

The compositions (corrected for PEC) of major and volatile elements are applied to models to facilitate connection between melt inclusion data and crystallisation temperatures, magma storage depths, degassing systematics and magmatic CO₂ budgets.

5.1.1 Geothermometry

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).

Figure 9

Figure 9. Pressurevs temperature graph for inclusions at SHV showing the magmatic conditions at the time of plagioclase and orthopyroxene crystal growth. Grey horizontal dashes represent the depths of magma storage based on deformation modelling (Foroozan et al., 2010). Error bars represent the uncertainty in depths based on uncertainties in total CO₂. Error bars for orthopyroxene-hosted melt inclusions are not shown as the uncertainty is smaller than the size of the symbol.

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).

5.1.2 Saturated pressures and associated depths

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} = \frac{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).

Figure 10

Figure 10. Schematic of the plumbing system beneath Soufriere Hills Volcano, exhibiting mush properties similar to other arc volcanoes. This is a vertically extensive transcrustal magmatic system consists of a melt-rich zone with high H₂O and CO₂ contents and a crystal rich zone dominated by calcic plagioclase and orthopyroxene spanning ∼5–17 km, and is periodically flushed with CO₂ from deeper in the system. Melt inclusions from Phases 1 and 5 plot in the mid crustal region whereas Phases 2 and 4 span the upper to mid crustal regions.

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.

5.2 CO₂ budget

The maximum measured CO₂ concentration in melt inclusions glass and bubble, representing the least degassed CO₂ trapped in melt inclusions is compared to estimates of CO₂ degassed through the plume using the petrological method (Devine, et al., 1984), assuming that the concentration of volatiles in melt inclusions represents the concentration of volatiles dissolved in the magma at depth. A total of 4.5 Mt of CO₂ (M_dv) is calculated to have been contained in SHV magma during Phases 1, 2, 4, and 5 based on Equation 3:

M_{d v} = X_{M I} V ρ (3)

where X_MI is the maximum concentration of CO₂ measured in melt inclusions (3230 ppm) in the studied Phases, V is the volume of erupted magma (Wadge et al., 2014), and ρ refers to the magma density of 2400 kg m⁻³ Dense Rock Equivalent (Melnik and Sparks, 2002). Phases 1 and 2 held more dissolved CO₂ in their magmas with maxima of 1.2 and 2.6 Mt respectively, exceeding the corresponding figures for Phases 4 and 5 (0.25 and 0.44 Mt respectively). This total of 4.5 Mt of CO₂ dissolved in the magma is however much less than the 7.7–12 Mt total plume-emitted CO₂ estimated by Edmonds et al. (2014) during the period 1995–2009, using an average molar CO₂/SO₂ ratio of 5.1. However, the ratio of 5.1 used in their estimation introduces a large error, as it is derived from 5 days of measurements during Pause 3, and assumes a constant ratio over 14 years of activity. This mismatch of petrological estimates in comparison to direct measurements is often seen in SO₂ fluxes, and is a feature common to arc settings where magmas are more oxidized, and is not generally observed in non-arc magmas (e.g., Sigurdsson, et al., 1990; Westrich and Gerlach, 1992; Gerlach et al., 1994; Blake, 2003; Sharma et al., 2004).

Applying the 5 days of CO₂ measurements (Edmonds et al., 2014) in the plume taken during July 2008 to a shorter time period of Phase 4 (July 2008-January 2009), we calculate that 0.24 Mt was released in the plume over 158 days. Therefore 96% of the CO₂ dissolved in Phase 4 magma was emitted during the eruption. Applying the calculations to Phase 5 (October 2009-January 2010), an average of 2297 tonnes/day of CO₂ emitted from June-November 2010 (Edmonds et al., 2014) amounts to 0.29 Mt or 66% of the CO₂ dissolved in the magma being released over 125 days. While these estimates are termed ‘maxima’ based on the maximum CO₂ considered for each Phase, it is noted that masses calculated can be underestimated due to 1) unerupted volumes of magma which are not taken into account, and 2) entrapment of already degassed melt (Moune et al., 2007). The percentage of CO₂ released during Phases 1 and 2 is unable to be estimated due to lack of published CO₂ fluxes.

Original magmatic CO₂ contents can also be estimated using the CO₂ flux derived by combining SO₂ fluxes measured using UV spectrometer and CO₂/SO₂ ratio acquired via Multisensor Gas Analyser System (MultiGAS) over the 5 days in July 2008 and 15 days during June to October 2009. The highest CO₂ flux reported (Edmonds et al., 2014) of 5494 tonnes/day measured on 1 October 2010 is used to calculate a maximum mass of CO₂ for a period of 125 days during Phase 5, amounting to 0.69 Mt, and magma masses are calculated from reported magma volume estimates and magmatic density of 2400 kgm⁻³. We can therefore use Equation 3 to estimate the required melt inclusion CO₂ concentration to match the measured CO₂ emissions. This value of 4107 ppm is 877 higher than the initial CO₂ from the least degassed melt inclusions measured across the eruption. This value correlates to a pressure of 8090 bars or 31.0 km depth assuming a H₂O content of 6.39 wt%, the highest of the Phase 5 inclusions. As a result, we attribute this high CO₂ flux to CO₂ flushing (e.g., Blundy et al., 2010; Caricchi et al., 2018) from degassing deeper in the system, a process that is commonplace in the arc, and occurs in all island groups (Metcalfe et al., 2023). The related depth would extend into the Moho situated ∼30 km under Montserrat, based on crustal xenolith petrology and seismic data (Melekhova et al., 2019). Based on fluxes generated during June to October, and compared to Phase 5, a CO₂ flux greater than ∼1500–1750 tonnes/day cannot be accounted for by the dissolved CO₂ in this batch of magma, and therefore relate to CO₂ flushing from mantle depths.

5.3 CO₂ controls on explosivity and implications for the future

Quantification of dissolved magmatic gases as well as their emissions is essential for understanding the influence of volatiles on volcanic system dynamics and the related hazards. The dissolved magmatic CO₂ concentration can, alongside more abundant H₂O, drive eruptions, dictate eruptive style or impact both in volcanic settings (e.g., Burton et al., 2007; Allard, 2010; Allison et al., 2021). The eruptive style at SHV can be described as transitioning, with both effusive and explosive eruptions occurring during each Phase, producing lava domes and explosions (e.g., Druitt et al., 2002; Edmonds and Herd, 2007; Ryan et al., 2010). However, at SHV, microlite petrology from Phase 3 samples has demonstrated that shallow conduit processes determine eruptive style, with similar decompression pathways within the deeper magmatic system resulting in both lava dome growth and explosions (Murch and Cole, 2019). Nevertheless, melt inclusion analysis conducted in this study shows that while each Phase produced VEI 2-3 events, their SiO₂ content, H₂O content, CO₂ content, duration, and the number of explosions differ. Particularly, Phase 1 and 5 are similar in silica content, H₂O and also have similar pressures and depths based on solubility modeling of melt inclusion H₂O and CO₂ concentrations. While these Phases exhibit similar chemistry, their explosive behaviours have differed, with Phase 5 producing more Vulcanian explosions within a shorter timeframe (125 days) than Phase 1 (846 days). This is likely due to the amount of CO₂ in the system during these Phases as Phase 1 had ∼1,000 ppm less CO₂ dissolved within its melt and Phase 5 was influenced by CO₂ flushing. However, due to the limited availability of CO₂ flux data at SHV, it is unclear if CO₂ flushing occurred during either Phases.

The monitoring of volatile release at SHV has primarily focused on SO₂ flux (e.g., Young et al., 1998; Carn and Prata, 2010; Nicholson et al., 2013; Christopher et al., 2015), however trends in SO₂ have not correlated with the eruptive Phases, but instead occur on multi-year (∼2 years) and multi-week cycles (∼50 days and ∼5 months), independent of lava extrusion and dome building (Nicholson et al., 2013; Christopher et al., 2015), upon which Phases are determined. We recommend that CO₂ should also be routinely monitored as it provides insight into movement of magma in the deep system, and also has the ability to drive eruptions through processes such as CO₂ flushing, which can alter the physical properties of the magma - increasing magma volume, leading to overpressure and potentially eruptions (Caricchi et al., 2018). While data are severely limited, instances of CO₂ flushing at SHV are evident at the transition between Pause 3 and Phase 4, and ∼6 months after the cessation of Phase 5. As the crystallinity of SHV magma, at 30%–45%, is at or beneath the threshold at which magmas reach their rheological locking point (40%–50%), flushing at SHV can potentially trigger eruptions. This further highlights the importance of monitoring of CO₂ flux at SHV, especially during the current period of unrest.

6 Conclusion

We have provided the most comprehensive study published to date of magmatic carbon in the Soufrière Hills Volcano magmatic system, providing the first measurements to include bubble-hosted CO₂ in melt inclusions and the first pre-eruptive CO₂ measurements for the most recent eruptive phases. H₂O and total CO₂ contents for melt inclusions hosted in plagioclase from pumice erupted during four phases of eruption are 2.47–6.40 wt% H₂O and 131–3230 ppm CO₂. Those hosted in orthopyroxene are 5.15–7.26 wt% and 199–736 ppm CO₂. The sequestration of CO₂ to the bubble has proven to be an important post-entrapment modification at SHV, with 5%–99% (median 90%) of CO₂ contained in bubbles and thus excluded in earlier analyses. Our measurements therefore provide unprecedented insight into the magmatic CO₂ concentrations within the SHV system and indicate higher concentrations than previous studies. Our volatile solubility-based estimates of magma storage depths (5.7–17 km) indicate that a transcrustal magmatic system was persistent throughout the eruption, consistent with geodetic modelling (Alshembari et al., 2024).

We use melt inclusion compositions to estimate magmatic temperatures at entrapment, based on liquid-mineral pairs. Temperature results (828°C–916°C) are similar to published estimates for Phase 1 (812°C–891°C) and are the first temperature estimates for Phases 2, 4, and 5.

The total amount of CO₂ contained in magma from Phases 1, 2, 4, and 5 is ∼4.5 Mt, using the petrological method. We calculate that 96% of the dissolved CO₂ in Phase 4 magma was released to the atmosphere during Phase 4, and 66% of the dissolved CO₂ in Phase 5 was released. Our results indicate that ∼1500–1750 tonnes/day CO₂ can be emitted from the magma, and thus any CO₂ fluxes higher than 1750 tonnes/day are likely the result of CO₂ flushing from deeper in the magmatic system. Pauses 4 and 5 therefore experienced CO₂ flushing, with up to 5494 tonnes/day being released (Edmonds et al., 2014), and are the only Phases with reported CO₂ fluxes.

Our new melt inclusions dataset covers four of five Phases of eruptive activity at SHV, and provides the first account of total (glass + bubble hosted) CO₂. Our data indicates that the magmatic system feeding SHV is a vertically elongated mush system, as corroborated by geophysics and geochemistry.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

TH: conceptualization, writing–original draft, writing–review and editing, investigation, methodology, formal analysis, visualization, funding acquisition. TC: conceptualization, writing–review and editing, formal analysis, supervision. SM: conceptualization, writing–review and editing, formal analysis, supervision. HT: conceptualization, writing–review and editing, formal analysis, supervision, funding acquisition. FS: investigation, methodology, formal analysis, writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This study was supported by a Government of Montserrat PhD scholarship awarded to TH. Funding for ionprobe analysis was provided by NERC IMF grant IMF746/1122 awarded to HT and TH. Raman and whole rock analyses were graciously covered by the Volcanology team at Laboratoire Magmas et Volcans. Article Processing Fees were provided by Lancaster University.

Acknowledgments

We thank the Montserrat Volcano Observatory for granting use of samples. Many thanks to Cees-Jan de Hoog for assistance with SIMS, Iris Buisman for EPMA assistance and Claire Fonquernie for whole rock analysis. We are grateful to David Neave and Margaret Hartley for discussions on melt inclusions and in particular their reconstruction.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2024.1509409/full#supplementary-material

References

Allard, P. (2010). A CO₂-rich gas trigger of explosive paroxysms at Stromboli basaltic volcano, Italy. J. Volcanol. Geotherm. Res. 189 (3-4), 363–374. doi:10.1016/j.jvolgeores.2009.11.018