Estimation of eruption source parameters for the 2021 La Soufriere eruption (St. Vincent): implications for quantification of eruption magnitude on volcanic islands

Eruption source parameters (ESPs) used to characterize explosive eruptions are estimated from tephra deposit data using different models (statistical or numerical) and inversion approaches. The ESPs thus derived are subject to substantial uncertainties when the bulk of the tephra deposit, including information about its full spatial extent and spatial variation in grain size distribution is missing due to geographic and environmental conditions. We use an advection-diffusion model coupled with a Bayesian inversion and uncertainty quantification algorithm to investigate how ESPs can be robustly estimated given such conditions. The 2021 eruption of La Soufrière Volcano (St. Vincent and the Grenadines) is our case study. An inversion is conducted for the first two explosive phases of this eruption (U1 and U2). We estimate: an erupted mass of 3.3 x 10 10  1 x 10 10 kg for U1 and 3.1 x 10 10  1.9 x 10 9 kg for U2, with an average particles release height of ~13.5 km a.s.l.  0.5 km for both phases. Given the efficiency of the proposed approach and the plausibility of the


Introduction
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 to 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;Bonadonna and Costa, 2013;Volentik et al., 2010;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., Pyle, 1989;Carey and Sparks, 1986;Fierstein and Nathenson, 1992;Bonadonna et al., 2005;Bonadonna and Costa, 2012;Pyle, 2015;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 favor 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., 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 ( Constantinescu et al., 2022;White et al., 2017;Connor and Connor, 2006).Volcanological investigations often rely solely on the proximal (e.g., < 6 km from the vent) and very proximal deposits (e.g., < 4km 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., Constantinescu et al., 2022;White et al., 2017).
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., Cole et al., 2023;Horváth et al., 2022;Taylor et al., 2022;GVP, 2021a, b; Figure 1 in Supplementary Material).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).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 (< 4km 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 km and 6 km from the vent , with few measurements on the south-eastern 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 advectiondiffusion model used to calculate tephra mass accumulation per unit area (kg m -2 ), 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 grainsize 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
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 south-west of the volcano (Joseph et al., 2022).The eruption began on December 27, 2020 with a new lava dome forming inside the crater, located in the west south-west 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 April 9, 2021 and continued intermittently throughout the next ~ 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 Report, 2021b;Joseph et al. 2022).
Between April 10 and April 11, 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 Report, 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 April 12, 2021, intense explosions were associated with ash venting, creating plumes up to 12 km a.s.l and the generation of ash-rich PDC deposits.Discrete explosions continued to occur by the end of April 14, 2021 during which PDC deposits were emplaced in the valleys on the S and SW 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 April 15 -April 22, 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 favored tephra dispersal mainly east over the Atlantic Ocean.The low-level winds (0 -4 km a.s.l.) mostly blew to the west and south-west while the mid-level winds (8 -18 km a.s.l.) blew towards the east between 30 0 and 120 0 with a wind shear between 0 to 6 km a.s.l.Above 18 km the winds blew toward the west.Wind velocities for the low-level winds ranged between 0 -5 m s -1 while the mid-level winds reached up to 15 m s -1 .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 south-east, east and north-east with significant deposition on the windward (East) side of SV (Fig. 1 ).Low-level winds favored the dispersion of tephra to the north-west, west and south-west on the leeward (West) side with wind shear allowing the fine ash to reach as far south as Kingstown (Fig. 1 in Supplementary Material).Fine ash from the most vigorous plumes reached as far east as Barbados (> 150 km) with some trace ash reported in St. Lucia.
By April 22, 2021, tephra fallout from the 32 explosions produced significant damage to the northern side of the island.After the last explosion on April 22, 2021, the seismicity levels dropped and remained low, marking the end of the eruption (GVP Report, 2021c).

The tephra fall deposit
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 April 9 -10, 2021 and it reflects a north-east dispersal axis.The maximum measured thickness of this unit is ~ 20 cm close to the crater rim (Cole et al., 2023).On the windward coast the thickness ranges between 3 cm and 6.5 cm and on the leeward side between 0.2 cm and 5 cm (Figure 2 in Supplementary Material).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 April 10, 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 (~ 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 (Figure 2 in Supplementary Material).
Distinctive eruptive pulses associated with longer repose intervals between 18:30 UTC and 21:20 UTC on April 10, 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 millimeters 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 April 10, 2021 and 15:09 UTC on April 22, 2021 (Cole et al., 2023).U4 is an ash layer with accretionary lapilli (~8 mm) that thins from ~ 12 cm in the proximal facies to a few mm at > 4km away from the vent.An uneven lapilli layer, U5, reaches between 2 cm and 5 cm in the proximal locations on the south-western and south-eastern flanks.This lapilli layer is overlain by the fine-ash sequence of U6 which is found only on the south-western flanks of the volcano.U7 is identified only at a few localities on the south-western flanks of the volcano and consists of 3 to 4 very thin finegrained lapilli layers (Cole et al., 2023).

Field data collection and analysis
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 (~ 150 km East from the vent).The measured localities generally follow the shorelines on both coasts in 10 km to 15 km long transects and on the south-eastern flank (Fig. 1).The thickness of these units varies between 0.1 cm and 23 cm and are consistently thicker on the windward coast along the observed main axis of dispersion (Fig. 2 and Table 1 in Supplementary Material).Eleven samples were collected for grain size analysis (Fig. 1; Table 1 in Supplementary Material).Samples were dry-sieved down to 7 ϕ (31 μm) at 1 ϕ intervals and show a distribution between -3 ϕ and 7 ϕ.

Prior ESP estimates and the Inverse Modelling Approach
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., Connor et al., 2011;Bonadonna et al., 1998).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 highdimensional, ill-posed, non-linear inverse problems, including reconstructing tephra dispersion and sedimentation (e.g., Mingari et al., 2022a, b;White, 2018;Emerick and Reynolds, 2013;Crestani et al., 2013).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 modeled deposit.
We invert tephra mass per unit area (kg m -2 ) to estimate the erupted mass, average particle release height and total grainsize distribution.Mass load at each sampled location was obtained from thickness data using an average bulk deposit density of 1500 kg m -3 (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 km a.s.l. and 18 km a.s.l.(e.g., Horváth et al., 2022;GVP Report, 2021a).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 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 Table 2 of the Supplementary Material.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 ~70% and ~90% of the maximum plume tops (e.g., Horváth et al., 2022).The grainsize 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 grainsize 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 grainsize 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 modeled 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 (Constantinescu et al., 2022;White et al., 2017).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-behavioral 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 was downloaded from the NOAA NCEP/DOE Reanalysis 2 database (https://psl.noaa.gov/data/gridded/data.ncep.reanalysis2.html)(Kanamitsu et al., 2002).This data, as a daily mean, is 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.

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

Erupted mass
For the erupted mass of U1, the uncertainty in the posterior estimates was reduced by ~ 60% when compared to the prior distribution (Fig 3a; Table 1).A mean erupted mass was estimated at 3.3 x 10 10 kg (or a volume of 2.2 x 10 7 m 3 for 1500 kg m 3 bulk density) with an uncertainty range between 2.3 x 10 10 and 4.4 x 10 10 kg (Fig. 3a).The mean erupted mass for U2 was estimated at 3.1 x 10 10 kg (uncertainty range between 2.8 x 10 10 and 3.3 x 10 10 kg).The uncertainty in the posterior estimate was reduced by ~ 90% when compared to the prior (Fig. 3d).

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

TGSD mean
Figures 3c and 3f show 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 ϕ --0.15 ϕ) shows an uncertainty reduction of ~ 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 grainsize analysis.

Discussion and final remarks
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 uncertaint ies 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 (~ 150 km away from the vent).The major source of uncertainty in the case of the 2021 SVG eruption is the missing data from the medial and distal deposit facies.In general, this observational uncertainty along with deposit variability, contribute significantly to the estimates of erupted volume and plume height, regardless of the model (statistical or numerical) used for these determinations (e.g., Primerano et al., 2021;Connor et al., 2019;Bonadonna et al., 2015;Engwell et al., 2013).
The quality of the models' output relies on inversion of thickness (or mass load) and/or grainsize 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 grainsize information is available from across the deposit (e.g., Connor et al., 2019;White et al., 2017;Costa et al., 2016;Mele et al., 2020;Volentik et al., 2010;Bonasia et al., 2010;).
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., Constantinescu et al., 2022;White et al., 2017).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 x 10 9 -1 x 10 11 kg for U1 and 10 10 -10 11 for U2.The posterior uncertainty was reduced by ~60% (U1) and ~90% (U2), with a best-fit (i.e., mean of the posterior distributions) of 3.3 x 10 10 kg for U1 and 3.1 x 10 10 kg for Assuming a 1500 kg m -3 bulk density of the tephra deposit, these erupted masses indicate volumes of ~2.2 x 10 7 m 3 and ~2.1 x 10 7 m 3 for U1 and U2 respectively, comparable and within the uncertainty range of those obtained by Cole et al. (2023) using isopach interpolation (i.e., ~2.1 x 10 7 m 3 and ~2.4 x 10 7 m 3 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 x 10 10 kg (1.6 x 10 7 m 3 ) and the posterior uncertainty was reduced by ~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., Constantinescu et al., 2022;Connor et al., 2019;White et al., 2017;Magill et al., 2015;Connor and Connor, 2006) 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 wellconstrained 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)  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 km a.s.l. and 8 km a.s.l.(GVP Report, 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 ~ 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., Connor et al., 2019;Bonadonna et al., 2015;Bonadonna and Costa, 2012).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 grainsize data measured at each sampled location in addition to deposit thickness or mass (e.g., Connor et al., 2019;White et al., 2017;Volentik et al., 2010).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 grainsize 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 grainsize analysis of 11 samples taken from the island but this grainsize data was 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 grainsize classes at these locations.We acknowledge the fact that the grainsize 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.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 modeled deposit when different combinations between particle size -wind speed -release height are used (e.g., Constantinescu et al., 2022;White et al., 2017;Bonasia et al., 2010).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 is 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 (Fig. 3 in Supplemental Material).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 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 data set 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 x 10 10 kg and 3.5 x 10 10 kg), when compared with the results obtained by inverting only the field data (3.3 x 10 10 kg and 3.1 x 10 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 ://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 15, 2023 ://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 15, 2023 quantifying the ESPs uncertainty that inevitably results from missing deposit information and choice of forward model.
://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 15, 2023 ://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 15, 2023 from https://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 15, 2023 from https://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 15, 2023 ://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 15, 2023 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 ://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 15, 2023 ://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 15, 2023 ://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 15, 2023 ://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 15, 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.://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 15, 2023 ://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 15, 2023 ://www.lyellcollection.orgby University of Plymouth, Charles Seale-Hayne Library on Jun 15, 2023 may be further reduced if the wind field is parameterized and treated as unknown, and subsequently estimated through inversion along with the ESPs.

Figure Captions Figure 1 .
Figure Captions

Figure 2 .
Figure 2. Typical tephra sequence associated with the explosive phase of the LSV 2021 eruption consisting of several thin, alternating layers of ash and lapilli.The location of this stratigraphic section is represented by the green dot in Figure 1.

Figure 3 .
Figure 3. Prior and posterior kernel density plots for the erupted mass, average particle release height and TGSD mean for U1 and U2 units.The prior distributions are shown by the grey solid lines and the posterior estimates are shown by the solid blue lines.The uncertainty range is reported as 2 standard deviations and is shown by the vertical dashed lines colored accordingly for the prior and posterior distributions.

Figure 4 .
Figure 4.The estimated total grain size distribution for U1 and U2.

Figure 5 .
Figure 5. Isomass maps (units in kg m -2 ) for U1 and U2 as calculated by running the forward model with the best-fit set of ESPs obtained from the inversion analysis (UTM Zone 20).

Figure 6 .
Figure 6.Map showing the locations of the synthetic 0 kg measurements added to the original dataset.These locations (orange dots) are assumed to have received zero tephra during the eruption.This outline was based on available satellite imagery of the eruptive plumes.