Petrology of the 2020–21 effusive to explosive eruption of La Soufrière Volcano, St Vincent: insights into plumbing system architecture and magma assembly mechanism Open Access

After more than 40 years of dormancy, La Soufrière volcano on the island of St Vincent, Eastern Caribbean, started erupting on 27 December 2020 (Fig. 1). The eruption commenced with extrusion of a highly viscous lava dome onto the floor of the main crater, adjacent to the lava dome that had been emplaced in 1979. Over the next three months, lava formed an elongated coulee within the crater with a total volume of c. 18 × 10⁶ m³ (Joseph et al. 2022; Stinton et al. 2023). The eruption entered a new phase of activity on 9 April 2021, when it produced a series of explosive eruptions that generated both tephra-fall and pyroclastic density currents (Cole et al. 2023). The timely evacuation of about 16 000 people (Joseph et al. 2022) was a great success, especially in the context of the 1902/03 eruption of La Soufrière St Vincent that caused the tragic death of more than 1500 people (Pyle et al. 2018).

Renewed activity of La Soufrière provides the opportunity to use petrological indicators to investigate the structure and dynamics of the subvolcanic magma feeding system, its controls on the observed activity, and its relation to previous activity. Mobilization processes of magma from storage zones in the Earth's crust have been reconstructed at numerous arc volcanoes. Many studies have used mineral chemistry and textures to link the input of fresh magma pulses into pre-existing reservoirs (i.e. crystal mushes) to eruption initiation (e.g. Tepley et al. 2000; Humphreys et al. 2006; Ruprecht and Wörner 2007; Andrews et al. 2008; Kent et al. 2010; Druitt et al. 2012; Bouvet de Maisonneuve et al. 2015; Wanke et al. 2019; Weber et al. 2020). Yet, examples without clear evidence for a relationship between magma recharge and eruption triggering are not uncommon; in such cases, exsolution of dissolved volatiles is frequently invoked in the mobilization of resident magmas (e.g. Blake 1984; Kent et al. 2007; Cassidy et al. 2016; Stock et al. 2016; Budd et al. 2017; Andersen et al. 2018; Caricchi et al. 2018). External forcing of magmatic systems, for example by large earthquakes, may also contribute towards destabilization of magma feeding systems and initiation of volcanic eruptions (e.g. Walter and Amelung 2007; Bebbington and Marzocchi 2011; Nishimura 2017; Watt 2019; Seropian et al. 2021). Moreover, while some volcanoes show persistent evidence for specific reactivation processes such as mafic magma recharge (e.g. Nevado de Toluca: Weber et al. 2020; Mount Hood: Kent et al. 2010), it is common for both the depth from which magma ascends and the types of mobilization processes to vary. This can be observed even during different phases of the same eruption or different events at the same volcano (e.g. changes in reservoir depth at Parinacota: Ginibre and Wörner 2007; vapour phase accumulation v. silicic magma recharge at Mt St Helens: Kent et al. 2007; Cashman and Blundy 2013). It is therefore important to reconstruct magma storage and mobilization conditions for individual eruptions to link specific petrological processes to geophysical monitoring signals, which can in turn provide valuable insights for future volcanic hazard mitigation.

Based on analysis of earlier eruptive products, current petrological models of the magmatic plumbing system beneath La Soufrière consider it to be characterized by a vertically extensive and compositionally diverse plexus of storage zones (Fedele et al. 2021). The recent eruptions of La Soufrière, St Vincent provide an excellent opportunity to compare mineralogically derived pressure estimates with observed depths of unrest and activity recorded by contemporaneous geophysical detection and related surface activity. Pressure and temperature reconstructions employing mineral–melt equilibria are challenging at La Soufrière due to the high crystallinity and groundmass microlite content of the erupted products. However, such estimates are required to constrain the depths from which magma was sourced, to understand the reactivation process that primed the feeding system for eruption and to explain why La Soufrière eruptions often switch between explosive and effusive phases (Robertson 1995).

Here we present detailed petrography and mineral chemistry from the now destroyed 2020–21 lava dome and subsequent explosive phase of the 2020–21 eruption of La Soufrière St Vincent. We use a petrological modelling approach based on melt inclusion chemistry to reconstruct a wide range of magma system variables. Based on these data, we propose a conceptual model of the magma feeding system and mobilization processes for the 2020–21 eruption that we compare with the geophysical indicators of melt movement and degassing obtained during the course of the eruption. This paper complements that of Frey et al. (2023), which focuses on the evolving compositions and textures of microlites and their relative abundances during the effusive and explosive sequence.

La Soufrière, St Vincent (13° 20′ 20″ N, 61° 10′ 43″ W) is one of the most active and hazardous stratovolcanoes in the Eastern Caribbean (Robertson 1995). Located in the southern part of the Eastern Caribbean Volcanic Arc, it is part of an 800 km-long chain of subaerial and submarine volcanoes fuelled by hydrous basaltic melt produced in response to westward subduction of the North and South American plates beneath the Caribbean plate (Macdonald et al. 2000; Cooper et al. 2020). Throughout its history, La Soufrière has produced basalts and basaltic andesites of monotonous transitional calc-alkaline to tholeiitic affinity (Heath et al. 1998a). While the eruptive products are restricted in geochemistry, the volcano exhibits a range of eruptive styles that includes lava flows and domes as well as explosive paroxysms that often feature pyroclastic density currents and lahars (Robertson 1995; Cole et al. 2019). In historical time, La Soufrière has produced at least four explosive eruptions (Cole et al. 2019), including a large paroxysm in 1902/1903 classified as a Volcanic Explosivity Index (VEI) of 4 (Pyle et al. 2018). Between 1903 and the 2020–21 eruption, the volcano had two eruptive episodes: a purely effusive event in 1971–72 that produced a lava dome in the water-filled central crater (Aspinall et al. 1973) and a series of explosions in 1979 followed by the extrusion of a large lava dome (Shepherd et al. 1979; Graham and Thirlwall 1981; Fiske and Sigurdsson 1982).

Since the 1980s, petrological, experimental, geochemical and isotopic studies have been used to develop conceptual models of the structure and dynamics of the igneous plumbing system beneath St Vincent. The generation of primary basaltic melts, in particular, has received much attention through petrological experiments (Pichavant et al. 2002; Pichavant and Macdonald 2007; Melekhova et al. 2015, 2019), whole-rock and mineral chemistry (Heath et al. 1998a), and melt inclusion studies (Bouvier et al. 2008, 2010). High-MgO basalts from St Vincent are found to represent primary mantle melts generated by partial melting of a MORB-like lherzolitic mantle source previously metasomatized by slab-derived fluids (McDade et al. 2003). The last equilibration conditions of these basalts with their mantle source lie between 11.5 and 16 kbar and 1185 to 1235°C under relatively oxidizing conditions above the quartz–fayalite–magnetite buffer (QFM + 1–QFM + 2 log units) and with water contents of 1.5 to 4.5 wt% (Heath et al. 1998b; Pichavant et al. 2002; Pichavant and Macdonald 2007). Using primitive melt inclusions containing 1.5 to 4 wt% H₂O, Bouvier et al. (2008, 2010) argued for a slightly more restricted (but consistent) range of mantle melting conditions, 13–14.5 kbar and 1220–1190°C.

Basaltic andesites from La Soufrière show mantle-like isotopic signatures, consistent with an interpretation that crustal melting does not play a major role in the petrogenesis of rocks at St Vincent (Graham and Thirlwall 1981; Heath et al. 1998a; Tollan et al. 2012). Eruptive products are typically highly porphyritic; phase equilibria considerations suggest that basaltic andesites are derived from high-MgO basalt parents through fractional crystallization at pressures of 7 to 10 kbar and temperatures between 1050 and 1130°C (Heath et al. 1998a; Melekhova et al. 2015). Trace element data show, however, that basaltic andesites from recent eruptions (1971, 1979) may not represent simple derivatives of each other, but instead derive from separate magma batches (Graham and Thirlwall 1981). Further inferences on magma differentiation processes have been drawn from detailed studies of coarse-grained ‘plutonic’ gabbroic to ultramafic xenoliths, which are ubiquitous in La Soufrière magmas (Wager 1962; Arculus and Wills 1980; Tollan et al. 2012; Melekhova et al. 2015, 2019).

Insight into eruptive processes comes from crustal residence times of basaltic andesite magma batches of tens of thousands of years as determined from U–Th disequilibria (Heath et al. 1998b). Comparison of published experimental data and SEM-EDS mineral compositions from the historical and prehistoric (pre-AD 1700) activity of La Soufrière suggests that eruptions are initiated at shallow crustal levels (c. 1–4 kbar) by volatile accumulation pressure but may also involve episodic rejuvenation of shallow crustal basaltic andesite intrusions with more mafic magma (Fedele et al. 2021). Diffusion modelling of Sr concentrations in plagioclase crystals indicates shallow magma residence times of about 150 years (Zellmer et al. 1999).

In summary, previous constraints based on petrological experiments on the subvolcanic plumbing system of St Vincent indicate deep generation of primitive basaltic magmas in the mantle at 11–16 kbar and 1185 to 1235°C. These parental magmas ascend and fractionate to produce basaltic andesite at c. 7–10 kbar and temperatures of 1050–1130°C, before being stored and erupted from shallow (1–4 kbar) crustal levels. Plutonic xenoliths capture snapshots of this process throughout the crust (Melekhova et al. 2019). In this study we provide the first detailed mineral–melt thermobarometry for La Soufrière, St Vincent to better constrain the magma source depth and evaluate eruption initiation processes and timescales prior to the 2020–21 eruption.

The samples analysed in this study were collected during two field campaigns in January and April 2021, targeting the effusive and explosive phases of the eruption respectively. A sample of the 2021 lava dome collected on an active lobe of the advancing flow on 16 January 2021 (13° 20′ 1.71″ N, 61° 11′ 10.36″ W) was split into 8 subsamples. This represented a discrete subsample of one part of the dome at a relatively early stage of its growth. Importantly, because of safety and logistical concerns, spatial and temporal heterogeneity could not be assessed. Juvenile scoria samples were collected from the upper parts of the April pyroclastic deposit, where clasts were largest. These clasts represent scoria generated during the later explosions of the sequence and were sampled in three separate locations in the east and NE of the island. We analyse 14 samples of the explosive phase that were very likely erupted on 10/11 April 2021. These samples are from clasts sampled in April and May 2021 at five locations around the volcano, representing later phases of the explosive eruptions (U3 to U7 Cole et al. 2023). The PDC clast samples were quenched on sampling so assumed to represent fragmented magma incorporated into the PDC. The sampled materials and locations, GPS coordinates and performed analyses are listed in Supplementary Table S1.

Five petrographic thin sections of the 2021 lava dome sample, three of the late-stage scoria (U3 and U5 unit; Supplementary Table S1), were examined using polarized light microscopy. To establish the proportions of individual phases, we obtained full photographic scans of all thin sections in plane and crossed polarized light. Modal percentages were determined on the scans by random grid-point counting using the software package JMicroVision version 1.3.1 (Roduit 2007). To establish robust modal percentages, a grid of random locations (points) was allocated on each thin section. For each random point, the phase was identified and counted until the phase proportions converged to constant values, which was readily achieved with more than 500 counted points. All modal percentages were calculated on a vesicle free basis, which required high-resolution back-scatter electron (BSE) imaging to resolve scoria micro-vesicularities. BSE images were collected using a JEOL JXA-8200 electron microprobe at the Research Laboratory for Archaeology and History of Art, University of Oxford.

The major and minor element compositions of selected bulk-rock samples were determined by x-ray fluorescence analysis (XRF) at the University of Plymouth (UK) and by inductively coupled plasma optical emission spectrometry (ICP-OES) at AcmeLabs, a division within Bureau Veritas (US). For the whole-rock major element chemical analyses by ICP-OES, powdered samples were prepared with a lithium borate flux followed by a highly aggressive dissolution. Detection limits of these analyses are <0.01% for most oxides. For XRF analysis, samples were mixed with lithium tetraborate and lithium metaborate flux and converted to glass beads. Measurements were carried out using a Panalytical Axios Max instrument.

The major and minor element composition of crystals was determined at the University of Oxford by electron probe microanalysis (EPMA), using a JEOL JXA-8200 at the Research Laboratory for Archaeology and History of Art and a Cameca SX-5 FE in the Department of Earth Sciences. Mineral analyses were obtained with a focused beam, 15 kV acceleration voltage and 15–20 nA beam current. Primary calibration was carried out on in-house standards: Albite (Na, Si, Al), Sanidine (K), TiO₂ (Ti), elemental Mn, Cr and Ni, and Andradite (Fe, Ca). Peak counting times for mineral analysis were as follows: 20 s for Na; 30 s for Mg, Si, Al, K, Fe, Ca; and 40 s Ti, Mn Cr, Ni. A subset of plagioclase analyses was carried out with increased beam current of 40 nA and peak counting times of 180 s to measure TiO₂ contents. Silicate mineral analyses that returned analytical totals outside the range 98.5 to 101.5% were discarded. At the beginning and end of each analytical session, we measured several secondary standards for quality control (Supplementary Table S10).

The bulk microlite-rich groundmass of scoria samples was analysed in 22 randomly distributed spots with a 20 µm defocused beam, 15 kV acceleration voltage and beam current of 10 nA. To mitigate diffusive alkali loss during analysis, Na and K measurements were placed first in the sequence with reduced peak counting times of 15 and 20 s respectively. Peak counting times for all other elements were 30 s. Analytical totals for interstitial groundmass analyses ranged between 98.6 to 99.9% indicative of low dissolved volatile contents.

Melt inclusions hosted in plagioclase and pyroxene crystals from the explosive phase of the eruption were analysed by EPMA using a JEOL JXA-8350 F in the School of Earth Sciences at the University of Bristol. Analysis conditions included a 5 nA beam current, a 10 µm spot size and an accelerating voltage of 15 kV; K, Ca, Si, Na, and Al were analysed first for 10 s, Ti and Mg for 60 s and Fe, Mn, P, S, and Cl for 50 s. Barite was used as calibration standard for S measurements. Water contents were estimated using an iterative method for obtaining volatiles by difference (in Probe for Windows; see Roman et al. 2006). This method corrects matrix effects that can arise in the presence of H₂O in glasses by iteratively adjusting water contents and recalculating the glass composition. Typical uncertainty of the method is less than 0.5 wt% (Roman et al. 2006).

Polarized light and BSE images show great similarity in textures and mineralogy between the 2020–21 lava dome and scoria samples (Fig. 2). Both lava and scoria samples are characterized by a high crystal content (45 ± 1.9 vol%; 1 standard deviation) and a mean mode of plagioclase (26.1 ± 1.9%), orthopyroxene (8.1 ± 1%), clinopyroxene (5.9 ± 1.7%), titanomagnetite (4.5 ± 1.1%), and olivine (0.6 ± 0.4–2.2 ± 0.5%) (Figs 2a, b & 3). Given the variation between the different thin sections counted (dome: n = 5, scoria: n = 5; Fig. 3; Supplementary Table S2), crystallinity and macrocrystal modes for the dome and scoria samples analysed in this study are not distinct except for olivine, which is more abundant in the scoria. Noticeable differences also exist in olivine textures, which show reaction coronas comprised of orthopyroxene and Fe–Ti oxides in lava dome samples (Fig. 2a) compared with their more pristine but subhedral embayed appearance in scoria samples (Fig. 2b). Plagioclase crystals record strong normal zonation pattern in BSE images, in some cases weak oscillatory zoning, and generally euhedral crystal habits (Fig. 2c–e). Pyroxene crystals from the explosive phase of the eruption appear to be in textural equilibrium, with most grains showing facetted crystal faces and appearing compositionally homogeneous in BSE images. Similarly, pyroxenes in the lava dome are mostly euhedral but often show feathery reaction textures on crystal rims (Fig. S1c). Titanomagnetites frequently occur as inclusions in pyroxenes, are sub- to euhedral in shape, and show exsolution textures in lava dome samples but are homogeneous in samples from the explosive phase of the eruption.

Similar to historical eruptions (Fedele et al. 2021), plutonic xenoliths are abundant in the studied explosive samples, spanning a size range from centimetres to the sub-millimetre scale (Fig. 4). We distinguish two types of xenoliths based on their mineralogy and textures: (1) Troctolite xenoliths comprise unzoned plag–ol assemblages with minor titanomagnetite and some interstitial vesicularity (Fig. 4a). (2) Gabbronorite xenoliths (plag–opx–cpx–tmt ± ol) feature poikilitic titanomagnetite in pyroxenes, olivine rimmed by orthopyroxene reaction coronas, and strongly zoned plagioclase (Fig. 4b–d). Both xenolith types are hosted in the groundmass without chilled margins. All major phases also occur as glomerocrysts of highly variable modal mineralogy with strong resemblance in plagioclase zonation pattern to free-floating macrocrystals.

In both eruptive phases the groundmass has a hyalopilitic (microlite-rich) texture comprised of plagioclase, clinopyroxene, orthopyroxene, titanomagnetite and sparse olivine microlites with small (<10 µm) pockets of interstitial glass. The microlite population in the lava dome shows abundant reaction textures, whilst appearing pristine in explosive samples. Both eruptive phases have plagioclase crystals in the groundmass that can be subdivided into unzoned microlites (<30 µm in length) and strongly zoned microphenocrysts (>30 µm). A detailed petrological and textural analysis of the microlites reveals distinct variations of groundmass phase proportions in stratigraphically constrained samples (Frey et al. 2023). As shown by Frey et al. (2023) and Morrison-Evans et al. (2023), reaction textures in dome samples typically comprise <1 µm Fe–Ti oxide halos surrounding orthopyroxene crystals, symplectite textures in rare olivine macrocrysts and exsolution textures in Fe–Ti oxide macrocrystals. Further difference in the microlite population of dome and scoria samples are the greater abundance of orthopyroxene microlites in the lava dome and the presence of olivine microlites exclusively in explosive samples (Frey et al. 2023).

All analysed bulk-rock samples (n = 26) classify as basaltic andesite with SiO₂ contents ranging from 52.8 to 55.4 wt% (Fig. 5). Scoria samples of the explosive phase (n = 18) are slightly more mafic with mean SiO₂ contents of 54.6 wt% (±0.6 1SD) compared to dome samples (n = 8) whose mean SiO₂ is 55.2 wt% (±0.1 1SD). Bulk-rock data follow a linear trend in total alkalis v. silica space and in components such as MgO and CaO, consistent with the trajectories defined by historical eruptions of La Soufrière (Fedele et al. 2021). Based on available data (Fedele et al. 2021), whole-rock compositions of the 2020–21 eruption are on average slightly more mafic than samples from historical (AD 1718/1812; SiO₂ wt% range: 53.6–56.2) and prehistoric (AD 1440: 55.3–56.2 SiO₂ wt%; AD 1580: 54.2–56.5 SiO₂ wt%) eruptions, but do not extend to such mafic compositions as the 1902/1903 eruption (50.1–55.5 wt% SiO₂). Overall, lava dome samples from the 2020–21 eruption are slightly more silicic than the 1979 dome and scoria (SiO₂ wt% range: 54–55). The scatter in the 2021 scoria data is most likely related to different amounts of entrained xenolith material that contribute to the bulk-rock composition.

Plagioclase (n = 373 spot analyses from 52 crystals) shows a wide range of compositions with molar percentage anorthite (An; CaAl₂Si₂O₈) ranging from 48 to 96% and FeO_t contents varying between 0.42 and 0.98 wt% (Fig. 6). Crystals from lava dome and scoria samples are chemically indistinguishable (Fig. 6a). Most crystal cores define a slightly inverse sloping trend from An₈₀ to An₉₅ in An–FeO_t space with considerable variation in FeO_t contents between 0.4 and 0.9 wt% (Fig. 6b). Plagioclase macrocryst rims cluster around a median of An₆₁ and FeO_t of 0.56 wt%. The strong bimodality of high An plagioclase core and lower An rim compositions produces distinctive abrupt changes in BSE greyscale intensity (Fig. 2e, g). Microlites (<30 µm) are distinct from macrocrystal rims, with elevated FeO_t contents (0.78–0.98 wt%) and median of An₅₂. A broader range of FeO_t contents in plagioclase between 0.6 and 1.6 wt% was reported by Frey et al. (2023), but the majority of their analyses are consistent with the range reported here. The outermost microphenocryst rim compositions are similar in An to microlites but typically lower in FeO_t; An₇₀–An₉₁ cores of microphenocrysts have 0.60–0.88 wt% FeO_t. Xenolith and glomerocryst plagioclase crystals have compositions of An₆₀–An₉₄ and FeO_t = 0.49–0.92 wt%.

Orthopyroxene (opx) and clinopyroxene (cpx) macrocrystal cores and rims have a restricted range in major element composition (Fig. 7). Enstatite (En; Mg₂Si₂O₆) contents in cpx (n = 240 spot analyses) range from En₃₇–En₄₅ (Mg# 65–78) with a median of En₄₂, whereas opx (n = 321 spot analyses) compositions vary between En₆₁ and En₇₃ (Mg# 64–75) with a median of En₆₆. Al₂O₃ varies between 0.6 and 5.8 wt% (median: 1.9 wt%) in cpx and –0.5 and 2.2 wt% (median: 1.1 wt%) in opx and is correlated with TiO₂ contents, which are 0.3–1.3 wt% (median = 0.5 wt%) in cpx and 0.2–0.5 wt% (median: 0.3 wt%) in opx. Cpx Na₂O shows little variation around a median of 0.3 wt%. There is no compositional difference in pyroxenes from lava dome and scoria samples (Fig. 7a). Furthermore, outermost rim and core compositions, as well as different crystal types (i.e. macrocysts, glomerocrysts, plutonic xenoliths) are similar in dome and scoria samples, except that microlite compositions span a larger compositional range (En₃₇ to En₇₂; Mg# 59–78) (Fig. 7b, c). Frey et al. (2023) report a similar microlite opx Mg# range of 59–74 and a Mg# of 67–89 in cpx microlites.

Forsterite (Fo; Mg₂SiO₄) contents in olivine phenocrysts (n = 193 analyses) range from 62 to 82 mol% and show a bimodal distribution with average values of Fo₇₆ and Fo₆₃ for the two groups. Most analyses (88%) belong to the high Fo cluster (Fo₇₂–Fo₈₂), corresponding the bulk interior of most olivine crystals. Scarce symplectic olivine grains in the dome lava (n = 7 analyses) record Fo contents between 74 and 82 mol%. Texturally, high Fo cores are found in scoria olivine macrocrysts and in xenolith grains that have not been in direct contact with the carrier melt (Fig. 8). Crystal rims (both macrocrysts and xenolith grains) that have been exposed to the carrier melt have outermost rim Fo contents of 66 to 77 mol%, likely reflecting varying extents of diffusive re-equilibration with the melt. Two analysed remnant olivine cores surrounded by coronas of peritectic orthopyroxene have relatively constant and evolved compositions of Fo₆₅ and Fo₆₂, respectively. MnO (wt%) is negatively correlated with Fo, with higher MnO values (0–59–0.83 wt%) in the low Fo group and lower values (0.30 and 0.58 wt%) in the high Fo group. CaO contents of 0.10–0.26 wt% are indistinguishable between low and high Fo groups. Olivine crystals in xenolith interiors, remote from the carrier melt, are compositionally homogeneous (Fo_76–77). No melt inclusions were observed in the studied olivine crystals.

Fe–Ti oxides (n = 177 analyses) classify as titanomagnetites and have a restricted compositional range in scoria samples, with ulvöspinel (Fe₂TiO₄) molar fractions of 0.37 to 0.41. An exception is a single grain hosted in a xenolith with an MgO content (3.43 ± 0.25 wt%) that is slightly higher than the average of 3.27 ± 0.14 wt%. No significant variation was observed between crystal cores and rims. While Fe–Ti oxides in the explosively erupted samples are homogeneous, lava dome samples show exsolution lamellae of Ti-rich compositions hosted in Fe-rich grains (Supplementary Fig. S1; Morrison-Evans et al. 2023). No discrete ilmenite grains were found in any of the samples.

Groundmass glasses (n = 10 spot analyses), interstitial between microlites, were only analysed in the lava dome sample as sufficiently large glass patches (>10 µm) were not found in the scoria samples. Lava dome interstitial glasses are rhyolitic in composition, ranging in SiO₂ from 70.80–78.61 wt% (median: 77.97 wt%). Frey et al. (2023) describe more mafic (andesitic–dacitic) interstitial glasses in the scoria. Melt inclusions (n = 18) were analysed only in rapidly quenched scoria samples, hosted in plagioclase and pyroxene crystals. The inclusions show rounded morphologies, were glassy and devoid of bubbles. The analysed inclusions are andesitic to dacitic in composition (58.66–67.20 wt% SiO₂) with median SiO₂ of 64.06 wt%. No compositional difference is observed between plagioclase and pyroxene hosted inclusions.

The concentration of pre-eruptive magmatic volatile species (H₂O, Cl, S) was constrained by EPMA analysis of melt inclusions (n = 18) hosted in plagioclase and pyroxene crystals from the explosive phase of the eruption (Supplementary Table S8b). Water contents, calculated using an iterative approach to the volatiles by difference method with typical uncertainty of less than 0.5 wt% (Roman et al. 2006), range from 1.5 to 3.9 wt% H₂O, which is similar to but slightly lower than previously published melt inclusion data for more primitive St Vincent magmas (0.8 to 5.2 wt% H₂O; Bouvier et al. 2008). Cl contents range from 0.21 to 0.28 wt% and show a weak tendency towards increasing values with decreasing water content (R² = 0.32). No correlation is evident for Cl or H₂O contents with S in melt inclusions, which range from 44 to 266 ppm, although these low S concentrations have large errors when measured by EPMA (average c. 60%).

The composition of the carrier melt in which macrocrystals are embedded is difficult to measure directly because of the high groundmass crystallinity. For this reason, the carrier melt composition, here defined as groundmass glass + microlites, must be reconstructed using different methods. Given that these reconstruction methods have strengths and limitations that may impact thermobarometric estimates, we propagate the uncertainty resulting from different approaches using a threefold strategy to determine the carrier melt chemistry: (1) mass balance calculations based on phase proportions, bulk-rock and mineral chemistry (described in the Supplementary materials); (2) averaging of bulk groundmass composition as determined by 22 defocused beam analysis; and (3) pyroxene-hosted melt inclusion analysis.

As we demonstrate below, all reconstructed permissible carrier melt compositions indicate that the crystal cargo of the 2020–21 eruption was brought to the surface in an andesitic to dacitic melt during both the effusive and explosive eruptive phases (Fig. 9; Table S8a–c). Importantly, all reconstructed carrier melt compositions are consistent with the melt evolution trend defined by published melt inclusion data for historical eruptions from La Soufrière. Melt compositions for five thin sections reconstructed by mass balance give broadly consistent results for all major element oxides with SiO₂ contents of 61.3–63.0 wt% for the 2020–21 lava dome and slightly more mafic chemistry (59.6–61.1 wt% SiO₂) for scoria samples. Melt inclusions from the 2020–21 scoria record slightly more silicic melt compositions (58.7 to 67.2 wt% SiO_2; average = 63.8 wt%) but may have been modified by post-entrapment crystallization, as there is considerable scatter in the data. Alternatively, they may represent melt trapped at previous evolutionary stages resulting from crystallization processes rather than carrier melt itself (Kilgour et al. 2013). Groundmass analyses with a defocused beam are slightly lower compared to melt inclusions (average 62.0 wt% SiO₂) with a range of 59.3 to 64.3 wt% SiO₂ resulting from scatter of individual defocused beam analyses. The range of estimates for the scoria samples relates in part to the spread in bulk-rock compositions, where higher olivine contents derived from plutonic xenoliths may shift the reconstructed melt to more mafic compositions. Given the differences in sample size, we average the results of each reconstruction approach independently. This suggests that the carrier melt of the 2020–21 eruption was andesitic to dacitic in composition with silica contents of 60 to 64 wt%.

Reliable temperature and pressure reconstructions based on cpx chemistry require knowledge of the melt composition in equilibrium with the crystals, which is not known a priori. The restricted range of cpx major and minor element variation indicates, however, that the crystals have grown from mostly similar melts that may be represented, as a first approximation, by either the bulk-rock or the range of estimated andesitic to dacitic carrier melts (Fig. 9a). Without further testing, however, it cannot be presumed that either the bulk-rock or the groundmass composition (i.e. carrier melt) represents a plausible equilibrium liquid. To test if equilibrium conditions between the observed cpx compositions and the permissible melts prevailed, we calculated Fe–Mg exchange coefficients (i.e. $KDcpx−liqFe−Mg$⁠), which at equilibrium should vary between 0.24 and 0.3 (Putirka 2008). As a further test for equilibrium, we also calculated the mismatch between predicted and observed Diopside–Hedenbergite components in cpx (hereafter called ΔDiHd; Scruggs and Putirka 2018). The predicted DiHd component in cpx was calculated based on cation fractions in the liquid phase following Putirka (1999). Observed DiHd components were calculated using a normative scheme (Putirka et al. 1996). To establish whether these two parameters effectively constrain equilibrium conditions between cpx and melts similar to those erupted from La Soufrière, St Vincent, we compared experimental cpx–melt pairs from St Kitts in the northern Lesser Antilles (Melekhova et al. 2019; Fig. 10a, b) with known pressure and temperature conditions to thermobarometric reconstructions using the model of Putirka (2008) (their equations 30 and 33). As shown in Figure 10b, reliable pressure and temperature estimates are retrieved for experimental cpx–liquid pairs that show $KDcpx−liqFe−Mg$ in the range 0.24 to 0.3 and ΔDiHd between −0.1 and 0.1, giving confidence that these filtering criteria are suitable to recover reliable pressure and temperature estimates.

Fe–Mg exchange coefficients and ΔDiHd were determined for a range of permissible melt compositions (Fig. 10a), including the median bulk-rock composition of the lava dome and scoria, reconstructed groundmass compositions from defocused beam analysis (median, 1st and 3rd quartile) and mass balance (5 thin sections for scoria and dome each), as well as the melt inclusion compositions (median, 1st and 3rd quartile). All calculations were performed with a melt Fe²⁺/(Fe²⁺ + Fe³⁺) ratio of 0.68 to 0.78. This is equivalent to an oxygen fugacity (Gualda et al. 2012) of 1 and 2 log units above the quartz–fayalite–magnetite buffer (QFM + 1 and QFM + 2), consistent with previous estimates for La Soufrière (Heath et al. 1998a), and water contents of 3.9 wt% based on our estimate described above. Using average bulk-rock compositions as liquid yields $KDcpx−liqFe−Mg$ variations between 0.4 and 0.5 with only 2% of the cpx–liquid pairs at permissible equilibrium conditions (Fig. 10c). Andesitic to dacitic carrier liquid compositions based on mass balance and defocused beam analysis yield equilibrium conditions for 23% and 27% of the cpx–liquid pairs, respectively. The overall best fit is observed using melt inclusion compositions (median, 1st and 3rd quartile) as liquid, which results in a 48% match. Note that these percentages incorporate all cpx data (see also Supplementary Fig. S4), while Figure 10c shows a subset for which measurements of distance from crystal rims are available.

The results indicate that only a small number of the clinopyroxenes analysed could have been in equilibrium with a liquid composition equivalent to the basaltic andesite bulk rock. We therefore conclude that the bulk-rock does not sufficiently approximate the melt composition from which most of the cpx crystals formed. Clinopyroxene crystals, in contrast, may have precipitated from a liquid similar to the carrier melt. Andesitic to dacitic liquids as represented by melt inclusion compositions and groundmass analysis cannot explain the full dataset, however, adding uncertainty to pressure and temperature estimates for magma crystallization in the volcanic plumbing system. Given that glass inclusions from historical eruptions of La Soufrière show a well-defined major element oxide melt evolution trend that diverges from the array of whole rock compositions in various components (Fig. 9), potential equilibrium melts may have different compositions to the explored range. Consequently, we also used a modelling approach to reconstruct the most likely melts in equilibrium with cpx based on the relation provided by the melt inclusion trend. This modelling approach explores a much wider compositional spectrum compared to the reconstructed melts. Major element oxides were expressed as a function of silica content using three exponential models (central, lower, and upper prediction intervals); fits correspond to 50 and 65% prediction levels, to account for variability in the data (red dashed curves in Fig. 9). An exponential rather than linear model was chosen because of curvature in liquid lines of descent, although the scatter in melt inclusion data limits the importance of the choice of exponential over linear fit. Melts were calculated in increments of 1 wt% between 50 and 70 wt% SiO₂ and subsequently used to calculate Fe–Mg exchange coefficients and ΔDiHd for all combinations of model melt and cpx analyses. The resulting distribution of $KDcpx−liq(Fe−Mg)$ and ΔDiHd is shown in Figure 11a, colour contoured for melt SiO₂, and illustrates the impact of melt chemistry on cpx–liquid equilibria for La Soufrière compositions. Considering all permitted equilibrium pairs of cpx and melt calculated using this approach, a maximum is observed for melt compositions between 55 and 60 wt% SiO₂ (Fig. 11b). Depending on the model choice, possible equilibrium liquids show a maximum at basaltic andesite (upper model) to andesites equivalent to the carrier melt (lower model; Fig. 11b). From this we conclude that part of the cpx crystal population of the 2020–21 eruption originally precipitated from basaltic andesite liquids that were similar but not identical to the whole-rock composition. Equilibrium crystallization of such basaltic andesite liquids may explain the observed variation in slow diffusing minor elements such as Al₂O₃ in cpx crystals and is consistent with the opx reaction coronas found around xenolith olivines (Fig. 4b–d).

Pressures and temperatures of cpx crystallization were calculated using matching groundmass and melt inclusion liquids (n = 858), as well as the entire range of modelled equilibrium melt–cpx pairs (n = 4183), thus incorporating the uncertainty stemming from melt composition. Calculations were performed with 2 to 4 wt% H₂O in the melt, consistent with melt inclusion data for the 2020–21 eruption and with previous estimates of the water content of St Vincent magmas (Heath et al. 1998a; Bouvier et al. 2010). We used three recent thermobarometric models based on either conventional linear regression of high pressures–high temperature experimental data (Putirka 2008; Neave and Putirka 2017) or a machine-learning approach that uses an extremely randomized trees algorithm (Petrelli et al. 2020). The latter technique predicts pressures and temperatures by generating a large number of hierarchical charts (regression trees) that compare an observed cpx population with experimental data (for details see Petrelli et al. 2020).

Temperatures calculated using groundmass cpx pairs and equation 33a of Putirka (2008) indicate crystallization conditions of 962 ± 15°C (median ± 1 sd) compared with the temperature derived from machine-learning thermometry (1015 ± 10°C; Fig. 12a). These results agree with each other within relative mean standard error (RMSE) of 40–60°C (errors resulting from analytical uncertainty are much smaller) for both methods and together yield an average pre-eruptive temperature of $997−35+18∘C$ (median, ±IQR: inter quartile range, 1 sd = 29°C). If we use the modelling approach described above to allow for a greater variability of melt compositions through the modelling, we recover a median temperature of 1027 ± 38°C. Cpx crystallization pressures calculated using Neave and Putirka (2017) range from −5 to 1.5 kbar but recover mostly negative pressures and are therefore not considered further. The barometer of Putirka (2008; eqn 30) recovers a median pressure of 3.4 kbar with most estimates between 2.8 and 4 kbar (IQR). Pressures calculated using the machine-learning algorithm of Petrelli et al. (2020) indicate a narrower IQR between 2.5 and 2.9 kbar with most estimates around a median of 2.7 kbar (Fig. 12b). The larger range of the Putirka (2008) barometer likely relates to the greater uncertainty of his method (RMSE of 4.5 kbar) compared to the machine-learning algorithm (RMSE = 2.9 kbar). Taken together, these approaches indicate pressures of pyroxene crystallization c. 2 to 3 kbar. Considering the full array of possible equilibrium melts expands the range of temperatures in both approaches but recovers similar median pressures and temperatures of 2.6 kbar and 1027°C. No difference is found in thermobarometric reconstructions for cpx from the dome and scoria samples, within the range represented by the variance we observe and the uncertainty of the modelling approach. Model water contents affect the temperature and pressure estimates of the Putirka (2008) model only within a range of ±10°C and ±0.7 kbar for 1 wt% change in H₂O. The model of Petrelli et al. (2020) is not sensitive to water contents. In summary, the clinopyroxene crystal cargo of the 2020–21 eruption of La Soufrière records crystallization temperatures of around 1000°C (RMSE ±40°C) and pressures between 2 and 3 kbar, with little textural evidence of significant time to re-equilibrate at different conditions.

The cpx–opx thermometer of Putirka (2008) was used to test for equilibrium conditions and to calculate crystallization temperatures (Fig. 13). Reported uncertainty for this calibration is on the order of ±38°C (Putirka 2008). Calculated temperatures range from 894 to 1023°C with slightly higher temperatures for the scoria samples (median: 969°C) compared to the dome (median: 943°C). However, calculated Fe–Mg exchange coefficients lie mostly outside the accepted range of 1.09 ± 0.14 (Putirka 2008) for crystallization under magmatic conditions (grey shaded field in Fig. 13a), instead having low values (median ± 1 sd: 0.86 ± 0.09) more consistent with sub-solidus systems (red shaded field in Fig. 13a). Low $KDcpx−opx(Fe−Mg)$ values are observed for random combinations of outermost cpx and opx rims that were in contact with the carrier melt, touching opx–cpx grains in glomerocrysts and xenoliths, as well as directly touching profiles across opx–cpx grain boundaries (Fig. 13b). Given that pyroxene compositions in xenoliths, glomerocrysts and macrocrystals are identical, this does not imply that the grains have grown from a liquid of different composition to the carrier melt but could reflect that the Fe–Mg exchange coefficient of 1.09 ± 0.14 is not applicable in this case. This may reflect uncertainty in the estimated Fe³⁺ of the pyroxenes, which in turn affects $KDcpx−opx(Fe−Mg)$⁠. Generally, the calculated two-pyroxene temperatures broadly agree with reconstructions based on cpx–liquid thermometry but are not considered further due to the failed equilibrium test.

The plagioclase–melt thermometer of Putirka (2005) was used to further constrain pre-eruptive temperatures. Calculations using all plagioclase analyses (rims and cores) and the bulk rock composition of the dome and explosive phase show that, as for cpx, an assumed $KDplag−liq$ (Ab–An) of 0.1 ± 0.05 indicates that only about 15% of the plagioclase data could be in equilibrium with a liquid of such composition. Within this range of permissible equilibrium pairs, no clear relation is evident between textural position (high An cores v. rims), indicating that the crystals did not originally grow from a liquid represented by the bulk rock composition. To test the hypothesis that the plagioclase rims were in equilibrium with the andesitic–dacitic carrier melt and that the cores grew from more mafic liquids, we used an analogous approach to the cpx–melt model described above. The major element composition of the melt was systematically varied from 50 to 70 wt% SiO₂ using the three exponential models (Fig. 9) and the resulting $KDplag−liq$ (Ab–An) was tracked as an indicator of potential equilibrium pairings. All calculations were performed with either 4 or 2 wt% H₂O in the melt, the latter resulting in temperature estimates typically 60°C higher. These model calculations show that plagioclase rims (An_55–67) of the 2020–21 La Soufrière ejecta are consistent with the reconstructed andesite carrier melt as an equilibrium liquid for melt water contents of 4 wt%, while calculations with 2 wt% H₂O would indicate melt compositions with median of 53 wt% SiO₂ (Fig. 14a). Plagioclase cores (An>70) are consistent with crystallization from basaltic-to-basaltic andesite melts for both 2 and 4 wt% H₂O. Modelled equilibrium pairs for plagioclase cores at 4 and 2 wt% H₂O yield median temperatures of 1018°C and 1040°C, respectively (Fig. 14b). Plausible plagioclase rim-carrier melt equilibration temperatures for 4 wt% H₂O show an interquartile (IQR) range between 969 and 993°C and are cooler than for 2 wt% H₂O, which yields an IQR of 1062–1071°C. However, calculation results at 2 wt% H₂O are not consistent with the andesitic–dacitic composition of the carrier melt. No systematic differences are found between crystal types or between dome and scoria samples. In summary, plagioclase rim crystallization at 4 wt% H₂O and temperature of $982−13+11∘C$(median ± IQR) matches the andesitic to dacitic composition of the reconstructed carrier melt and maximum H₂O content determined by melt inclusion analysis, and is fully consistent with the pre-eruptive temperature obtained from clinopyroxene thermometry (⁠$997−35+18∘C$⁠).

Diffusion timescales were calculated for four profiles on three olivine crystals (sample LS21–38) from the explosive phase of the eruption (Fig. 15). Most crystals in our sample set exhibited prominent embayment textures or had thin evolved rims, indicating that olivine composition was exchanging Fe and Mg by diffusion with the carrier melt prior to eruption. Pre-diffusion initial conditions in the crystal were modelled as flat profiles in accordance with core compositions at some distance from the rim. To fully assess the uncertainty, modelling was done at three different temperatures for each profile corresponding to the median, 5th and 95th percentile of the plagioclase rim thermometry distribution (⁠$982−31+23∘C$⁠), which is in good agreement with the temperature range indicated by cpx thermometry. As the crystallographic orientation of grains was not constrained, Fe–Mg exchange for each analysed profile was modelled using the coefficients applicable for diffusion along the c-axis [001], as well as the a- [100] and b-axes [010] (Chakraborty 2010), thereby fully capturing the 6-fold uncertainty resulting from diffusion anisotropy. The best-fit models were calculated based on minimization of squared residuals between the observed and modelled profiles. Timescale estimates range from seven months to eight years with good correspondence between two macrocrysts and profiles on the opposite sides of the same crystal (Fig. 15a–c). An olivine crystal located at the margin of a xenolith (Fig. 15d) records timescales of only 1 to 14 days. The reconstructed timescales are maxima given that crystal growth may have contributed to observed zonation features.

Thermobarometric reconstructions of pre-eruptive magmatic variables by means of mineral–melt equilibrium modelling allow us to make inferences about the architecture of the magmatic feeding system prior to the 2020–21 eruption of La Soufrière. Although the governing relations of glass and melt inclusion chemistry for St Vincent permit the presence of a range of basaltic andesite to andesite liquids in the sub-volcanic system, our use of all permissible equilibrium cpx–melt pairs returns strikingly consistent thermobarometric results when compared to estimates based on reconstructed groundmass compositions. Machine-learning thermobarometry suggests that most cpx crystals last equilibrated with an andesitic–dacitic carrier melt at $T=997−35+18∘C$ and P = 2–3.5 kbar. The density–depth model of Melekhova et al. (2019) for St Vincent (Fig. S3) shows that these pressures are equivalent to 8–13 km depth, which we interpret as the location of a major intrusive complex beneath the volcano. This depth range lies between the magma accumulation pressures of 4 and 1–2 kbar (Fedele et al. 2021) estimated for historical eruptions beneath La Soufrière volcano by comparing xenolith and bulk-rock compositions with petrological experiments (Pichavant et al. 2002; Pichavant and Macdonald 2007; Melekhova et al. 2015, 2019). The mineral compositions retrieved from our 2020–21 sample set are similar to those representing historical activity, suggesting that no major changes in the architecture of the plumbing system have occurred at least over the last few hundred years within the resolution of the different thermobarometric reconstruction approaches.

Extreme chemical zonation in plagioclase crystals, a typical feature of both St Vincent and the Lesser Antilles in general (Melekhova et al. 2017; Fedele et al. 2021; Higgins et al. 2021), suggests that the feeding system is characterized by multi-level magma storage. Plagioclase crystals from the 2020–21 eruption show two distinct compositional groups, represented by crystal rims (An_55–67) and cores of An_80–95, which can be linked to shallower and deeper crystallization levels respectively. Plagioclase–melt modelling (Fig. 14a) shows that crystal core compositions are consistent with derivation from basalt or basaltic andesite liquids at temperatures of 1020–1050°C, slightly higher than inferred from cpx thermometry. Importantly, although estimates of the water content using both the volatile by difference method and H₂O analysis of primitive melt inclusions (Bouvier et al. 2010) point towards 4 wt% for La Soufrière magmas, H₂O contents during mid-crustal magma crystallization may have been significantly higher if volatile saturation was achieved at greater depth than pre-eruptive magma storage. Specifically, high water contents in crystallizing melts and resulting suppression of plagioclase crystallization are required to explain the trend of Al₂O₃ enrichment observed for lavas from La Soufrière (Melekhova et al. 2015). Petrological experiments at 4 kbar on basaltic andesite compositions from Montagne Pelée (Martinique) (Pichavant et al. 2002) show that An contents >85 typically require H₂O contents of >5 wt% at temperatures of 1025°C. Similarly, the Al₂O₃-based glass hygrometer of Parman et al. (2011) suggests a minimum H₂O content of 5.75 wt% for a maximum Al₂O₃ content of 20 wt%. Given that An content in plagioclase is more sensitive to pH₂O than total pressure (Frey and Lange 2011; Cashman and Blundy 2013), and considering typical values for water solubility in basaltic melts (e.g. Newman and Lowenstern 2002), both mid crustal (c. 4 kbar) and shallow crustal (c. 2 kbar) crystallization may explain the high An cores. BSE images, however, show that plagioclase cores typically have sharp facetted boundaries rather than progressive crystallization to lower An contents in rims, which suggests a two-stage crystallization process. Taken together, these observations suggest growth of high An cores in basaltic andesite liquids in a mid-crustal magma accumulation zone and formation of less calcic rims upon ascent and emplacement at 2–3 kbar. The consistently lower An contents of rims compared to crystal cores is inconsistent with growth from a more primitive melt.

Relatively high Ti contents in plagioclase rims compared to cores (Fig. 6c) point to either equilibrium crystallization in a melt of similar composition (Bindeman et al. 1998; Dohmen and Blundy 2014) or rapid disequilibrium growth induced by high rates of decompression or cooling (Iezzi et al. 2014; Mollo and Hammer 2017). The latter effect is clearly observed in the microlite population of the 2020–21 eruption in our dataset (see also Frey et al. 2023), which shows higher FeO_t and TiO₂ and lower An contents than measured in crystal rims (Fig. 6). In either case, the strong normal plagioclase zonation likely formed in response to changes in H₂O content resulting from decompression of ascending magma rather than drastic changes in melt composition. Strong zonation textures in plagioclase contrast with the homogeneous pyroxene population, which, based on our thermobarometric estimates and textural relation with plagioclase rims, have grown only in the shallow (2–3 kbar) part of the subvolcanic plumbing system. Conversely, olivine reaction rims and coronas, as well as touching relations of plagioclase cores and pristine olivine grains in troctolite xenoliths (Fig. 4), indicate that olivine was unstable in the andesitic to dacitic melt of the upper storage region and probably originated at the depths of plagioclase core formation, in which more mafic melt compositions prevailed. This interpretation is consistent with the findings of Frey et al. (2023), suggesting that high An plagioclase and olivine microlite assemblages in the explosive deposits correspond to a deeper magma source.

Preceding dome growth, low-level seismic signals were detected starting in November 2020, a few weeks before the dome-forming eruption began (Joseph et al. 2022). In 2021, two swarms of volcano-tectonic earthquakes detected on 23–24 March had located depths mostly shallower than 5 km. On 5–6 April 2021, shortly before the sequence of explosive eruptions began, the earthquake locations abruptly transitioned to depths of c. 10 km (Joseph et al. 2022). This seismic transition was followed by a rapid, but unquantified, acceleration in dome growth and the occurrence of banded tremor (i.e. tremor bursts separated by periods of quiescence; Cabras et al. 2012) that lasted from 8 April until the first explosion at 08:51 local time on 9 April (Joseph et al. 2022; Stinton et al. 2023). The observed progression of earthquake depths to c. 10 km observed prior to the explosive phase of the eruption agrees well with our clinopyroxene thermobarometry estimates of the topmost part of the magma source region that fed the eruption (2–3 kbar; $997−35+18∘C$⁠). We therefore suggest that the seismic swarm preceding the explosive phase of the eruption was likely an expression of increased fracturing of the overlying rock burden and continued opening of a new conduit prior to fast magma ascent to the surface. This is corroborated by the observed increase in heat flux in the lead-up to the explosive eruption (Thompson et al. 2022).

Several distinct types of ground deformation were detected by GPS and interferometric synthetic-aperture radar (InSAR) both before and during the eruption and modelled to constrain source depths by Camejo-Harry et al. (2023). Deep crustal inflation, modelled with a Mogi-source depth of 18 km, was detected in the six months before dome extrusion began and continued after lava effusion commenced. In the first few days of explosive activity, strong subsidence was observed and modelled with a Mogi-source at 6 km depth (Camejo-Harry et al. 2023). These estimates independently corroborate the interpretation based on plagioclase zonation textures and chemistry that the magma feeding system of La Soufrière is multi-layered. The inferred subsidence source depth preceding the explosive phase is shallower compared to our depth estimate for the magma source based on cpx barometry (mean: 9.5 km depth; preferred range: 8 to 13 km) but overlaps with the deflation source within uncertainty of both methods. Although petrological analysis combined with experimental data (Pichavant et al. 2002) permits only an approximation of the depth estimate for the deeper magma storage region, the duality of crustal deformation signals hints towards a possible link to the strong bimodality in plagioclase zonation that is characteristic of La Soufrière, St Vincent magmas. A Mogi-source depth of 18 km for pre-eruption inflation translates to c. 5 kbar pressure (Melekhova et al. 2019). We suggest that this depth corresponds to the lower storage region where plagioclase cores and troctolite xenoliths originate. This comparison shows the power of combining petrological and geophysical studies, although further work is required to better resolve the geometry of the La Soufrière, St Vincent magmatic feeding system.

Mafic magma recharge, remobilization of previously intruded magma, reactivation of the system by volatile fluxing and incoming felsic melt could all be invoked to explain the renewed activity at La Soufrière. These different scenarios can be tested using our mineral chemical data. Recharge of magma bodies with (typically) more primitive melt is frequently observed at arc volcanoes and widely recognized as a process that can generate sufficient overpressure to initiate volcanic eruptions (e.g. Morgan et al. 2006; Ruprecht and Wörner 2007; Kent et al. 2010; Degruyter et al. 2016; Morgavi et al. 2017; Ubide et al. 2019; Weber et al. 2020; Caricchi et al. 2021). As discussed above, however, there is no evidence in crystals of a temperature increase or introduction of mafic melt that would suggest an important role of magma recharge prior to the 2020–21 eruption (Fedele et al. 2021).

Alternatively, leftover degassed magma from eruptions in the 1970s may have been remobilized. This hypothesis is an attractive explanation of the initial lava coulée formation, particularly as, in contrast to the explosive phase, only low levels of seismic and degassing unrest were detected prior to dome extrusion that might indicate the formation of a new conduit and syn-eruptive gas emissions were low (Joseph et al. 2022; Thompson et al. 2022; Esse et al. 2023). Similar scenarios of magmatic plug remobilization within a shallow conduit have been proposed at other basaltic andesite volcanoes (Utami et al. 2022). A dynamic flow model for St Vincent (Stinton et al. 2023) suggests viscosities of c. 10¹⁰–10¹¹ Pa s for the lava dome and 10⁶–10⁷ Pa s for the explosive magma. The higher viscosities were interpreted as remnant degassed and partially crystallized magma still present in the shallow conduit pushed out during the effusive phase by the more gas-rich magma that caused the explosive eruption. However, our petrological data show that the plagioclase and pyroxene macrocrystal cargo of 2020–21 eruptive products from both phases of the eruption formed at similar storage conditions and do not show any textural signs of re-equilibration at shallow depth. Remobilization of the dome-forming magma therefore requires formation of the microlite-rich groundmass of the dome magma prior to remobilization. In such a scenario, groundmass crystals may have shielded macrocrystal rims from remnant melt pools, effectively preventing re-equilibration. Independent of the triggering mechanism, remobilization of vestigial 1979 magma (as opposed to solidified rock) in the shallow volcanic conduit has thermal constraints.

Typical volcanic conduit diameters are on the order of 40–60 m with the largest conduits reaching 100–200 m (Kolzenburg et al. 2019). However, assuming a cylindrical geometry, diameters of 40–60 m would involve unreasonable conduit lengths of 14.4 and 6.4 km to satisfy the 18 × 10⁶ m³ of extruded dome lava (Joseph et al. 2022). Larger conduit diameters of 80–100 m are more compatible with the volumetric constraint and the shallow crustal seismicity (<5 km depth) observed during lava effusion (Joseph et al. 2022). Sustaining a sufficiently large volume of magma (T > 950°C) in the shallow crust (geothermal gradient = 50–150°C km⁻¹) over a period of 40 years since the 1979 La Soufrière eruption requires initial conduit diameters >140 m and length of 3–4 km (Supplementary Fig. S3). Although dyke-like geometries may result in shorter conduit length, cooling takes place mostly through the sides, which requires similar dyke thicknesses to the cylindrical conduit diameters invoked above. Remobilization of leftover magma from previous eruptions would therefore require either significant undercooling of interstitial melt below the solidus or a large (>140 m) conduit. In either case, our petrological analysis neither refutes nor supports this hypothesis but adds the conditional constraint that the dome magma must have been extruded with an interstitial melt fraction of 10 to 20% to explain the lack of textural re-equilibration of crystal rims that would be expected in protracted storage and subsequent remobilization.

The strong volcano-tectonic earthquake swarm propagated to c. 10 km depth (Joseph et al. 2022); pre- and syn-eruptive deformation occurred at about 18 and 6 km (Camejo-Harry et al. 2023). These precursory signals point towards deeper level magmatic processes as important in the initiation of the explosive eruption. From a petrological perspective, the lack of evidence for pre-eruptive re-heating or interaction of crystal cargo with a compositionally contrasting melt in the storage region (8–13 km) indicates eruption initiation by gas- or buoyancy driven internal processes that were intrinsic to the andesitic–dacitic carrier melt. The abundance of antecrystic material in the form of glomerocrysts, and plutonic xenoliths additionally indicates the importance of magma–mush interaction. We propose that a plausible mechanism for the renewed volcanic activity is mobilization of crystal-rich basaltic andesite magma at about 10 km depth, likely in response to destabilization or pressurization of a deeper level andesitic melt layer (Fig. 16). Entrainment of troctolite xenoliths in the explosive phase of the eruption with unzoned high An crystals that are similar to macrocrystal plagioclase cores, as well as the greater abundance of olivine and high An microlites in the initial explosive phase (Frey et al. 2023), suggests that there is a deep component in this phase of the eruption that may correspond to the 18 km deformation source (Camejo-Harry et al. 2023). The presence of deep-derived melt and troctolite together with plag–opx–cpx–tmt mush fragments from the storage level at c. 10 km depth, indicates an important role of melt–mush interaction.

The andesitic–dacitic carrier melt can originate by compaction-driven extraction of interstitial melt during crystallization and interaction of recharging melts or volatiles (Dufek and Bachmann 2010; Parmigiani et al. 2014; Pistone et al. 2017; Hartung et al. 2019), percolative reactive flow through crystal mush (Solano et al. 2012; Jackson et al. 2018) or partial melting of previously intruded gabbro (Melekhova et al. 2015). Regardless of the formation mechanism, however, buoyant melt-rich layers can develop gravitational instabilities in crystal mushes that cause melt ascent to shallower levels or, ultimately, to the surface (Seropian et al. 2018; Sparks et al. 2019). For andesitic melt lenses in crystal mushes with viscosity of 10¹⁴ to 10¹⁵ Pa s, instability timescales of several decades to centuries can be expected (Seropian et al. 2018), consistent with the eruptive behaviour observed at St Vincent over historical time. In summary, intrusion of deeper-derived andesite melt accompanied by disintegration of crustal cumulates and crystal mush is consistent with the petrological characteristics of the 2020–21 eruption of La Soufrière, St Vincent.

Diffusion modelling of olivine in the 2020–21 scoria reveals that crystals were entrained in the carrier melt both prior to and during eruption. Additionally, although our preliminary calculations are subject to large uncertainties, modelling results indicate that individual crystals spent different periods of time in the carrier melt prior to eruption. Timescales on the order of days to weeks are likely related to the breakup of olivine-bearing xenoliths entrained by the andesitic–dacitic melt during the eruption, as previously described from St Vincent (Arculus and Wills 1980; Tollan et al. 2012; Melekhova et al. 2019). Calculated timescales of months to years extend the interaction timescale between carrier melt and olivine cargo to the pre-eruptive period, and before the detection of seismic unrest prior to the explosive eruption. These timescales suggest that mush–carrier melt interaction is a driver, rather than a consequence, of the eruption.

Less intense seismic flurries also occurred during the pre-eruptive episode at La Soufrière (Thompson et al. 2022). Indeed, plagioclase residence times modelled for the 1979 eruption suggested entrainment of crystals in the melt 100 to 450 years prior to eruption (Zellmer et al. 1999). While these timescales are longer than our estimates, these are maxima due to uncertainty in the initial plagioclase trace element concentrations. ²³⁸U–²³⁰Th isochrons from St Vincent cumulate xenoliths suggest crystal residence times on the order of tens of thousands of years (Heath et al. 1998b). Old radiometric ages and short diffusion timescales are commonly observed in magmatic systems (Cooper and Kent 2014). Importantly, they are not mutually exclusive as in-growth of radioisotopes and diffusive re-equilibration are activated by different processes (gradients in chemical potential/temperature v. radioactive decay). For this reason, it is not uncommon for arc magmas to have crystals with old U–Th ages and young rims (Cooper and Kent 2014; Cashman et al. 2017). Although not their favoured explanation, Heath et al. (1998b) showed that remobilization of old intruded crystals by a ‘zero-age’ melt can explain the U–Th isochron data. A more detailed study of olivine diffusion timescales in all phases of the recent eruption is required to assess whether olivine age distributions can be related to large regional earthquakes that may have contributed to the destabilization of the magma feeding system and mobilization of melts and cumulate piles beneath St Vincent.

Petrological evidence for the transition from effusive to explosive activity during the 2020–21 eruption is remarkably sparse in the analysed sample set. Our finding that macrocrysts mineral chemistry, pressure and temperature estimates are indistinguishable for the dome and scoria samples, together with the similar macrocryst contents of the erupted rocks, suggest no major difference in storage conditions prior to the final ascent during the explosive eruptions. Given these similarities, it seems unlikely that differences in magma viscosity at depth or initial volatile content caused the effusive–explosive transition at La Soufrière. However, as noted above, the observed eruption transition may record when dome magma in the shallow conduit at high crystallinity, and in a mostly degassed state, was extruded by deeper gas-rich magma, which caused a transition in eruptive style once the deep magma reached the surface (Sparks et al. 2023; Stinton et al. 2023). However, the age of this earlier magma is not constrained by the petrological data and may have been intruded subsequent to the last eruption in 1979.

Alternatively, gas loss during slow aseismic ascent of the dome magma from the storage zone at 8–13 km, followed by rapid (gas-retaining) ascent could explain the common characteristics of microcrystal assemblages in the effusive and explosive phase. Here the contrasting intensity of seismic unrest prior to both phases of the eruption is explained by differing magma ascent rates, as evidenced by the rapid increase in lava extrusion rate prior to the onset of explosive activity (Stinton et al. 2023) and differences in the microlite populations of rocks erupted during different phases of the 2020–21 eruption (Frey et al. 2023). Limited seismicity during both the precursory and initial stages of lava extrusion is still consistent with ascent of the dome magma from a crustal reservoir, given observations of aseismic conduit formation at other arc volcanoes (Roman and Cashman 2018). This poses the question of why ascent rates of magmas with equivalent storage conditions would be different. We suggest that differences in ascent rates could be related to the architecture of the multi-layered plumbing system, as evidenced by the presence of deep-derived troctolite xenoliths with unzoned high An plagioclase in explosive products and deformation sources at c. 18 and 6 km depth (Camejo-Harry et al. 2023). Pressurization of the deep source may have mobilized the shallow dome-forming magma; subsequent rapid ascent of deeper gas-rich andesitic–dacitic melt could have produced the later explosive activity. This two-stage model is similar to remobilization of leftover magma (from the 1970s or later) in the conduit but does not require formation of groundmass microlites prior to lava mobilization. Using the microlite assemblage from earlier explosions in the 2021 sequence, Frey et al. (2023) investigated late-stage variations in storage and ascent that could have led to these explosions.

Crucially, factors unrelated to the degassing behaviour of the magma are hard to reconcile with the observations. Although the erupted bulk-rock compositions of the explosive phase are slightly more mafic than the dome, much of this variation is controlled by the presence of partially disintegrated olivine-rich xenoliths in the latter and may simply represent sampling bias. Bulk-compositional control on eruptive style also seems unlikely because samples from the 1979 eruption record a similar chemical spread but effusive activity produced more mafic bulk-rock compositions. Melt compositional controls on viscosity (Di Genova et al. 2017) are also unlikely in crystal-rich basaltic andesite magmas. In contrast, textural differences between dome and explosive samples (Fig. S1), including reaction rims on olivine and pyroxenes, reacted microlites in groundmass, and exsolution of Fe–Ti oxides in the lava dome samples, can be attributed to slow cooling and oxidation during dome emplacement (Saito et al. 2007; Scott et al. 2012; Morrison-Evans et al. 2023). Given the similarity of macrocrystal compositions in effusive and explosive samples, it seems plausible that the transition in eruptive style is related to differences in ascent rates and outgassing efficiency (e.g. Cassidy et al. 2018). For example, if ascent is buoyancy-driven from a source region at c. 10 km depth, then progressive thickening of this layer will increase buoyancy forces and increase ascent rates. We speculate that thickening of a buoyant layer by injection of andesitic melt in late 2020 and early 2021 may have played a ‘bottom-up’ role in driving the transition in eruptive style, although we cannot rule out a ‘top-down’ process such as a decompression wave propagating down the conduit in response to continued extrusion of degassed lava. Indeed, limited olivine diffusion timescales of <14 days bracket the seismic precursors to the explosive event; improved spatial and temporal sample density could further test and constrain these models. In summary, differences in magma ascent rates, likely reflecting the multi-layer architecture of the plumbing system, are consistent with the effusive to explosive transition of crystal-rich basaltic andesite magma of similar composition and physical properties. Future work could focus on better defining the relationship between olivine diffusion timescales and geophysical monitoring signals, as well as constraining the budget of volatiles such as the CO₂ content for St Vincent magmas.

1. The 2020–21 eruption of La Soufrière, St Vincent is characterized by both slow extrusion and explosive eruption of crystal-rich basaltic andesite magma. The eruptive products are similar in bulk-rock composition and mineral chemistry to previous activity of the volcano, indicating that La Soufrière volcano is currently in a steady-state regime.

2. The eruption was fed by magma at about 1000°C from a reservoir between 8 and 13 km depth. Plagioclase zonation patterns indicate discrete transfer between deeper (i.e. about 18 km) and shallower levels (c. 8–13 km depth) of the magmatic system and show no evidence for chemical interaction or mixing of mafic and evolved melts. The shallow part of the intrusive complex beneath St Vincent is likely comprised of predominantly sub- or near-solidus material.

3. The crystal cargo was brought to the surface in an andesitic–dacitic carrier melt with water content of ≥4 wt%. Disruption of the shallow mush reservoir by the incoming andesite is suggested by the common occurrence of xenoliths and glomerocrysts. Seismic activity extending to 10 km depth may document breakup of previously intruded and mostly solidified magma batches followed by fracturing of the overlying rock burden.

4. Olivine crystals resided in the melt for months to years prior to the eruption and were entrained syn-eruptively by disintegration of xenoliths and show that intrusion of andesite melt into the shallow part of the feeding system occurred shortly prior to the eruption. The impact of andesite melt on the buoyancy of the system may have played an important role in its destabilization.

5. Lava dome and explosively erupted samples are nearly identical in phase proportions and mineral chemistry with no petrological evidence for protracted sub-surface storage and remobilization of dome-forming magmas. Differences in eruptive style at St Vincent likely reflect outgassing efficiency during magma ascent occurring at differing rates over the course of the eruption, possibly in response to changes in buoyancy forces in the source region at depth.

We would like to thank S. Sparks for valuable suggestions and discussion. DMP acknowledges support from the NERC Centre for Observation and Modelling of Earthquakes, Volcanoes and Tectonics (COMET). The authors are grateful to R. Robertson and the UWI Seismic Centre for discussion and support in the field. We thank V. Smith and A. Matzen for help with electron microprobe analysis, M. Beverley-Smith for rapid thin section preparation, and R. Galbraith for preliminary petrography. We also thank B. Andrews (Smithsonian) for discussion.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

GW: conceptualization (lead), data curation (lead), formal analysis (lead), investigation (lead), methodology (lead), software (lead), visualization (lead), writing – original draft (lead), writing – review & editing (lead); JoB: conceptualization (supporting), funding acquisition (equal), investigation (supporting), writing – review & editing (equal); JeB: conceptualization (supporting), funding acquisition (equal), investigation (supporting), project administration (equal), writing – review & editing (equal); DMP: conceptualization (supporting), funding acquisition (lead), investigation (supporting), project administration (lead), writing – review & editing (supporting); PC: conceptualization (supporting), funding acquisition (equal), investigation (supporting), writing – review & editing (supporting); HF: conceptualization (supporting), investigation (supporting), writing – review & editing (supporting); MM: investigation (supporting), writing – review & editing (supporting); BVD: investigation (supporting); KC: investigation (supporting), writing – review & editing (equal).

This study was supported by the UK Natural Environment Research Council (NERC) through an Urgency Grant (NE/W000725/1) to DMP. GW acknowledges funding through an Early Postdoc. Mobility Fellowship from the Swiss National Science Foundation (SNSF). JB acknowledges funding from a Royal Society Research Professorship (RP\R1\201048).

All data are available in the electronic supplementary materials to this article.