Estimation of eruption source parameters for the 2021 La Soufrière eruption (St Vincent): implications for quantification of eruption magnitude on volcanic islands Open Access

Supplementary material: Additional figures and data tables in support of this study are available at https://doi.org/10.6084/m9.figshare.c.6693413

Airborne tephra and the associated ground deposits represent the widest-reaching volcanic hazard associated with explosive eruptions (e.g. Bonadonna et al. 2021). Tephra particles are injected by volcanic plumes high into the atmosphere from where they are advected by wind and transported vast distances (e.g. Sparks et al. 1997; Rose and Durant 2009; Bonadonna et al. 2015; Pyle 2016). Particles fall to the ground according to their settling velocities and atmospheric conditions; the resulting deposits show thinning and fining trends with increasing distance from the vent (e.g. Sparks et al. 1997; Volentik et al. 2010; Bonadonna and Costa 2013; Houghton and Carey 2015). The rate of change in deposit thickness and grain-size distribution allows estimation of the tephra volume, plume height, umbrella cloud radius and total grain-size distribution (i.e. eruption source parameters (ESPs)) using various statistical or numerical methods (e.g. Carey and Sparks 1986; Pyle 1989, 2015; Fierstein and Nathenson 1992; Bonadonna et al. 2005; Bonadonna and Costa 2012; Rossi et al. 2019; Constantinescu et al. 2022). Subsequently these ESPs are used to characterize the eruption using magnitude and intensity scales (e.g. Newhall and Self 1982; Pyle 2015; Constantinescu et al. 2021).

Access to well-preserved deposits encompassing the full spatial extent of the tephra fallout is necessary for the inversion methods to yield precise/highly resolved estimates of the ESPs. Unfortunately, few deposits are well preserved and mapped in sufficient detail due to geographical and environmental conditions that do not favour their preservation (e.g. Pyle 2016). For most prehistoric eruptions, deposits are often buried, subjected to post-deposition modifications, or eroded, resulting in a geological record biased toward large eruptions (e.g. Engwell et al. 2013; Kiyosugi et al. 2015; Pyle 2016). Deposits of ongoing eruptions can be easily eroded by heavy rainfall or strong winds, prompting volcanologists to conduct immediate sampling campaigns in highly hazardous conditions. As many subaerial active volcanoes are located near oceans or on volcanic islands, tephra fallout can occur mostly offshore, rendering the deposit almost completely inaccessible (e.g. Bonadonna et al. 2002; Houghton and Carey 2015; Pyle 2016; Primerano et al. 2021). As a result, critical information about deposit geometry (i.e. extent, thinning rate) and the total grain-size distribution is lost, yielding an ill-posed inverse problem when trying to estimate ESPs. This is evidenced by different combinations of ESPs yielding (nearly) the same simulated tephra deposit (Connor and Connor 2006; White et al. 2017; Constantinescu et al. 2022). Volcanological investigations often rely solely on the proximal (e.g. <6 km from the vent) and very proximal deposits (e.g. <4 km from the vent) that present high variability due to the near-vent plume dynamics, compaction and burial by subsequent deposits, syn-eruptive wind and water erosion (e.g. Engwell et al. 2013; Buckland et al. 2020). The variability of proximal deposits and the lack of measurements in the medial and distal facies contribute to the aleatoric and observational uncertainty (e.g. Bonadonna et al. 2002; Engwell et al. 2013; Connor et al. 2019; Constantinescu et al. 2022). Owing to their simplifications (e.g. simplified physical models) and assumptions (e.g. extrapolated thinning trends), statistical or numerical models used for ESP estimation contribute with epistemic uncertainty to the results (e.g. Biass et al. 2019; Mannen et al. 2020; Yang et al. 2021; Constantinescu et al. 2022). The potential for large uncertainties in estimated ESP values may affect eruption characterization and classification, and ultimately impact future hazard assessment, therefore robust uncertainty quantification (UQ) is of paramount importance (e.g. White et al. 2017; Constantinescu et al. 2022).

The question arises: how well can erupted mass/volume and related ESPs be estimated, with uncertainty quantification, given the realities of the limited sample distribution constrained by geography and potentially corrupted by post-depositional processes (e.g. compaction, erosion) and complex near-vent physical processes not included in the models? To investigate these effects, the 2021 eruption of La Soufrière volcano (LSV) (Saint Vincent and the Grenadines (SVG)) is used as a case study. Owing to the prevailing winds, the tephra fallout occurred predominantly eastward, over the Atlantic Ocean (e.g. GVP 2021a, b; Horváth et al. 2022; Taylor et al. 2022; Cole et al. 2023; Supplementary material Fig. S1). On land, the bulk of the tephra deposit is confined to the northern side of the island, between the vent and the shoreline (<6 km distance; Fig. 1). The first thickness measurements and sampling of the tephra deposit occurred immediately after the explosions ended (late April–early May 2021) in order to maximize the chance of sampling the deposit before it was eroded or compacted by the heavy rainfall. At the time, no measurements or samples were collected from the very proximal facies (<4 km from vent) given the potential for unexpected volcanic activity and inaccessibility of the terrain. A few months after the eruption ended, and after a rainy season (i.e. January 2022), several deposit thickness measurements on the flanks of the volcano were added to the dataset. The spatial distribution of tephra thickness measurements is constrained to a narrow band near the shoreline, generally between 4 and 6 km from the vent, with few measurements on the southeastern flank of the volcano (Fig. 1). As described previously, we expect this tephra dataset to provide incomplete information in the inverse problem to strongly condition/inform all ESPs. In other words, the available data may not be sufficient to constrain the posterior parameter ESP estimates with reasonable uncertainty.

Here we couple Tephra2 (Bonadonna et al. 2005; Connor and Connor 2006), an advection–diffusion model used to calculate tephra mass accumulation per unit area (kg m⁻²), with a highly efficient Bayesian inversion and UQ algorithm in the form of a model-independent iterative ensemble smoother tool (White 2018), to estimate the posterior probability density functions of ESPs that best describe the observed variations in the accumulation of the tephra associated with the first two explosive phases of this eruption (i.e. U1 and U2; Cole et al. 2023). We invert deposit thickness/mass load to estimate the erupted mass, average plume height and the total grain-size distribution (TGSD) for each of the two phases. Apart from plume heights and mass discharge rates that were directly measured during the eruption, the erupted mass and the TGSD do not have direct estimates and must therefore be determined through inversion. Both erupted mass and TGSD rely on the observation of the thinning and fining rates of the deposit as a function of distance away from the vent and this missing information (i.e. missing thickness and grain-size data from the medial and distal deposit) represents the main source of uncertainty in the posterior estimates, along with the assumptions made in the forward model. The procedure adopted here is effective for quantifying the uncertainty of the ESPs that inevitably results from missing deposit information and choice of forward model.

St Vincent (SV) is a volcanic island located in the southern part of the c. 750 km long Lesser Antilles Volcanic Arc, developed by the subduction of the South American Plate under the Caribbean Plate starting from the Late Oligocene (c. 25–20 Ma) (Fig. 1) (e.g. Rowley 1978; Macdonald et al. 2000; Robertson 2005; Germa et al. 2011; Cole et al. 2019; Fedele et al. 2021). La Soufrière (1230 m a.s.l) is the only active volcano on the island of SV. Recent activity (<600 years BP) includes both effusive (i.e. basaltic andesite dome emplacement) and explosive (i.e. pyroclastic eruptions) behaviour. Seven confirmed eruptions are documented, some of which severely affected SVG and other surrounding islands: the 1440, 1580, 1718, 1812, 1902–03, 1971 and 1979 CE events (e.g. Anderson 1908; Aspinall et al. 1973; Shepherd et al. 1979; Shepherd and Sigurdsson 1982; Robertson 2005; Pyle 2017; Pyle et al. 2018; Cole et al. 2019; Fedele et al. 2021). After the 1979 eruption, the volcano showed no sign of eruptive activity until 2020.

The first signs of reactivation at LSV magmatic system occurred in early November 2020 with the onset of anomalous seismic activity (Joseph et al. 2022). This precursory seismic activity increased, and a thermal anomaly was observed by satellites in late December 2020 and greyish-white emissions from the summit area were seen by communities on the SW of the volcano (Joseph et al. 2022). The eruption began on 27 December 2020 with a new lava dome forming inside the crater, located in the WSW sector of the crater floor, adjacent to the 1979 dome (Joseph et al. 2022). Effusive activity continued for three months with an increase in seismic activity in late March 2021 and a rapid increase in effusion rate in early April leading to the explosive phase of the eruption. The first explosion occurred at 12:41 UTC on 9 April and was followed by a period of sustained, but pulsing, explosive activity from 16:00 UTC on 9 April to 06:00 UTC on 10 April (Joseph et al. 2022).

Explosive activity continued for 13 days, with a series of explosions grouped in four stages based on the intensity and duration of the individual explosive events and estimated average magma flux (e.g. Sparks et al. 2023). The first stage of explosive activity started in the morning of 9 April 2021 and continued intermittently throughout the next c. 20 hours. Several explosions created volcanic plumes that rose to 4–16 km a.s.l and tephra fallout was reported across SV and as far as Barbados (Fig. 1) (GVP 2021b; Joseph et al. 2022).

Between 10 and 11 April, the second stage comprised a series of nine strong explosions. The volcanic plumes rose to 12–16 km a.s.l. and tephra fallout occurred island-wide, reaching SSW to the Grenadines and to St Lucia to the NNE (GVP 2021b; Joseph et al. 2022).

The third stage corresponded with a change in seismicity pattern, fewer explosions with an increasing duration of repose intervals between them and decelerating magma discharge rates (Cole et al. 2023). On 12 April 2021, intense explosions were associated with ash venting, creating plumes up to 12 km a.s.l and the generation of ash-rich pyroclastic density currents (PDC) deposits. Discrete explosions continued to occur until the end of 14 April 2021 during which PDC deposits were emplaced in the valleys on the southern and southwestern flanks of the volcano and tephra fallout continued. The last stage was characterized by weaker plumes and longer repose intervals with only three explosions occurring between 15 and 22 April 2021 (Cole et al. 2023).

According to the National Oceanic and Atmospheric Administration (NOAA) database the wind profiles at the time of the explosive phase of the eruption favoured tephra dispersal mainly east over the Atlantic Ocean. The low-level winds (0–4 km a.s.l.) mostly blew to the west and SW while the mid-level winds (8–18 km a.s.l.) blew towards the east between 30° and 120° with a wind shear between 0 and 6 km a.s.l. Above 18 km the winds blew toward the west. Wind velocities for the low-level winds ranged between 0 and 5 m s⁻¹ while the mid-level winds reached up to 15 m s⁻¹. As most of the volcanic plumes reached up to 16 km a.s.l. (e.g. Horváth et al. 2022) the mid-level winds dispersed tephra predominantly to the SE, east and NE with significant deposition on the windward (east) side of SV (Fig. 1). Low-level winds favoured the dispersion of tephra to the NW, west and SW on the leeward (west) side with wind shear allowing the fine ash to reach as far south as Kingstown (Supplementary material Fig. S1). Fine ash from the most vigorous plumes reached as far east as Barbados (>150 km) with some trace ash reported in St Lucia.

By 22 April 2021, tephra fallout from the 32 explosions produced significant damage to the northern side of the island. After the last explosion on 22 April 2021, the seismicity levels dropped and remained low, marking the end of the eruption (GVP 2021c).

The explosive sequence and the resultant pyroclastic deposits have been described in detail by Cole et al. (2023); here we summarize the characteristics of the tephra fall deposit. The stratigraphy of the tephra fall deposits consists of several distinctive units of interbedded layers of lapilli and ash (Fig. 2). Due to the very short repose intervals between some explosions, their pulsating character and varying wind conditions, individual very thin layers within these units are visible in some locations but indistinguishable in others. At some locations, the sequence is incomplete due to the variable intensity of the plumes and variability of the wind profile. In general, the lowermost tephra units (U1 and U2) are the thickest, best preserved and are found on both windward and leeward coasts of SV.

The lowermost unit (U1) resulted from several explosions between 9 and 10 April 2021 and it reflects a NE dispersal axis. The maximum measured thickness of this unit is c. 20 cm close to the crater rim (Cole et al. 2023). On the windward coast the thickness ranges between 3 and 6.5 cm and on the leeward side between 0.2 and 5 cm (Supplementary material Fig. S2). The U1 unit consists mainly of fine- to medium-grained scoria lapilli, coarse ash and interspersed dense lithics up to 1 cm. Layering and grading can be discerned where the deposit is thickest.

The second unit (U2) is rich in medium to coarse ash with relatively uniform thickness on both sides of the island and is associated with a series of explosions between 09:35 UTC and 16:20 UTC on 10 April 2021 (Cole et al. 2023). At this time, fine ash fell as far south as Kingstown and even reached Barbados. This ash-rich unit, interbedded with coarse scoria lapilli and lithics (c. 0.5 cm), comprises several brownish and greyish layers (7–8 layers in proximal locations) and reaches a maximum of 22 cm closer to the crater rim and 6.5 cm in thickness in measured locations on the windward side (Supplementary material Fig. S2).

Distinctive eruptive pulses associated with longer repose intervals between 18:30 UTC and 21:20 UTC on 10 April 2021 are recorded in the alternating ash and lapilli layers of the third unit (U3). This sequence is identifiable as a ‘couplet’ with individual thin ash and fine lapilli layers between two layers of coarser lapilli. Individual layers are millimetres thick with the whole unit reaching a maximum of 2.5 cm.

The uppermost layers (U4 through U7) are associated with the explosive phases between 23:00 UTC on 10 April 2021 and 15:09 UTC on 22 April 2021 (Cole et al. 2023). U4 is an ash layer with accretionary lapilli (c. 8 mm) that thins from c. 12 cm in the proximal facies to a few millimetres at >4 km away from the vent. An uneven lapilli layer, U5, reaches between 2 and 5 cm in the proximal locations on the southwestern and southeastern flanks. This lapilli layer is overlain by the fine-ash sequence of U6 which is found only on the southwestern flanks of the volcano. U7 is identified only at a few localities on the southwestern flanks of the volcano and consists of 3 to 4 very thin, fine-grained, lapilli layers (Cole et al. 2023).

Thickness measurements for the two units were conducted at 64 and 78 locations, for U1 and U2 respectively, on the leeward and windward coasts, with very few measurements inland given the proximity to the active cone or due to the inaccessible rough topography. Two measurements for U2 were taken in Barbados (c. 150 km east from the vent).

The measured localities generally follow the shorelines on both coasts in 10 to 15 km long transects and on the southeastern flank (Fig. 1). The thickness of these units varies between 0.1 and 23 cm and they are consistently thicker on the windward coast along the observed main axis of dispersion (Supplementary material Fig. S2 and Table S1). Eleven samples were collected for grain-size analysis (Fig. 1; Supplementary material Table S1). Samples were dry-sieved down to 7ϕ (31 µm) at 1ϕ intervals and show a distribution between −3ϕ and 7ϕ.

Tephra2, an advection–diffusion model for tephra sedimentation from subvertical plumes is used as a forward model to estimate tephra mass accumulation on the ground (Bonadonna et al. 2005; Connor and Connor 2006). We refer the reader to Connor et al. (2019) and Bonadonna et al. (2005) for detailed descriptions of the model. In summary, the model simulates the release of particles of different sizes from a subvertical source above the vent that resembles the volcanic plume (Connor et al. 2011). The shape of the plume and the distribution of the mass of particles within the plume are controlled by a Beta distribution (Connor et al. 2019) so that the model can simulate a wide range of configurations of the vertical distribution of tephra within the source area above the vent. Once released, the particles fall according to their settling velocity within the specified vertically-discretized wind profile (e.g. Bonadonna et al. 1998; Connor et al. 2011). Their final (x, y) position on the ground is controlled by particle characteristics (e.g. size, shape, density) and atmospheric properties (e.g. atmospheric diffusion, wind profile). For simplicity, several assumptions are made: particles are spherical and fall on a flat topography given as an average altitude around the volcano; the horizontal atmospheric diffusion is uniform with no vertical diffusion.

The forward tephra fallout model is coupled with pestpp-ies, a Bayesian inversion and UQ algorithm (White 2018). Pestpp-ies is a model-independent iterative ensemble smoother using an ensemble to approximate the Jacobian matrix in the Gauss–Levenberg–Marquardt (GLM) scheme. The ensemble smoother algorithms have been used extensively for parameter estimation and UQ in various problems across the natural sciences and are particularly well suited to coping with high-dimensional, ill-posed, non-linear inverse problems, including reconstructing tephra dispersion and sedimentation (e.g. Crestani et al. 2013; Emerick and Reynolds 2013; White 2018; Mingari et al. 2022a, b). Here we chose to use pestpp-ies for its ability to efficiently quantify the posterior uncertainty in models with a nonlinear relation between ESP quantities and the resulting simulated spatial distribution of tephra. This is achieved by the inversion through minimizing the error between the observed and modelled deposit.

We invert tephra mass per unit area (kg m⁻²) to estimate the erupted mass, average particle release height and total grain-size distribution. Mass load at each sampled location was obtained from thickness data using an average bulk deposit density of 1500 kg m⁻³ (Cole et al. 2023). Average particle release height is preferred to plume height, since U1 and U2 layers represent deposition from several plumes and explosive pulses that reached different heights between 12 and 18 km a.s.l. (e.g. GVP 2021a; Horváth et al. 2022). TGSD is a necessary input parameter for tephra dispersal models (e.g. Volentik et al. 2010; Bonadonna and Costa 2013) and, to our knowledge, no prior TGSD estimation for U1 and U2 was available at the time of this study. We therefore estimate TGSD through the inversion method and acknowledge the associated uncertainty.

The inversion procedure estimates the posterior uncertainty in each parameter using an ensemble of ESP values. The process begins by generating a Prior ESP ensemble, sampled by the user-defined prior probability density function (PDF) of each ESP. The prior ESP PDF is informed by expert knowledge of the eruption and constrained to physically plausible values (e.g. Constantinescu et al. 2022), but still contains considerable uncertainty in as much as each ESP is unknown prior to inversion. In other words, we guide the inversion to search for optimal solutions by specifying a search domain bounded by a minimum and maximum value we assume a parameter can take. For each ESP, an initial assumed value represents the mean of the prior PDFs and the distance between the bounds is proportional to the standard deviation.

The prior ESP PDFs for the erupted mass, TGSD mean and average particle release heights of both U1 and U2 units are reported in Table 1 and a list of all inputs used in the inversion is provided in Supplementary material Table S2. The boundaries for the prior distributions of the erupted masses were set to include the tephra volumes obtained by Cole et al. (2023) but wider in order to account for the uncertainty introduced by the simplifications assumed in the forward model (e.g. spherical particles, 2D wind field, no vertical diffusion). The average particle release height was set between 12 and 16 km a.s.l., assuming sedimentation occurs at between c. 70% and c. 90% of the maximum plume tops (e.g. Horváth et al. 2022). The grain-size analysis of the 11 samples collected from SV was used to define a prior PDF for the mean grain size of the deposit. We highlight that this grain-size information is used only to define the prior PDF and is not used explicitly in the inversion analysis (i.e. we do not invert the mass of individual grain-size classes at each sampled location), which is only conducted on deposit thickness/mass load.

In practice, a non-linear search process iterates an ensemble of combinations of ESPs from the parameter prior distributions through repeated approximate linearizations of the GLM algorithm until an acceptable misfit between the modelled and observed mass loading is found. A criterion established a priori, the weighted Residual Sum of Squares (RSS), is used to evaluate the model's goodness-of-fit (White et al. 2017; Constantinescu et al. 2022). In an effort to increase the efficiency of the inversion (i.e. minimize the error between the observed and calculated deposit), an additional rejection criterion was applied for each iteration of the inversion algorithm to remove non-behavioural realizations of ESP values (i.e. realizations that fall outside an acceptance ratio of 1.25 times the mean RSS of the entire ensemble). Conceptually, the ensemble analysis quantifies the posterior uncertainty in the ESPs that remains after assimilating the (uncertain) tephra thickness data. In this way, the posterior ESP ensemble accounts for both prior uncertainty in the ESPs, as well as the expected noise in the tephra thickness measurements.

The wind data used for the inversion were downloaded from the NOAA NCEP/DOE Reanalysis 2 database (https://psl.noaa.gov/data/gridded/data.ncep.reanalysis2.html (Kanamitsu et al. 2002). These data, as a daily mean, are provided in height, velocity, direction format for 17 atmospheric levels. We acknowledge that by treating the wind data as a known quantity (i.e. wind speed and direction vary continuously whereas these measurements represent a point in time) may contribute to the biases in the final ESP estimates (e.g. White et al. 2017); future work will explore joint inversion of ESP and wind profile quantities.

The results from the inversion analysis are reported in Figure 3 and Table 1 for erupted mass, average particle release height and TGSD mean. The results are evaluated by how much the posterior uncertainty in parameter estimates is reduced when compared with the uncertainty in the prior distributions. The uncertainty range is reported as ±2 σ for both prior and posterior distributions.

Generally, the relation between prior and posterior ensembles is coherent with the expected information content of the tephra data used in the inversion, and the posterior PDFs for each ESP are within regions of substantial prior probability mass.

For the erupted mass of U1, the uncertainty in the posterior estimates was reduced by c. 60% when compared to the prior distribution (Fig. 3; Table 1). A mean erupted mass was estimated at 3.3 × 10¹⁰ kg (or a volume of 2.2 × 10⁷ m³ for 1500 kg m³ bulk density) with an uncertainty range between 2.3 × 10¹⁰ and 4.4 × 10¹⁰ kg (Fig. 3).

The mean erupted mass for U2 was estimated at 3.1 × 10¹⁰ kg (uncertainty range between 2.8 × 10¹⁰ and 3.3 × 10¹⁰ kg). The uncertainty in the posterior estimate was reduced by c. 90% when compared to the prior (Fig. 3).

The posterior estimates for the average particles release height yielded a mean of c. 13.6 and c. 13.3 km a.s.l for U1 and U2 respectively, with uncertainty ranges between c. 12.7 and c. 13.9 km a.s.l for U1 and between c. 13.1 and 14.1 km a.sl. for U2 (Fig. 3).

Figure 3 shows the estimated mean of the TGSD for both U1 and U2. The best-fit model for U1 indicates an estimated TGSD mean of −0.55ϕ with a standard deviation of 2.31ϕ. The uncertainty range (−0.96ϕ to −0.15ϕ) shows an uncertainty reduction of c. 70% in the posterior when compared to the prior. The 1.14ϕ TGSD mean for U2 unit indicates a finer deposit, consistent with observations. The standard deviation calculated for this unit is 2.36ϕ. Figure 4 shows the TGSD for each of these two units as estimated by the inversion procedure. The simulated phi intervals between −3ϕ and 7ϕ was set based on the 11 samples that were collected for grain-size analysis.

Eruption source parameters such as erupted mass/volume, plume height and umbrella cloud radius are essential criteria used to quantify the sizes of explosive volcanic eruptions and, absent of any direct observations, are usually estimated from tephra deposit data, (e.g. Newhall and Self 1982; Carey and Sparks 1986; Pyle 1989; Bonadonna and Costa 2013; Constantinescu et al. 2021). Sampling the tephra deposit thoroughly while it is still pristine may reduce the uncertainties associated with ESP estimates; a general rule of thumb being that a uniform spatial distribution of measurements across the proximal, medial and distal facies produces more accurate volume estimates (e.g. Primerano et al. 2021). In this case study, we rely only on measurements taken predominantly on a crosswind profile from the proximal deposit; measurements from the medial and distal facies are completely missing for both units, except for two measurements taken for U2 in the distal deposit, in Barbados (c. 150 km away from the vent). The major source of uncertainty in the case of the 2021 SVG eruption are the missing data from the medial and distal deposit facies. In general, this observational uncertainty along with deposit variability, contributes significantly to the estimates of erupted volume and plume height, regardless of the model (statistical or numerical) used for these determinations (e.g. Engwell et al. 2013; Bonadonna et al. 2015; Connor et al. 2019; Primerano et al. 2021).

The quality of the models’ output relies on inversion of thickness (or mass load) and/or grain-size information from across the deposit. For example, studies showed that advection–diffusion models may perform better and parameters such as plume height and TGSD are better conditioned when grain-size information is available from across the deposit (e.g. Bonasia et al. 2010; Volentik et al. 2010; Costa et al. 2016; White et al. 2017; Connor et al. 2019; Mele et al. 2020). Nevertheless, deposit thickness/mass load is the most commonly used source of information (e.g. Constantinescu et al. 2022). At SVG we rely only on deposit thickness/mass load data from the proximal deposit, and this leads to an ill-posed inverse problem, where some parameters are better conditioned than others.

The inversion method adopted here relies on the definition of a parameter search domain in the form of a prior PDF for each ESP. The model results are sensitive to the choice of prior ESP PDFs (and the associated lower and upper limits of the ESPs (e.g. White et al. 2017; Constantinescu et al. 2022). For this analysis, the prior ESP PDFs were informed by direct observations of the eruption (i.e. for average particle release height) and previous estimates of 12–18 km a.s.l. (e.g. Horváth et al. 2022; Cole et al. 2023). The erupted mass was given a prior range between 5 × 10⁹–1 × 10¹¹ kg for U1 and 10¹⁰–10¹¹ for U2. The posterior uncertainty was reduced by c. 60% (U1) and c. 90% (U2), with a best fit (i.e. mean of the posterior distributions) of 3.3 × 10¹⁰ kg for U1 and 3.1 × 10¹⁰ kg for U2. Assuming a 1500 kg m⁻³ bulk density of the tephra deposit, these erupted masses indicate volumes of c. 2.2 × 10⁷ m³ and c. 2.1 × 10⁷ m³ for U1 and U2 respectively, comparable and within the uncertainty range of those obtained by Cole et al. (2023) using isopach interpolation (i.e. c. 2.1 × 10⁷ m³ and c. 2.4 × 10⁷ m³ for U1 and U2 respectively) (Table 2). The mean posterior erupted mass values may be underestimated considering the lack of samples from the medial and distal deposit and the difficulty in extrapolating the deposit limits (e.g. Bonadonna and Costa 2012; Connor et al. 2019). One way to investigate whether the erupted mass varies significantly when distal samples are missing is to execute the inversion model for U2 without considering the two samples from Barbados. The remaining samples are strictly from SVG island, from the proximal and very proximal deposit (Fig. 1). Without using the distal data, the estimated erupted mass decreased slightly to 2.4 × 10¹⁰ kg (1.6 × 10⁷ m³) and the posterior uncertainty was reduced by c. 90% when compared to the prior. This indicates that: (i) the estimated erupted mass is linearly dependent on the mass in the forward model (e.g. Connor and Connor 2006; Magill et al. 2015; White et al. 2017; Connor et al. 2019; Constantinescu et al. 2022) and (ii) the erupted mass estimated with advection–diffusion models may be underestimated when information from the medial and distal deposit is absent. However, when the posterior PDFs of erupted mass are considered, the erupted mass estimated by these two inversion experiments are statistically very similar, in that the posterior PDFs overlap.

Overall, acknowledging the major source of uncertainty, the erupted mass seems well constrained for both units and is in accordance with the values estimated using other methods (Table 2; Cole et al. 2023; Sparks et al. 2023). In Table 2, we report the volumes of U1 and U2 converted to dense rock equivalent (DRE) using a conversion factor of 0.6 (Cole et al. 2023). The volumes estimated by Cole et al. (2023) are based on isopach interpolation while those obtained by Sparks et al. (2023) use a technique based on RSAM data and satellite observation of erupted plume heights.

For many contemporary eruptions, visual, radar or satellite observations offer rapid and accurate estimation of plume height (e.g. Arason et al. 2011; Petersen et al. 2012; Mori et al. 2016). During the LSV 2021 eruption, the observed plumes generated by the explosions reached between 12 and 16 km a.s.l., with very few short-lived plumes between 4 and 8 km a.s.l. (GVP 2021a). Subsequent analysis of satellite imagery indicated that most of the plumes spread around 15–18 km a.s.l. with very few plumes reaching up to 18–20 km a.s.l. (e.g. Horváth et al. 2022). We use 12–16 km a.s.l. as the input parameter range for the inversion analysis, assuming the maximum average particle release height would be between 70% and 90% of the observed plume tops. For both units, the mean average particle release height estimated by the inversion analysis is c. 13.5 km a.s.l. ±0.5 km. These values indicate sedimentation from a height of roughly 85% of the observed average plume tops.

Tephra sedimentation is controlled by particle fall time and atmospheric properties; therefore, the TGSD is a crucial parameter for advection–diffusion models (e.g. Bonadonna and Costa 2012; Bonadonna et al. 2015; Connor et al. 2019). Previous studies using Tephra2 combined with an inversion algorithm showed that the uncertainty in estimated ESPs can be further reduced if the inversion is conducted on the grain-size data measured at each sampled location in addition to deposit thickness or mass (e.g. Volentik et al. 2010; White et al. 2017; Connor et al. 2019). In most cases this information may not be available and ESP estimation relies only on deposit thickness or mass-load measurements. Here we aim to see how well the inversion algorithm can estimate TGSD from thickness measurements only. In the absence of grain-size data at each sampled location, the posterior uncertainty reduction is dependent on the prior distributions decided by expert knowledge of the eruption and on the limited information content of the tephra thickness data. For SVG, the prior distribution for TGSD mean were informed by the grain-size analysis of 11 samples taken from the island but these grain-size data were not used explicitly in the inversion problem. In other words, the inversion was conducted on the cumulative mass load of particles of all sizes at each sampled location and not on the mass load of individual grain-size classes at these locations. We acknowledge the fact that the grain-size distributions of the 11 samples may not be representative for the whole deposit but they can be informative. The inversion model yielded a mean grain size of −0.55ϕ and +1.14ϕ for U1 and U2 respectively, with standard deviations of 2.31 and 2.26 (Fig. 4). These values indicate coarser- to finer-ash deposits consistent with observations of the on-land deposit and the eruptive plumes.

Tephra deposit geometry is largely controlled by particle size, release height and atmospheric conditions. In simulations, advection–diffusion models can yield equally good solutions for the modelled deposit when different combinations between particle size–wind speed–release height are used (e.g. Bonasia et al. 2010; White et al. 2017; Constantinescu et al. 2022). A portion of the posterior uncertainty in the estimated ESPs may be explained by the correlation between particle position on the ground and its release point in the eruption column. As both U1 and U2 deposits are the result of sedimentation from several plumes, particles of the same size may have sedimented from different heights adding to the estimated uncertainty. In addition, the lack of distal data can lead to poorer constraints on the estimated particle release height, which rely solely on the information from the proximal deposit.

Treating the wind field as a known quantity may interact with the posterior uncertainty ESP estimates. In reality the wind profile may vary throughout the eruption. Due to the limitations of the forward model, the wind data are used in the inversion as the mean wind direction and speed for the time interval corresponding to the first two explosive stages; i.e. no time variation is assumed for the wind profile. The wind field used for the first two eruptive phases indicates a general dispersion of tephra towards NE for U1 and predominantly towards E for U2 and it agrees with observations (Supplementary material Fig. S3). We suggest that: (i) the wind field along with (ii) the assumptions made in the prior PDFs for TGSD means and average particle release heights; (iii) the spatial distribution of the samples and, (iv) the simplifications assumed in the forward model, lead to the small difference between the estimated ESP for U1 and U2. The uncertainty in the ESP estimates may be further reduced if the wind field is parameterized and treated as unknown, and subsequently estimated through inversion along with the ESPs.

The ESPs estimated in this analysis generally reproduce the observed tephra thickness measurements, as well as the results obtained by other models (Cole et al. 2023). The isomass maps calculated by running the forward model with the best-fit set (posterior means) of ESPs for both units are in agreement with the observed extension of the deposit and of the eruption plumes (Fig. 5). The footprint of the calculated isomass maps for U1 and U2 are similar and comparable to the isopach maps presented by Cole et al. (2023).

In an effort to better understand what data could be collected to reduce posterior ESP uncertainty, we attempted a different strategy to execute the inversion model. To increase the information available for conditioning the ESPs, an alternative dataset containing 24 synthetic measurements of zero mass were added to the existing field data. In other words, if the inversion model cannot be informed by locations where tephra was measured, we augment our limited tephra dataset with locations where we know, or assume based on available satellite imagery, that tephra did not reach. This is an attempt to ‘bound’ the area where we expect tephra to be (Fig. 6). The inversion yielded a higher mean posterior erupted mass for both U1 and U2 units (4.1 × 10¹⁰ kg and 3.5 × 10¹⁰ kg), when compared with the results obtained by inverting only the field data (3.3 × 10¹⁰ kg and 3.1 × 10¹⁰ kg). No significant change for the particle release height and mean TGSD is found. This crude strategy can be explored in more detail. Will fully enclosing the area where we believe tephra was sedimented with measurements of zero mass loading compensate the lack of field measurements from medial and distal deposit and inform the inversion to better constrain ESP estimation? This bounding may be achieved by using satellite imagery and expert knowledge of the eruption; however, extrapolating a line of zero deposit thickness may further contribute to uncertainty in the estimated ESPs. This approach is considered in on-going research for its potential to find a way to improve ESP estimation for similar eruptions where the bulk of the tephra deposit is missing or inaccessible.

Overall, our analysis shows that inversion methods paired with uncertainty analysis can yield relatively well-constrained ESPs for explosive eruptions for which the bulk of the deposit is missing or is inaccessible. UQ provides a better understanding of the variations in ESPs and therefore eruption characteristics and its classification. Further research is recommended to improve such inversion-based approaches to ESP estimation.

The authors wish to thank S. Sparks for the valuable discussions on the study and the manuscript. We wish to thank our reviewers and handling editor for the suggestions that improved the manuscript.

The authors declare no competing interests.

RC: conceptualization (equal), formal analysis (equal), investigation (lead), methodology (equal), software (equal), writing – original draft (lead), writing – review & editing (equal); JTW: formal analysis (equal), investigation (equal), methodology (equal), software (equal), writing – original draft (equal), writing – review & editing (equal); CC: conceptualization (equal), formal analysis (equal), funding acquisition (lead), investigation (equal), methodology (equal), software (equal), writing – original draft (equal), writing – review & editing (equal); PC: conceptualization (supporting), data curation (lead), formal analysis (equal), investigation (supporting), methodology (supporting), writing – original draft (supporting), writing – review & editing (equal); KF: conceptualization (supporting), formal analysis (equal), investigation (equal), supervision (equal), writing – original draft (supporting), writing – review & editing (equal); JB: conceptualization (supporting), data curation (equal), formal analysis (equal), investigation (equal), writing – original draft (supporting), writing – review & editing (equal); RR: conceptualization (supporting), data curation (equal), formal analysis (supporting), investigation (equal), writing – original draft (supporting), writing – review & editing (equal).

RC and CC acknowledge partial support from the NSF grant #1841928. JB and PC acknowledge support from the NERC Urgency grant NE/W000725/1 and Royal Society APEX Award APX/R1/180094.

The data and the models developed and used in this study can be found at https://zenodo.org/record/8127789 or by request to the corresponding author.