Geochemistry, Geophysics,
Geosystems®

RESEARCH ARTICLE

10.1029/2025GC012498

**Key Points:**

* We adapt a 1D model of volcanic column to estimate the mass eruption rate from the evolution of the geometry of downwind plumes

* We perform a sensitivity analysis and validate the model using GOES-16 satellite images of the 2021 La Soufrière St Vincent eruption

* A new tool in the HOTVOLC system allows estimating source parameters in near real-time for Mt Etna eruptions

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

A. Michaud-Dubuy,
audrey.michaud-dubuy@uca.fr

**Citation:**

Michaud-Dubuy, A., & Gouhier, M. (2025). Validation of near real-time retrieval of plume mass eruption rates: The case of the 2021 eruption of La Soufrière, St Vincent. *Geochemistry, Geophysics, Geosystems, 26*, e2025GC012498. https://doi.org/10.1029/2025GC012498

Received 5 JUN 2025
Accepted 12 SEP 2025

# Validation of Near Real-Time Retrieval of Plume Mass Eruption Rates: The Case of the 2021 Eruption of La Soufrière, St Vincent

A. Michaud-Dubuy¹ and M. Gouhier¹

¹Université Clermont Auvergne, CNRS, IRD, OPGC, Laboratoire Magmas et Volcans, Clermont-Ferrand, France

<figure><img src="image_1.png" /></figure>

**Abstract** Satellite observations have enabled an important advance in the near-real-time quantification of the dynamic parameters of the volcanic plume spreading in the atmosphere. However, the link between these observations and the estimation of eruption source parameters, such as the mass eruption rate (MER), remains a scientific obstacle to be overcome. The previously developed methods to estimate the MER are less efficient for weak eruptions and/or occurring under strong wind conditions, which are the most frequent. Here, we update a 1-D volcanic column model for the estimation of the MER based on satellite measurements of wind-impacted plumes. The new model allows predicting the plume geometry as seen from space, and thus linking the source MER to the geometry far from the source. We find that the predictions mostly depend on the wind speed and the MER. We test the model using measurements made on GOES-16 images during the 2021 eruption of La Soufrière, St Vincent, and find a good agreement between our MER estimates and those found in the literature (with a mean ΔMER of -0.23). We finally test our ability in estimating the MER in near real-time using the HOTVOLC system and Meteosat-SEVIRI images of 10 paroxysms from Mt Etna. The GIS-based tool integrated to HOTVOLC allows easier measurements of the plume growth and will provide a robust tool for a rapid interpretation of satellite data in terms of source conditions, which are necessary inputs for tephra dispersion models, such as those used by the Volcanic Ash Advisory Centers.

**Plain Language Summary** Explosive volcanic eruptions produce ash columns that can rise in the atmosphere and be transported by winds hundreds of kilometers far from the volcano. As these clouds can have an impact on air traffic, health and crops, predictions of their paths in the atmosphere are crucial. The dispersal models used by monitoring centers to calculate their trajectories require several input parameters, such as the maximum column height and the flow rate feeding the eruption. While satellites now allow estimating the maximum column height, the flow rate remains more difficult to assess in real-time. Here, we adapted a volcanic column model to the calculation of the evolution of the cloud geometry even far from the volcano. The flow rates estimated from the comparison between our predictions and satellite measurements on clouds produced during the 2021 eruption of La Soufrière, St Vincent, compare well with those found in the literature. We finally demonstrate that our method could be used in near real-time thanks to the HOTVOLC platform, a satellite-based system developed at the Observatoire de Physique du Globe de Clermont-Ferrand (France) allowing simple and direct measurements on satellite images, and thus a rapid estimation of the flow rate during an eruption.

## 1. Introduction

Explosive volcanic eruptions stand as one of the most powerful natural phenomena on Earth. They are characterized by the formation of a volcanic column ejecting hot gas and pyroclasts at a large flow rate (>10⁴ kg s⁻¹). At the eruptive vent, this hot mixture of gas and tephra (i.e., volcanic particles of all sizes) has a bulk density greater than the ambient air and is therefore driven by its initial momentum only. As the column rises into the atmosphere, it experiences turbulent mixing with the cold atmospheric air, resulting in heating and rapid expansion of gas. If the turbulent mixing is efficient enough, the bulk density of the flow reduces to values lower than the atmospheric one, which allows the now positively buoyant volcanic column to rise up to several tens of kilometers in the atmosphere until it reaches a Neutral Buoyancy Height (NBH). On the other hand, if either the jet consumes its initial momentum before the buoyancy reversal or air entrainment is insufficient, the volcanic column collapses partially or totally and forms pyroclastic density currents.

The behavior of the buoyant volcanic column at the NBH depends on two main parameters: the eruption intensity (that can be represented by the mass eruption rate, MER in kg s⁻¹) and the wind velocity

© 2025 The Author(s). Geochemistry, Geophysics, Geosystems published by Wiley Periodicals LLC on behalf of American Geophysical Union. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schematic illustration of a weak/distorted volcanic plume, with both the volcanic column part framed in red and the downwind turbulent current framed in purple, showing some physical parameters used in this study (see text and Notation section for symbol description). The inset shows the top view of the plume, as seen on satellite images.</figcaption></figure>

(Bonadonna & Phillips, 2003). Theoretical and experimental studies indeed shown that for high eruption intensity and/or low wind speed, the volcanic column forms a “strong” plume largely unaffected by winds (Bonadonna & Phillips, 2003; Carazzo et al., 2014). In this case, the plume can overshoot the NBH and reach a greater maximum height, before collapsing back to the NBH and forming a circular turbulent current spreading radially. On the other hand, for low eruption intensity and/or high wind speed, the volcanic column bends over and forms a “weak” plume strongly controlled by winds (Bursik, 2001). In this case, the plume forms a lateral turbulent current spreading downwind around the NBH (Figure 1). Experimental studies also demonstrated the existence of a third, transitional regime called “distorted” plume, in which the plume slightly overshoot the NBH before spreading downwind (Carazzo et al., 2014). The strong/weak plume transition is thus governed by the ratio of wind velocity to the characteristic velocity scale of the plume, classically represented by a dimensionless number (Carazzo et al., 2014; Degruyter & Bonadonna, 2013; Devenish et al., 2010). Note that in this work, we use the word “column” for the subvertical/bended convective region from the eruptive source to the maximum height, “downwind turbulent current” refers to the part that spreads laterally around the NBH, while “plume” is a collective term used to describe the overall phenomenon (column + turbulent current, Figure 1).

Regardless of the plume regime, volcanic eruptions can spread various products (ash, SO₂, halogens,...) up to hundreds of kilometers from the volcano and lead to various types of damage, from air pollution (Stewart et al., 2022) to major air traffic disruption (Picquout et al., 2013), as well as global climate change (Sigl et al., 2015). The need to understand the dynamics of volcanic columns, the conditions leading to the formation of stable (or collapsing) plumes, and to accurately predict tephra dispersal extent led the volcanological community to develop increasingly sophisticated numerical models (Costa et al., 2016). 1-D models of volcanic columns (based on the work of Wilson (1976) and Woods (1988)) only describe the evolution of the plume from the vent to the maximum height. 3-D models are extremely useful to understand more precisely the physical processes or to incorporate rapid changes in wind direction and velocity, but their use comes with high computational costs and difficulties in using these models in real-time during a volcanic crisis. An intercomparison exercise performed by Costa et al. (2016) has shown reasonably consistent predictions between 1D and 3D models. 1D models were also proven to be capable of resolving the general behavior of the weak plume regime (Costa et al., 2016).

For operational purposes, the Volcanic Ash Advisory Centers (VAACs, in charge of providing expertise on volcanic ash for civil aviation) use volcanic ash transport and dispersion models requiring several Eruption Source Parameters (ESPs) to run. These ESPs include plume height, total erupted mass, MER, duration, and total grain-size distribution (TGSD). One of the current challenges in volcanology is thus to be able to provide
accurate ESPs as rapidly and precisely as possible during an eruption. As volcanic plumes are now routinely tracked in real time by remote sensing methods, satellite imagery has proven to be an effective way to provide ESPs quickly during an eruption. As volcanic plumes can reach the tropopause (9–17 km a.s.l.) in a few minutes, the high temporal resolution of Geostationary Earth Orbit platforms (such as GOES-16, METEOSAT, Himawari-9) is a valuable tool for the rapid detection of volcanic ash in the atmosphere, and the estimation of ash concentrations and maximum column height (Chazette et al., 2012; Gouhier et al., 2019, 2020; A. Prata, 1989; A. Prata & Prata, 2012). As an example, the HOTVOLC monitoring system developed at the Observatoire de Physique du Globe in Clermont-Ferrand (OPGC) uses METEOSAT satellite data received and processed in near real-time to detect volcanic ash for 50 volcanic targets, and provides the Toulouse VAAC (Météo-France) with ESPs (such as the maximum column height) so they can run their dispersal model MOCAGE-Accident (Gouhier et al., 2020).

Operational protocols strongly rely on the observed maximum column height. The latter is usually estimated from the Cloud-Top Temperature method, allowing to convert satellite-based measurements of the temperature at the top of the cloud to column height by comparison with atmospheric temperature profiles. The MER, a difficult parameter to measure in real-time, is in turn estimated from the maximum column height using scaling laws (OD models) based on observed eruptions, such as Mastin et al. (2009); Carazzo et al. (2014). Overall, this method relies on several assumptions and thus bears strong uncertainties (Woods & Self, 1992). Indeed, the cloud must be optically thick so that the cloud surface temperature is correctly estimated; the method is only applicable to tropospheric plumes since the atmospheric temperature profile above the tropopause is not monotonic (making impossible to obtain a single solution of maximum height); and the scaling law usually used to link the maximum height to the MER has its own uncertainties (reflecting those of field observations and simplified descriptions of the atmosphere on which the OD model is based) (Costa et al., 2016; Mastin, 2014; Mastin et al., 2009). The main goal of this study is to propose a robust alternative method to retrieve the MER that will be added to the already high value-added open-access products delivered by HOTVOLC, and thus improve the ash dispersal simulations performed by the VAAC.

Current efforts consist in coupling plume models and satellite observations to be able to infer other ESPs, such as the MER. Pouget et al. (2013) developed a method to estimate the source MER from the temporal evolution of the plume measured on satellite images. The longitudinal and lateral spreading of the downwind turbulent current indeed give major constraints on the wind speed and direction at the NBH and on the MER fueling the volcanic column, respectively. This model however does not take into account the column dynamics from the eruptive source. Ishii et al. (2022) recently developed the 1-D model of volcanic column NIKS-1D to estimate rapidly after an eruption the particle size distribution released and nourish tephra dispersal models. The model incorporates a 1-D model of volcanic column considering particle sedimentation and a downwind gravity current model around the NBH validated for weak plumes. Inspired by these works, we aim at updating the 1-D model of volcanic column PPM (Girault et al., 2014, 2016; Michaud-Dubuy et al., 2018, 2020) to model the evolution of the plume geometry from the source to the far field. The model is then used to infer the MER from the satellite observations. PPM is a good candidate for this study as it includes the combined effects of the reduced turbulent entrainment at the base of the eruptive column, the gas entrapment in large pumices during the fragmentation process, open porosity, and wind, key processes determining the plume regime (Aubry et al., 2017; Michaud-Dubuy et al., 2018). PPM also takes into account the sedimentation of particles and various types of grain-size distributions (that can have an impact on the maximum column height, Girault et al. (2016)), while having a low computational cost. We will focus in this study on weak and distorted plumes, as they occur more often than strong plumes. The recent eruption of La Soufrière, St Vincent in 2021 offered a unique occasion to test the ability of our model to estimate source MER, as the eruption was followed at close range and was well visible on GOES-16 images. Significant work has been devoted to investigating the dynamics of the eruption using variety of techniques (Robertson et al., 2023). This led to the availability of extensive data on ESPs (such as the MER) for multiple explosions observed during the explosive phase of the eruption in April 2021.

We first present the 2021 eruption of La Soufrière, St Vincent and the various data available for this study. We then introduce our new version of the PPM model linking the volcanic column dynamics to the turbulent current region. The results section presents a sensitivity analysis of the model, the application to La Soufrière, St Vincent eruption, and a test of the HOTVOLC system, using several paroxysms of Mt Etna as case studies. Finally, we discuss the limitations of our modelling approach, and the improvements compared to previous estimates.
## 2. The 2021 Eruption of La Soufrière, St Vincent, West Indies

St Vincent is a volcanic island located in the southern part of the Eastern Caribbean Island arc. La Soufrière is currently the only active volcano on the island, and one of the most active in the arc with five magmatic eruptions since 1718 (Robertson, 2005; Robertson et al., 2023). These eruptions include extrusions of lava domes (1784, 1971–1972, 1979, 2020–2021), as well as explosive sequences producing volcanic columns that can endure partial/total collapse and form pyroclastic density currents (1718, 1812, 1902–1903, 1979, 2020–2021; Robertson et al. (2023), Cole et al. (2023)). All explosive eruptions have generated extensive fallout that sometimes reached Barbados located 170 km to the east (Shepherd et al., 1979). After a 41 years hiatus, La Soufrière, St Vincent resumed its eruptive activity on 27 December 2020 by the extrusion of a basaltic andesite lava dome. This effusive activity followed a minimal period of precursory seismicity (Joseph et al., 2022; Robertson et al., 2023) and lasted for about 3 months. The new lava dome grew in the summit crater, next to the 1979 dome (Robertson et al., 2023), and its total volume was estimated at $16–19 \times 10^6 \text{ m}^3$ (Stinton, 2023). Following volcano tectonic earthquake swarms on 23 March and 5 April, and banded tremor on 8 April, explosive activity started on 9 April 2021 with a single explosion at 12:41 UTC (Joseph et al., 2022). By that time, evacuation orders had already been given by the Prime Minister (Robertson et al., 2023). The explosive eruptive sequence lasted from 9 to 23 April 2021, at which date the eruption was declared over.

The eruption was well-observed and monitored both at close quarters and remotely thanks to satellite and radar imagery. Various methods based on field deposits (Cole et al., 2023), seismicity (Sparks et al., 2023) or tephra dispersal modeling (Constantinescu et al., 2023), for example, were used to estimate ESPs, making this eruption a valuable case-study. A thorough description of the eruption and the numerous data and studies associated are given in Robertson et al. (2023). Based on Real-time Seismic Amplitude Measurement (RSAM) and thermal infrared imagery from GOES-16 geostationary satellite, Sparks et al. (2023) identified at least 32 discrete explosions between 9 and 22 April 2021. The measured brightness temperature at the top of the cloud was then used to approximate the maximum cloud heights of each event by comparison with atmospheric temperature profiles (Cloud-Top Temperature method). These maximum heights were in turn used to constrain peak intensities (in $\text{m}^3 \text{s}^{-1}$) from an empirical power-law relationship derived from the analysis of global data sets and depending on a scaling constant and an empirical exponent that can be adjusted (Mastin et al., 2009; Sparks et al., 1997, 2023). The analysis of the evolution of peak intensities, magma supply rates, and erupted volumes led the authors to divide the eruption into four stages (Sparks et al., 2023). Table 1 compiles the ID of each explosive event identified, their dates, the corresponding eruption stage, the time (UTC) at which the volcanic cloud was visible on satellite imagery, the plume maximum heights and peak intensities (converted into MER) estimated by Sparks et al. (2023).

The first stage, going from 9 April 12:50 UTC (the timing of eruptions are given in this study by the time at which the volcanic cloud was visible on satellite images) to 10 April 07:50 UTC, includes 14 explosions for which a peak intensity was estimated by Sparks et al. (2023) (see their Table 3). From 9 April at 19:10 to the end of stage 1, there was a period of continuous pulsating activity making it difficult to distinguish the explosions from one another based solely on satellite imagery, especially for the events 4 to PO16 (i.e., from 10 April 00:20 UTC to 04:30 UTC). For this reason, we regrouped these events in our Table 1 and calculated a mean column height and MER from the data of Sparks et al. (2023) (using their Tables 2 and 3. Stage 1 corresponds to the first lapilli-rich deposit unit, named U1 by Cole et al. (2023), whose volume was estimated at $10–13 \times 10^6 \text{ m}^3$ (Cole et al., 2023; Constantinescu et al., 2023; Sparks et al., 2023). The second stage, going from 10 April at 09:50 UTC to 11 April at 08:10 UTC comprises 12 more powerful (with higher MER) and discrete explosive events every ~1–2.5 hr. Four deposit units, alternating from ash-rich (U2, U4) to coarse-grained (U3, U5) layers, were attributed to this stage (Cole et al., 2023). This alternation is consistent with the instability of the volcanic column, starting on 10 April at 18:30 UTC with the production of pyroclastic density currents, that lasted almost until the end of the eruption (Cole et al., 2023). The tephra fall volumes of U2 and U3 were estimated using different methods at $9.6–14.5 \times 10^6 \text{ m}^3$ for U2 (Cole et al., 2023; Constantinescu et al., 2023; Sparks et al., 2023), and $3.8–4.7 \times 10^6 \text{ m}^3$ for U3 (Cole et al., 2023; Sparks et al., 2023). The 11 events of the third stage, from 11 April at 08:10 UTC to 14 April at 15:40 UTC, show decreasing intensity and longer repose periods between
<table><caption>**Table 1**
Timing and Eruption Source Parameters of Explosive Events Between 9 and 14 April 2021 at La Soufrière Volcano, St Vincent, Data From Sparks et al. (2023)</caption><thead><tr><th>Eruption stage</th><th>Event ID</th><th>Day, April 2021</th><th>Cloud time (UTC)</th><th>Plume max height (km)</th><th>MER (kg s<sup>-1</sup>)</th><th>Standard deviation</th></tr></thead><tbody><tr><td rowspan="10">1</td><td>1</td><td>09</td><td>12:50</td><td>13.0</td><td>4.08 E + 06</td><td>1.54 E + 06</td></tr><tr><td>2</td><td>09</td><td>19:10</td><td>15.5</td><td>8.60 E + 06</td><td>3.35 E + 06</td></tr><tr><td>PO1</td><td>09</td><td>20:20</td><td>15.5</td><td>8.60 E + 06</td><td>3.35 E + 06</td></tr><tr><td>PO2</td><td>09</td><td>21:30</td><td>14.5</td><td>6.50 E + 06</td><td>2.50 E + 06</td></tr><tr><td>3</td><td>09</td><td>22:40</td><td>14.3</td><td>6.13 E + 06</td><td>2.35 E + 06</td></tr><tr><td>PO7</td><td>09</td><td>23:40</td><td>15.7</td><td>9.08 E + 06</td><td>3.55 E + 06</td></tr><tr><td>4-PO16</td><td>10</td><td>00:20</td><td>15.3–15.9</td><td>9.00 E + 06</td><td>3.52 E + 06</td></tr><tr><td>5</td><td>10</td><td>06:50</td><td>10.0</td><td>1.04 E + 07</td><td>4.08 E + 06</td></tr><tr><td>6</td><td>10</td><td>07:50</td><td>14.0</td><td>5.63 E + 06</td><td>2.15 E + 06</td></tr><tr><td rowspan="11">2</td><td>7</td><td>10</td><td>09:50</td><td>16.0</td><td>9.83 E + 06</td><td>3.85 E + 06</td></tr><tr><td>8</td><td>10</td><td>11:00</td><td>16.0</td><td>9.83 E + 06</td><td>3.85 E + 06</td></tr><tr><td>9</td><td>10</td><td>12:10</td><td>16.0</td><td>9.83 E + 06</td><td>3.85 E + 06</td></tr><tr><td>10</td><td>10</td><td>13:00</td><td>16.0</td><td>9.83 E + 06</td><td>3.85 E + 06</td></tr><tr><td>11</td><td>10</td><td>14:40</td><td>16.0</td><td>9.83 E + 06</td><td>3.85 E + 06</td></tr><tr><td>12</td><td>10</td><td>16:30</td><td>16.0</td><td>9.83 E + 06</td><td>3.85 E + 06</td></tr><tr><td>13</td><td>10</td><td>18:50</td><td>15.2</td><td>7.93 E + 06</td><td>3.08 E + 06</td></tr><tr><td>14</td><td>10</td><td>21:30</td><td>16.0</td><td>9.83 E + 06</td><td>3.85 E + 06</td></tr><tr><td>15</td><td>10</td><td>23:10</td><td>16.2</td><td>1.04 E + 07</td><td>4.05 E + 06</td></tr><tr><td>16</td><td>11</td><td>01:00</td><td>16.1</td><td>1.01 E + 07</td><td>3.95 E + 06</td></tr><tr><td>17</td><td>11</td><td>02:50</td><td>16.1</td><td>1.01 E + 07</td><td>3.95 E + 06</td></tr><tr><td rowspan="11">3</td><td>18</td><td>11</td><td>05:10</td><td>15.7</td><td>9.08 E + 06</td><td>3.55 E + 06</td></tr><tr><td>19</td><td>11</td><td>08:10</td><td>14.8</td><td>7.08 E + 06</td><td>2.73 E + 06</td></tr><tr><td>20</td><td>11</td><td>10:50</td><td>14.9</td><td>7.30 E + 06</td><td>2.83 E + 06</td></tr><tr><td>21</td><td>11</td><td>13:40</td><td>14.5</td><td>6.50 E + 06</td><td>2.50 E + 06</td></tr><tr><td>22</td><td>11</td><td>18:20</td><td>13.6</td><td>4.95 E + 06</td><td>1.89 E + 06</td></tr><tr><td>23</td><td>11</td><td>20:20</td><td>14.0</td><td>5.63 E + 06</td><td>2.16 E + 06</td></tr><tr><td>24</td><td>12</td><td>00:50</td><td>14.0</td><td>5.63 E + 06</td><td>2.16 E + 06</td></tr><tr><td>25</td><td>12</td><td>08:10</td><td>14.3</td><td>6.13 E + 06</td><td>2.36 E + 06</td></tr><tr><td>26</td><td>12</td><td>21:30</td><td>12.2</td><td>3.13 E + 06</td><td>1.18 E + 06</td></tr><tr><td>27</td><td>13</td><td>10:40</td><td>14.0</td><td>5.63 E + 06</td><td>2.16 E + 06</td></tr><tr><td>28</td><td>14</td><td>02:40</td><td>13.5</td><td>4.80 E + 06</td><td>1.82 E + 06</td></tr><tr><td>29</td><td>14</td><td>15:40</td><td>13.0</td><td>4.08 E + 06</td><td>1.54 E + 06</td></tr></tbody></table>

Note. Peak intensities in m³ s⁻¹ were converted into mass eruption rates using a magma density of 2500 kg m⁻³.

explosions. The deposit units U5, U6 and U7 are attributed to this stage, but with larger uncertainty due to their thinness in the field (especially for U7) (Cole et al., 2023). The volumes of units 4 to 7 were estimated between 7.1–15.1 × 10⁶ m³ (Cole et al., 2023; Sparks et al., 2023). The fourth stage comprises the last three and much weaker explosions that occurred until 22 April, for which the estimated MER are very low (≈ 7 × 10⁴ kg s⁻¹) and subjected to higher uncertainties (Sparks et al., 2023). No tephra fall deposit could be associated with this stage in the field (Cole et al., 2023). In this study, we seek to estimate the source MER of each event presented in Table 1, thanks to a combination of downwind plume measurements on GOES-16 satellite images and theoretical predictions of the downwind plume geometry from an updated version of the 1-D model of volcanic column PPM, and compare them to the intensities estimated by Sparks et al. (2023).
3. Methods

3.1. The Volcanic Column Model PPM

3.1.1. Conservation Equations and Constitutive Laws

We build the new turbulent current model over the Paris Plume Model (PPM) of Michaud-Dubuy et al. (2020) that relies on a 1D steady-state “top-hat” formalism (i.e., assuming a conical jet in which all the dynamical variables are considered constant inside the jet and zero outside; Morton et al. (1956)) in a plume-centered coordinate system, where s denotes the curvilinear abscissa along the plume centerline and θ is the local inclination of the jet (Hewett et al., 1971; Hoult et al., 1969). For a particle-laden volcanic jet and considering the effect of wind, the conservation laws of mass, axial and radial momentum, and energy fluxes are written as (Bursik, 2001; Girault et al., 2016; Michaud-Dubuy et al., 2020):

$$
\frac{d}{ds}(\rho u r^2) = 2\rho_a r u_e + \sum_{\phi=1}^{N_{\phi}} \frac{dQ_{\phi}}{ds}, \qquad (1)
$$

$$
\frac{d}{ds}(\rho u^2 r^2) = (\rho_a - \rho) gr^2 \sin \theta + w \cos \theta \frac{d}{ds}(\rho ur^2) + u \sum_{\phi=1}^{N_{\phi}} \frac{dQ_{\phi}}{ds}, \quad (2)
$$

$$
(\rho u^2 r^2) \frac{d\theta}{ds} = (\rho_a - \rho) gr^2 \cos \theta - w \sin \theta \frac{d}{ds}(\rho ur^2), \quad (3)
$$

$$
\frac{d}{ds}(\rho c_p T u r^2) = 2\rho_a c_a T_a r u_e - \rho_a g u r^2 \sin \theta + c_p T \sum_{\Phi=1}^{N_{\Phi}} \frac{dQ_{\Phi}}{ds}, \quad (4)
$$

where *r* is the column radius, *u* is the average vertical velocity along the column, *w* is the wind speed, *g* is the acceleration of gravity, *ρ*, *c*<sub>*p*</sub> and *T* are the density, specific heat, and temperature at constant pressure of the bulk mixture, respectively, *ρ*<sub>*a*</sub>, *c*<sub>*a*</sub>, and *T*<sub>*a*</sub> are those of the atmosphere, *Q*<sub>*φ*</sub> = *x*<sub>*φ*</sub>*ρ* *u*<sup>2</sup>, where *x*<sub>*φ*</sub> is the mass proportion of *φ*-sized particles in the grain-size distribution, is the mass discharge rate of *φ*-sized particles (kg s<sup>-1</sup>) distributed within 20 classes of grain sizes (*N*<sub>*Φ*</sub>) ranging from *d*<sub>min</sub> = 10*φ* (1 µm) to *d*<sub>max</sub> = -9*φ* (0.5 m) with one *φ* intervals, and *u*<sub>*e*</sub> is the rate of turbulent entrainment of ambient air into the column.

The calculation of u_e is as follows (Hewett et al., 1971):

$$
u_e = \alpha|u - w \cos(\theta)| + \beta|w \sin(\theta)|, \quad (5)
$$

where β is the coefficient quantifying additional entrainment related to the presence of wind (Aubry & Jelli-
nek, 2018; Devenish et al., 2010; Mastin, 2014), and α is the “top-hat” entrainment coefficient of Morton
et al. (1956). Here, we take β = 0.1, a value consistent with the best range determined by Suzuki and Koya-
guchi (2015) and used by Ishii et al. (2022). As in Kaminski et al. (2005), we consider that α depends on the
buoyancy of the column relative to the ambient air, and can be expressed by:

$$
\alpha = 0.0675 + \left(1 - \frac{1}{A}\right) \mathrm{Ri}, \qquad (6)
$$

where Ri = g(ρ_a - ρ) r/ρ_a u² is the local Richardson number and A is a dimensionless parameter that depends on the flow structure and evolves as a function of the downstream distance from the source and buoyancy (see Appendix A of Michaud-Dubuy et al. (2018) for more details).

The set of conservation equations is complemented by constitutive equations giving the evolution of the ther-
modynamic properties of the column with altitude, assuming no thermal and mechanical disequilibrium between
the volcanic gas and the particles (Woods, 1988):
$$ \frac{1}{\rho} = \frac{(1 - x_g)}{\rho_p} + \frac{x_g R_g T}{P_a}, \qquad (7) $$

$$ x_g = 1 + (x_{g0} - 1) \frac{\rho_0 u_0 r_0^2}{\rho u r^2}, \qquad (8) $$

$$ R_g = R_a + (R_{g0} - R_a) \frac{(1 - x_g)}{x_g} \frac{x_{g0}}{(1 - x_{g0})}, \qquad (9) $$

$$ c_p = c_a + (c_{p0} - c_a) \left( \frac{1 - x_g}{1 - x_{g0}} \right), \qquad (10) $$

where $\rho_p = 2000 \text{ kg m}^{-3}$ is the average particle density, $x_g(z)$ is the effective gas mass fraction in the flow (i.e., the total fraction of gas minus the amount of gas trapped in bubbles in magmatic fragments and considering open porosity), $P_a(z)$ is the atmospheric pressure, $R_{g0} = 461 \text{ J kg}^{-1} \text{ K}^{-1}$ and $R_a = 287 \text{ J kg}^{-1} \text{ K}^{-1}$ are the bulk column and the air gas constant, respectively. The subscript 0 denotes a value at the vent. $P_a(z)$, $T_a(z)$, and $\rho_a(z)$ depend on the atmospheric conditions which we take as typical mid-latitude or tropical (Carazzo et al., 2008).

Furthermore, as we focus solely on plumes spreading downwind (i.e., whose dynamics are influenced by wind) in this study, we introduced in our model the dimensionless ratio proposed by Carazzo et al. (2014), based on laboratory experiments, and fully describing the transition from weak to distorted plumes, and from distorted to strong plumes:

$$ \frac{w^*}{u^*} = \frac{\alpha_p^{1/2} w}{(F_0 N)^{1/4}}, \qquad (11) $$

where $u^*$ is the plume velocity ratio, $w^*$ the wind velocity ratio, $\alpha_p = 0.125$ is the “top-hat” entrainment coefficient for a self-similar pure plume (Carazzo et al., 2006; Wang & Law, 2002), $F_0$ is the source buoyancy flux, and $N$ is the Brunt-Väisälä frequency defined as $N^2 = -\frac{g}{\rho} \frac{d\rho}{dz}$. Thus, a $w^*/u^* < 0.11$ characterizes a strong plume, whereas for $w^*/u^* > 0.37$ a weak plume develops. For intermediate conditions (i.e., $0.11 \le w^*/u^* \le 0.37$), the plume is distorted (Carazzo et al., 2014). This dimensionless ratio was calculated at each run to insure that we based our conclusions on weak and distorted plumes only.

### 3.1.2. Particle Sedimentation

We consider the mass loss of $\phi$-sized particles from the column edges to be proportional to the mass discharge rate of particles $Q_\phi$ and to the terminal fall velocity $V_\phi$, such as (Michaud-Dubuy et al., 2018; Woods & Bursik, 1991):

$$ \frac{dQ_{\phi}}{ds} = -p_s \frac{Q_{\phi}}{r} \frac{V_{\phi}}{u}, \qquad (12) $$

where $p_s = 0.27$ is a probability of sedimentation (Ernst et al., 1996) and $V_\phi$ is the terminal fall velocity. For a given particle size, $V_\phi$ is calculated using (Bonadonna et al., 1998):

$$ V_{\phi} = \begin{cases} \left[ \frac{3.1 g d_{\phi} \Delta \rho}{\rho_a} \right]^{1/2} & \text{for } \text{Re}_p \ge 500, \\ d_{\phi} \left[ \frac{4 g^2 \Delta \rho^2}{225 \mu_{\text{air}} \rho_a} \right]^{1/3} & \text{for } 0.4 \le \text{Re}_p \le 500, \\ \left[ \frac{d_{\phi}^2 g \Delta \rho}{18 \mu_{\text{air}}} \right] & \text{for } \text{Re}_p \le 0.4, \end{cases} \qquad (13) $$

where $\Delta\rho = (\rho_p - \rho)$, $\text{Re}_p = \rho d_\phi V_\phi / \mu_{\text{air}}$ is the particle Reynolds number, $d_\phi$ is the particle diameter of a $\phi$-sized particle, and $\mu_{\text{air}}$ is the dynamic viscosity of air.
### 3.1.3. Exit Conditions at the Vent

In this study, we consider mass eruption rates ranging from $10^3$ to $10^9$ kg s$^{-1}$ consistent with the range for historical explosive eruptions, and the total mass fraction of gas is set to $x_{tot} = 4$ wt% in all calculations. Field data and fragmentation experiments have shown that volcanic rocks fragment according to a power-law distribution:

$$N(R_p \ge r_p) = \lambda r_p^{-D}, \quad (14)$$

where $N(R_p \ge r_p)$ is the number of fragments of size $R_p$ larger than or equal to $r_p$, $\lambda$ is a normalization constant and $D$ is the power-law exponent. For a given set of minimum and maximum grain sizes, the exponent $D$ fully characterizes the power-law distribution (Girault et al., 2016; Kaminski & Jaupart, 1998). The effective amount of gas $x_{g0}$ in the turbulent jet at the vent is calculated by considering the amount of gas exsolved, trapped in large particles and released thanks to open porosity. For a given value of the power-law exponent $D$, $x_{g0}$ then corresponds to a proportion of $x_{tot}$ (see Figure 2b in Michaud-Dubuy et al. (2018)). Here we use in all simulations an open porosity of 65%, a value consistent with open porosity measurements made on pumices from explosive eruptions of various magma composition (Michaud-Dubuy et al., 2018). Three values of the power-law exponent $D$ are tested during the validation process: $D = 2.8$ (i.e., coarse TGSD), $D = 3.0$ (mixed TGSD), and $D = 3.3$ (fine TGSD). Following Michaud-Dubuy et al. (2018) and considering $x_{tot} = 4$ wt%, such values correspond to $x_{g0} = 2.96$ wt% (for $D = 2.8$), $x_{g0} = 3.38$ wt% (for $D = 3.0$), and $x_{g0} = 3.88$ wt% (for $D = 3.3$). For examples on how the value of $D$ can be interpreted in terms of TGSD, we refer the reader to Figure 2 in Girault et al. (2014).

The exit velocity $u_0$ is calculated assuming the hypothesis of a freely decompressing jet (Woods & Bower, 1995):

$$u_0 = 1.8 \sqrt{x_{g0} R_{g0} T_0}, \quad (15)$$

where $T_0 = 1200$ K is taken as the average temperature of andesitic magma.

### 3.1.4. Wind Conditions

The atmospheric conditions are taken in this study as typical mid-latitude or tropical, using the same parameterization as in Carazzo et al. (2008). The tropopause height is set at 11 km for the mid-latitude atmosphere, and 17 km for the tropical one. Our model also considers three types of wind speed profile as a function of the altitude (Girault et al., 2016; Michaud-Dubuy et al., 2020). In the “uniform” wind profile, we set a constant crosswind speed from the volcanic vent up to the maximum height of the column. The “linear” wind profile considers a linear increase of wind speed from no wind at $z = 0$ to a maximum wind speed (set by the user) at the tropopause height. The wind speed then remains constant above the tropopause and up to the maximum column height. The “complex” wind profile is a parabolic profile, taken as a first-order model for typical wind conditions in the atmosphere (Hedin et al., 1996). The wind speed is set at $5 \text{ m s}^{-1}$ at $z = 0$, increases up to a maximum value at the tropopause (defined by the user), then decreases above the tropopause until it reaches 2/3rd of its maximum value, and finally keeps this constant speed at higher altitude. For each type of wind profile, we test three maximum wind speeds (10, 20, and 30 m s$^{-1}$).

## 3.2. Downwind Turbulent Current Model

In the current version of the model, PPM allows calculating the bending of a plume subjected to winds in the atmosphere, and thus the turbulent current growth (i.e., the evolution of its width, $L$) as a function of the downwind distance $x$, but only until the plume reaches its maximum height, at which point the calculation stops. Even in case of strong winds (and thus strong bending), the volcanic column reaches its maximum height at a small horizontal distance from the vent (about 1 km), making the comparison with plume growth measurements on satellite images relatively limited. We thus seek to upgrade the PPM model by incorporating the turbulent current spreading dynamics in order to predict the plume growth rate far enough from the source and allow a direct comparison with satellite data.

After reaching the NBH, the plume is laterally advected by wind while also spreading orthogonally. In the literature, this behavior is referred as a downwind turbulent gravity current, meaning that the plume spreads as a
turbulent gravity current due to the difference in density between the upper and lower parts of the clouds and the stratified atmosphere at the considered altitude, while being advected in the downwind direction (Woods & Kienle, 1994). Classically, the characteristic turbulent current growth velocity can be written (Bursik et al., 1992; Sparks et al., 1997; Woods & Kienle, 1994):

$$U = \lambda_u NH, \qquad (16)$$

where $\lambda_u$ is the spreading constant of order unity depending on a flow geometry (Bursik et al., 1992), $N = 0.01$ s⁻¹ is the Brunt-Väisälä frequency for a typical atmosphere, and $H$ is the turbulent current thickness. This velocity can also be defined as the variation of the turbulent current width $L$ with time (the crosswind growth direction being the only one not affected by wind speed). In this case,

$$\frac{dL}{dt} = \lambda_u NH \qquad (17)$$

Injecting the volumetric flux of the turbulent current $V_1 = wLH$ into Equation 17 gives:

$$\frac{dL}{ds} = \frac{\lambda_u NV_1}{w^2 L} \qquad (18)$$

The analytical solution of Equation 18 can be written (Ishii et al., 2022):

$$L(s) = \frac{\sqrt{2\lambda_u NV_1(s - s_0) + w^2 L_0^2}}{w}, \qquad (19)$$

where $s_0$ marks the beginning of the downwind turbulent current (somewhere in the downwind direction from the eruptive vent). The source conditions of this downwind turbulent current model, that is, its volumetric flux $V_1$, its width $L_0$ at the beginning of its spreading, the wind speed $w$ at the NBH where the turbulent current is spreading, and $s_0$, are provided by PPM at the NBH. Following the approach of Ishii et al. (2022), we use an empirical parameter $\mu$ to link the plume volumetric flux $V_0$ at the NBH with the turbulent current volumetric flux $V_1$:

$$\mu = \frac{V_1}{V_0}, \qquad (20)$$

Because the entrainment of ambient air is greater at the top of the eruptive column than at its base, Ishii et al. (2022) suggest that $V_1$ would be substantially greater than $V_0$, and thus that $\mu \ge 1$. As we focus in this study on weak and distorted plumes with no (or very limited) overshoot region at the top of the column, we assume that the volumetric flux remains constant at the transition column/turbulent current. In this framework, we set the empirical parameter $\mu = 1$ for most of the calculations, and we test the effect of varying $\mu$ from 1 to 4 (Section 4.2) that is discussed in Section 5.1.

### 3.3. Measurements on Satellite Data (GOES-16 and MSG)

To test our predictions on the evolution of the plume geometry, we measured on satellite images the downwind plumes observed for 33 explosive events during the 2021 La Soufrière, St Vincent eruption: the 32 events identified by Sparks et al. (2023) between 9 and 14 April (see Table 1), and another event on 10 April at 05:20 UTC, not detected on the RSAM data of Sparks et al. (2023) but clearly visible on satellite images. Data were from GOES-16 satellite imagery combining the infrared anomaly with ash detected in visible light (“ABI RadF ash”), available as a supplementary file in Cole et al. (2023) (https://doi.org/10.6084/m9.fishshare.c.6474317). Using this video as an image sequence in the ImageJ software (Schneider et al., 2012), we measured at each time step (i.e., every 10 minutes) both the plume width $L$ and the corresponding downwind distance $x$ (i.e., the distance from the eruptive vent–La Soufrière crater–to the point where the plume width measurement was made). As an example, we show in Figure 2a the five plume growth measurements made on each image for event 6 (10 April 2021) between 07:50 and 08:30 UTC. We also estimated the average wind speed during each explosion from the measurements of the downwind distance $x$ as a function of time.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** (a) Example of plume growth measurements made on a GOES-16 satellite image showing the event 6 (10 April 2021, 07:50 UTC) of La Soufrière, St Vincent eruption. (b) Screenshot of the HOTVOLC interface showing the eruptive plume produced by Mt Etna during the 21 February 2022 eruption, and the plume growth measurements made on the satellite images (MSG-SEVIRI) using the new GIS-based tool (the toolbar on the left can be accessed by clicking the set square at the top right of the interface).</figcaption></figure>

Because of its strong operational potential, we wish to integrate in a near future the results of this study into the satellite-based system HOTVOLC developed at the OPGC. HOTVOLC (https://hotvolc.opgc.fr) is a Web-GIS volcano monitoring system developed at the OPGC in 2009 and using SEVIRI (Spinning Enhanced Visible and Infrared Imager) sensor on-board METEOSAT geostationary satellite (Gouhier et al., 2020). The SEVIRI sensor allows the characterization in real-time of volcanic ash, sulfur dioxide, and lava flow emissions for about 50 active volcanoes, with a frequency of 1 image/15 min and a pixel resolution of 3 km × 3 km. During an explosive eruption, HOTVOLC provides data such as the plume altitude to the Toulouse VAAC, which is required to promptly run their dispersal model and provide a first map showing the areas likely to be impacted by ash. A new GIS-based tool was recently added to the HOTVOLC system (Guéhenneux & Gouhier, 2024), allowing the user to directly measure the plume on the METEOSAT images, at each time step (i.e., every 15 min). To test whether our method could be used in near real-time, we estimated a posteriori the source MER for ten well-observed explosive eruptions of Mt Etna (Italy) between 2015 and 2022. Figure 2b shows a screenshot of the HOTVOLC web interface displaying the METEOSAT satellite images (MSG-SEVIRI) taken over Mt Etna during the 21 February 2022 eruption at 12:45 UTC. The plume produced by the explosive eruption was clearly visible on the images starting at 11:30 UTC, and we show in Figure 2b the six measurements of both the plume width *L* and the downwind distance *x* made on each image between 11:30 and 12:45 UTC.

### 3.4. Uncertainties

Uncertainties in our MER estimations can particularly arise from (a) the measurements of the plume geometry made on satellite data, and thus from (b) the comparison between these measurements and the model predictions. By repeating the measurements on GOES-16 images for La Soufrière several times, we estimated a variation of about 10% of the values of *L* and *x* between each measurement (see Figure S1 and Table S1 in Supporting Information S1). We thus calculated an uncertainty of the MER of about 5% when we overestimate the values of *L* and *x*, and of 20% when we underestimate them (see Figure S2 in Supporting Information S1). For Mt Etna, the new GIS-based tool added to the HOTVOLC system allows the measurements to be more precise, with a variation of only 5% between two measurements made on the same image (see Table S2 in Supporting Information S1). We thus estimated our MER to have an uncertainty of about 6.5% when we overestimate the values of *L* and *x*, and of 2.5% when we underestimate them (see Figure S3 in Supporting Information S1). Details about the method used to characterize these uncertainties are given in Text S1 and S2 in the Supporting Information S1.

### 4. Results

We tested our updated version of PPM by varying several input parameters in order to investigate their effects on the theoretical predictions of plume growth (i.e., the evolution of the plume width *L* as a function of the downwind distance *x*). We chose to vary four main input parameters likely to influence the evolution of the plume: the type of wind profile together with the maximum wind speed, the empirical parameter *μ* introduced by Ishii et al. (2022) and linking the predictions at the NBH with the plume geometry far from the source, the granulometry (i.e., the total grain-size distribution, TGSD), and the MER at the source of the eruption. All these parameters (exception made of the parameter *μ*) influence the total maximum height reached by the volcanic column (Girault et al., 2014, 2016). We also test the model's ability to estimate the MER at the source of an explosive eruption by adapting
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

**Figure 3.** Theoretical predictions for the evolution of the plume width, *L* (m) as a function of the downwind distance, *x* (m) when using (a) a uniform wind profile (shown in inset) with a wind speed of 10 m s⁻¹ (in blue) or 20 m s⁻¹ (in turquoise); (b) A linear wind profile (shown in inset) with a maximum wind speed of 10 m s⁻¹ (in blue), 20 m s⁻¹ (in turquoise) or 30 m s⁻¹ (in pink); (c) A complex wind profile (shown in inset) with a maximum wind speed of 10 m s⁻¹ (in blue), 20 m s⁻¹ (in turquoise) or 30 m s⁻¹ (in pink). The triangles on each curve show the downwind distance *x* at which the NBH was reached in the calculations. Note that the blue curve in (b) is the same as the blue curves in Figures 4 and 5, and the orange one in Figure 6.

PPM to the environmental conditions of La Soufrière, St Vincent and Mt Etna, and by using the measurements made on satellite images (Section 3.3).

## 4.1. Effect of Wind Profile

Figure 3 shows the evolution of the plume width as a function of the downwind distance (i.e., with time) when considering the uniform (Figure 3a), linear (Figure 3b), and complex (Figure 3c) wind profile (the wind profiles are given in the insets). Regardless of the wind type profile, we note that the plume growth predictions show two behaviors, with first a steady increase of the width with time (during the vertical ascent of the column toward the NBH), and then acceleration in the width evolution in the turbulent current region (at ≈ 1 km in horizontal distance from the eruptive vent). This effect is reinforced as the maximum wind speed is increased in the model, and can be explained by the greater contribution of *w* in Equation 19, compared with the evolution of the column width that relies mostly on the entrainment of ambient air. For a given downwind distance *x*, the plume width is increased by a factor of 1.8 when decreasing the maximum wind speed from 20 to 10 m s⁻¹, and by a factor of 2.6 when choosing 10 instead of 30 m s⁻¹. The predictions made when choosing a linear wind profile (Figure 3b) illustrate well this effect of the maximum wind speed at the tropopause (and from the eruptive vent to the turbulent current region, in contrast to the complex wind profile, Figure 3c), and are robust even for a strong wind speed (unlike the uniform wind profile, Figure 3a). We also note a crossover between the three curves in Figure 3c where the predictions made for a wind speed of 30 m s⁻¹ at the NBH cross the 10 m s⁻¹ curve.

## 4.2. Effect of Entrainment at Plume Top

Figure 4 shows the evolution of the plume width with the downwind distance from the eruptive source when varying the empirical parameter $\mu$ (see Equation 20) from 1 to 4. The parameter $\mu$ intervenes only in the calculations for the turbulent current region; its impact thus being solely seen after the NBH has been reached (i.e., at a downwind distance of about 1 km, shown by the triangle in Figure 4). The effect of this entrainment at the plume top is opposite to the wind speed effect, that is, the turbulent current width increases with $\mu$ (while it decreased when the maximum wind speed increased). For a given downwind distance $x$, the turbulent current width is increased by a factor of 1.5 when considering $\mu = 2$ instead of 1, and even doubled when considering $\mu = 4$. The effect of $\mu$ is thus critical on the model predictions. We will further discuss the implications of this impact in the discussion section (Section 5.1).

## 4.3. Effect of Granulometry

Given the strong influence of the TGSD on both the maximum column height (Girault et al., 2014) and column collapse (Michaud-Dubuy et al., 2018), we investigated its impact on the evolution of the plume geometry. The results are shown in Figure 5. Choosing different TGSDs at the eruptive vent (i.e., fine, mixed or coarse) induces very little change in the predictions. This is not surprising for the turbulent current region predictions (after about 1 km in horizontal distance) as most of the column is mostly composed of fine particle at the NBH (the coarser ones having sedimented before). It shows, however, that the TGSD has a limited impact on the entrainment rate in the column region, as the column width is not influenced by the choice of the power-law exponent *D*.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Theoretical predictions for the evolution of the plume width, *L* (m) as a function of the downwind distance, *x* (m) when using an empirical parameter *μ* of 1 (in blue), 2 (in turquoise), 3 (in pink), or 4 (in light pink). The triangle shows the downwind distance *x* at which the NBH was reached in the calculations.</figcaption></figure>

## 4.4. Effect of Mass Eruption Rate (MER)

We finally test the effect of the MER by calculating the evolution of the plume width as a function of the downwind distance for seven values of the MER, ranging from $10^3$ to $10^9$ kg s$^{-1}$ (Figure 6). The mass eruption rate is found to be the main parameter influencing the evolution of the plume geometry, with two orders of magnitude difference between the plume width calculated for $10^3$ kg s$^{-1}$ and the one calculated for $10^9$ kg s$^{-1}$ for a given value of *x*. The influence of the MER on the predictions also seems to increase with the MER. Indeed, the larger the MER, the larger the gap between predictions for two MER values. For example, for a given value of *x*, there is a factor 1.3 between the plume widths calculated for $10^3$ and $10^4$ kg s$^{-1}$, but there is a factor 6 between the plume width calculated for $10^8$ and $10^9$ kg s$^{-1}$. The predictions made for strong MER ($\ge 10^7$) should however be taken with caution as they are the only simulations for which we calculated $w^*/u^*$ ratios (Equation 11) corresponding to strong plumes. Equations 17–19 are however only valid for plumes spreading downwind, and not for radial spreading turbulent currents (for which the model developed by F. Prata et al. (2024), e.g., would be more appropriate).

In the following section, we show how the strong variations in the predictions for different values of the MER presented in Figure 6 allow us to effectively compare the measurements made on satellite images for observed downwind plumes with our calculations, and thus estimate the corresponding source MER for each explosive event.

## 4.5. Implications for La Soufrière, St Vincent Eruption

Figure 7a shows the theoretical predictions for three MERs ($10^5$, $10^6$, and $10^7$ kg s$^{-1}$), when the initial conditions of PPM are adapted to the 2021 eruption of La Soufrière, St Vincent. Based on the location of the Saint-Vincent island, we considered a tropical atmosphere, as well as a linear wind profile with a maximum wind speed of 20 m s$^{-1}$ at the tropopause. This is consistent with both the mean wind profile calculated for the duration of the eruption by Taylor et al. (2023) (their Figure 3b), and with the mean wind speeds we calculated from our measurements on satellite images (see Section 3.3, and Table S3 in Supporting Information S1). The three simulations correspond to the formation of weak to distorted plumes ($w^*/u^*$ ratios > 0.11, see Equation 11). We compared these predictions with the plume measurements that we made on the satellite images (GOES-16) for the 33 explosions that occurred

between 9 and 14 April 2021 (see Section 3.3). Three examples of these measurements are given in Figure 7a, for events 6 (10 April–07:50 UTC), 19 (11 April–08:10 UTC), and 28 (14 April–02:40 UTC). We see that the measurements collapse on the theoretical curves, suggesting that the model captures well the observations. We then estimated a MER at the source for each of the 33 explosions (Table 2). The min. MER and max. MER calculated using the uncertainties estimates in Section 3.4 are also given in Table 2.

Comparing our MER estimates with those of Sparks et al. (2023) (Figure 7b) provides a good agreement for the first two stages of the eruption, except for event 5 (10 April–06:50) for which we estimate a MER significantly below the one estimated by Sparks et al. (2023). As it is clearly visible on the satellite images that this event produced a rather weaker plume than the other explosions with a considerably smaller width, we are confident in our MER estimate (see Figure S4 in Supporting Information S1 for details). We suggest that this difference in the MER estimates could come from the fact that the rather powerful explosion that we see on satellite images shortly before (at 05:20 UTC, Figure S4 in Supporting Information S1) was not detected on RSAM data from Sparks et al. (2023), and that perhaps this additional event is included within the analysis of event 5 by Sparks et al. (2023) (although we note that the onset time for this event is posterior to our additional event in Sparks

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Theoretical predictions for the evolution of the plume width, *L* (m) as a function of the downwind distance, *x* (m) when using a fine TGSD (i.e., a power-law exponent *D* = 3.3; in blue), a mixed TGSD (i.e., *D* = 3.0; in turquoise), and a coarse TGSD (i.e., *D* = 2.8; in pink). The triangles on each curve show the downwind distance *x* at which the NBH was reached in the calculations.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Theoretical predictions for the evolution of the plume width, *L* (m) as a function of the downwind distance, *x* (m) when using a MER of 10³ (in blue), 10⁴ (in turquoise), 10⁵ (in pink), 10⁶ (in orange), 10⁷ (in purple), 10⁸ (in gray), and 10⁹ kg s⁻¹ (in light blue). The triangles on each curve show the downwind distance *x* at which the NBH was reached in the calculations. The dashed curves represent the only simulations corresponding to strong plumes and for which our calculation of the plume geometry could thus not be accurate.</figcaption></figure>

et al. (2023)). Our MER estimates also clearly separates from other estimates for the stage 3 (Figure 7b), although they remain in the same order of magnitude (with the exception of events 22 and 29). We thus test another model to understand this divergence for stage 3. Comparing the same plume measurements with the theoretical predictions from Pouget et al. (2013) (their Figure 7), we estimated the MER for each of the 33 explosive events of La Soufrière (see Table S3 in Supporting Information S1). Figure 8 shows the difference in MER predictions ΔMER ((MER - MER<sub>Sparks</sub>)/MER<sub>Sparks</sub>, where MER represents our estimates using either PPM or the predictions of Pouget et al. (2013), as a function of time. The predictions of Pouget et al. (2013) lead to underestimated values of MER for the stages 1 and 2 of the eruption, while PPM compares better overall with the data from Sparks et al. (2023) with a mean ΔMER of -0.23. This is especially true for stage 2, while we note for stage 1 that PPM shows greater overall deviations from the MER predictions of Sparks et al. (2023) than Pouget et al. (2013). For the stage 3, MER predictions using the model of Pouget et al. (2013) are closer to the estimations of Sparks et al. (2023). We discuss this discrepancy in Section 5.2.

Our results show that the model provides MER estimates in good agreement with field-based data. Our predictions are also consistent at first order with those made using previous similar models in the literature (Pouget et al., 2013), while providing a good refinement by taking into account the column dynamics from the eruptive source.

## 4.6. Real-Time MER Estimates Using the HOTVOLC Platform

We compile in Table 3 the 10 explosive events of Mt Etna for which we compared plume growth measurements with the theoretical predictions from PPM. Four examples are given in Figure 9, for the eruptions of 4 December 2015, 24 February 2021, 7 March 2021, and 21 February 2022. As for La Soufrière eruption, the data collapse well on the theoretical curves. To decide whether to use the theoretical predictions made for a linear wind profile with a maximum wind speed of 10 (in color) or 20 m s⁻¹ (in gray) to estimate the MER, we used both the mean wind speed during the explosion determined from our downwind distance measurements with time (see Table 3) and the mean wind speed calculated by Mereu et al. (2023) using the ERA5 reanalysis over a larger scale. The star in the first column of Table 3 indicates the events for which we used the 20 m s⁻¹ predictions. The estimated MER are given in Table 3, and the min. MER and max. MER calculated when considering the overall uncertainties are given in Table S3 of Supporting Information S1. Comparing the estimates of 8 of the 10 eruptions with MER calculated by Mereu et al. (2023) or Scollo et al. (2025) using radar data (Table 3) provides a good agreement between our estimates that remain within the same order of magnitude with a mean ΔMER of 0.6. This corresponds to a mean factor of 1.6 between our MER estimates and those from radar measurements (and a maximum factor of 2.5, Figure 10a). The comparison between our wind speed estimates and those from the ERA5 reanalysis by Mereu et al. (2023) are also in the same order of magnitude (Figure 10b). These results confirm that our methodology can be used for the rapid assessment of MER during future eruptive events.

# 5. Discussion

## 5.1. Relationship Between Volume Flux Ratio and the Plume Regime

In this study, we chose to use an empirical parameter μ linking the volumetric flux of the column at the NBH, $V_0$ with the volumetric flux feeding the downwind turbulent current, $V_1$ (Ishii et al., 2022) set to 1. Figure 4, however, shows that varying μ from 1 to 4 changes the predictions of the turbulent current geometry evolution, and thus the estimation of the MER at the source. Indeed, increasing μ up to 4 would reduce the estimated MER by one order of magnitude. In this study, we consider only weak and distorted plumes, with either no or limited overshoot height compared to the strong plume case. In the weak plume case, the transition from the bended column to the downwind turbulent current is smooth and we expect $V_0 = V_1$ (i.e., μ = 1). In the distorted plume case, the overshoot region being limited, there is no significant downflow that could increase the entrainment and thus the
<figure><img src="image_10.png" /><figcaption>**Figure 7.** (a) Evolution of the plume width, *L* (m) as a function of the downwind distance, *x* (m) measured on satellite images (GOES-16) for three explosive events of La Soufrière, St Vincent (blue diamond, event 6; orange circle, event 19; black square, event 28; see Table 1 for event dates) compared with PPM's theoretical predictions when considering a MER of 10⁵ (black solid line), 10⁶ (dashed black line), or 10⁷ kg s⁻¹ (dotted black line), a fine TGSD (i.e., power-law exponent *D* = 3.3), a linear wind profile with a maximum wind speed of 20 m s⁻¹, a tropical atmosphere, a gas content of 4 wt%, and an empirical parameter *μ* set at 1. (b) Comparison between the MER estimated by Sparks et al. (2023) (purple circles) and the ones estimated in this study using PPM's theoretical predictions (red squares), from 9 to 14 April 2021.</figcaption></figure>

volumetric flow feeding the turbulent current. Indeed, Suzuki and Koyaguchi (2009) suggested that the entrainment of air above the NBH mainly occurs at the very top of the plume (in the overshoot region due to the downflow from the maximum height toward the downwind turbulent current) rather than from the top/at the bottom of the downwind turbulent current. This assumption suggests that $V_0 \approx V_1$ (i.e., $\mu \approx 1$) for a distorted plume, while one can expect for strong plumes with large overshoot heights that $V_1 > V_0$ (i.e., $\mu > 1$) (Suzuki et al., 2005). This highlights the need to improve our knowledge about entrainment processes at the plume top region in order to apply our model to the specific case of strong plume (Ishii et al., 2022).

Our results suggest that $\mu = 1$ is a valid assumption for weak ($w^*/u^* > 0.37$, see Equation 11) and distorted plumes ($0.11 \le w^*/u^* \le 0.37$). However, the case studies investigated here are relatively similar (MER, and wind profile and speed), suggesting that the parameterization of $\mu$ could be improved by investigating further cases. Changing the MER and maximum wind speed values may indeed highlight a possible variation of $\mu$ with these parameters (i.e., with the plume regime). Laboratory studies simulating turbulent plumes in a windy environment (Carazzo et al., 2014; Michaud-Dubuy et al., 2020) may help to study the effect of the $w^*/u^*$ ratio on the effective value of $\mu$, and refine the model for the whole spectrum of volcanic plume behavior.
<caption>**Table 2**</caption>

Mass Eruption Rates Estimated Using Our New Version of PPM, and Their Uncertainties for Explosive Events Between 9 and 14 April 2021 at La Soufrière Volcano, St Vincent

<table><thead><tr><th>Event ID</th><th>Date and time (DD/MM/AA UTC)</th><th>MER (kg s<sup>-1</sup>)</th><th>Min. MER</th><th>Max. MER</th></tr></thead><tbody><tr><td>1</td><td>09/04/21 12:50</td><td>3.15 E + 06</td><td>2.99 E + 06</td><td>3.78 E + 06</td></tr><tr><td>2</td><td>09/04/21 19:10</td><td>5.02 E + 06</td><td>4.77 E + 06</td><td>6.02 E + 06</td></tr><tr><td>PO1</td><td>09/04/21 20:20</td><td>5.84 E + 06</td><td>5.55 E + 06</td><td>7.01 E + 06</td></tr><tr><td>PO2</td><td>09/04/21 21:30</td><td>6.14 E + 06</td><td>5.83 E + 06</td><td>7.37 E + 06</td></tr><tr><td>3</td><td>09/04/21 22:40</td><td>1.35 E + 07</td><td>1.28 E + 06</td><td>1.62 E + 07</td></tr><tr><td>PO7</td><td>09/04/21 23:40</td><td>1.06 E + 07</td><td>1.01 E + 07</td><td>1.27 E + 07</td></tr><tr><td>4-PO16</td><td>10/04/21 00:20</td><td>1.49 E + 07</td><td>1.42 E + 07</td><td>1.79 E + 07</td></tr><tr><td>/</td><td>10/04/21 05:20</td><td>7.09 E + 06</td><td>6.74 E + 06</td><td>8.51 E + 06</td></tr><tr><td>5</td><td>10/04/21 06:50</td><td>1.80 E + 05</td><td>1.71 E + 05</td><td>2.16 E + 05</td></tr><tr><td>6</td><td>10/04/21 07:50</td><td>5.56 E + 06</td><td>5.28 E + 06</td><td>6.67 E + 06</td></tr><tr><td>7</td><td>10/04/21 09:50</td><td>5.02 E + 06</td><td>4.77 E + 06</td><td>6.02 E + 06</td></tr><tr><td>8</td><td>10/04/21 11:00</td><td>6.14 E + 06</td><td>5.83 E + 06</td><td>7.37 E + 06</td></tr><tr><td>9</td><td>10/04/21 12:10</td><td>8.78 E + 06</td><td>8.34 E + 06</td><td>1.05 E + 07</td></tr><tr><td>10</td><td>10/04/21 13:00</td><td>1.06 E + 07</td><td>1.01 E + 07</td><td>1.27 E + 07</td></tr><tr><td>11</td><td>10/04/21 14:40</td><td>8.15 E + 06</td><td>7.74 E + 06</td><td>9.78 E + 06</td></tr><tr><td>12</td><td>10/04/21 16:30</td><td>1.01 E + 07</td><td>9.60 E + 06</td><td>1.21 E + 07</td></tr><tr><td>13</td><td>10/04/21 18:50</td><td>3.53 E + 06</td><td>3.35 E + 06</td><td>4.24 E + 06</td></tr><tr><td>14</td><td>10/04/21 21:30</td><td>9.73 E + 06</td><td>9.24 E + 06</td><td>1.17 E + 07</td></tr><tr><td>15</td><td>10/04/21 23:10</td><td>1.70 E + 07</td><td>1.62 E + 07</td><td>2.04 E + 07</td></tr><tr><td>16</td><td>11/04/21 01:00</td><td>5.56 E + 06</td><td>5.28 E + 06</td><td>6.67 E + 06</td></tr><tr><td>17</td><td>11/04/21 02:50</td><td>1.06 E + 07</td><td>1.01 E + 07</td><td>1.27 E + 07</td></tr><tr><td>18</td><td>11/04/21 05:10</td><td>7.09 E + 06</td><td>6.74 E + 06</td><td>8.51 E + 06</td></tr><tr><td>19</td><td>11/04/21 08:10</td><td>6.76 E + 06</td><td>6.42 E + 06</td><td>8.11 E + 06</td></tr><tr><td>20</td><td>11/04/21 10:50</td><td>2.23 E + 06</td><td>2.12 E + 06</td><td>2.68 E + 06</td></tr><tr><td>21</td><td>11/04/21 13:40</td><td>6.76 E + 06</td><td>6.42 E + 06</td><td>8.11 E + 06</td></tr><tr><td>22</td><td>11/04/21 18:20</td><td>7.85 E + 05</td><td>7.46 E + 05</td><td>9.42 E + 05</td></tr><tr><td>23</td><td>11/04/21 20:20</td><td>1.09 E + 06</td><td>1.04 E + 06</td><td>1.31 E + 06</td></tr><tr><td>24</td><td>12/04/21 00:50</td><td>3.23 E + 06</td><td>3.07 E + 06</td><td>3.88 E + 06</td></tr><tr><td>25</td><td>12/04/21 08:10</td><td>3.62 E + 06</td><td>3.44 E + 06</td><td>4.34 E + 06</td></tr><tr><td>26</td><td>12/04/21 21:30</td><td>1.48 E + 06</td><td>1.41 E + 06</td><td>1.78 E + 06</td></tr><tr><td>27</td><td>13/04/21 10:40</td><td>2.87 E + 06</td><td>2.73 E + 06</td><td>3.44 E + 06</td></tr><tr><td>28</td><td>14/04/21 02:40</td><td>1.27 E + 06</td><td>1.21 E + 06</td><td>1.52 E + 06</td></tr><tr><td>29</td><td>14/04/21 15:40</td><td>7.85 E + 05</td><td>7.46 E + 05</td><td>9.42 E + 05</td></tr></tbody></table>

Note. The table is separated in three parts representing the three first stages of the eruption (Sparks et al., 2023).

## 5.2. Comparison With MER Estimates From Other Methods

The comparison between different methods to retrieve the MER exhibits some discrepancies for both La Soufrière, St Vincent and Mt Etna eruptions. Thus, for the stage 3 of La Soufrière eruption, even if the predictions from Pouget et al. (2013) are closer to the estimates from Sparks et al. (2023) than PPM, there is still up to a factor 2 of difference. This discrepancy between the three methods suggests that we lack some understanding about the overall plume dynamics during this stage 3. The analysis of field deposits by Cole et al. (2023) show that the products of the stage 3 explosions, attributed to the deposit units U6 and U7, exhibit finer particle sizes compared to the other units. As theoretical predictions of 1D models of volcanic columns have shown that finer material can
<figure><img src="image_11.png" /><figcaption>**Figure 8.** ΔMER as a function of time, from first explosion on 9 April 2021 at 12:50 UTC, where ΔMER = (MER - MER<sub>Sparks</sub>)/MER<sub>Sparks</sub>, with MER representing our estimates using whether PPM (in red) or the predictions of Pouget et al. (2013) (in blue). A positive ΔMER means that the MER estimated in this study is greater than the MER estimated by Sparks et al. (2023) while a negative ΔMER means that the estimate of Sparks et al. (2023) is greater. The dotted lines show the mean values of ΔMER for each data set.</figcaption></figure>

lead to higher maximum column heights for the same MER at the source (Girault et al., 2014), we argue that the MER calculated from the maximum column heights by Sparks et al. (2023) could be overestimated for this final stage of the eruption.

Overall, our MER estimates are consistent with those of Sparks et al. (2023) with a good to very good agreement for the first two stages of the eruption, while remaining in the same order of magnitude for stage 3. The mean ΔMER calculated from the two sets of MER estimates compared with the ones from Sparks et al. (2023) (dashed lines in Figure 8) show that the predictions of PPM lead to an overall better consistency with previously published data. The evolution of the estimated MER with time is also consistent with the characteristics of each stage of the eruption, that is, a first stage with an increasing activity, a second stage with powerful and recurrent explosions, and a third stage with decreasing intensity and longer repose periods. These results allow us to validate this new version of PPM for the estimation of source MER from plume growth measurements on satellite images.

<caption>**Table 3**</caption>

Estimated Mass Eruption Rates and Mean Wind Speed, Compared to the Ones by Mereu et al. (2023) and Scollo et al. (2025), for Ten Explosive Events at Mt Etna Between 2015 and 2022

<table><thead><tr><th rowspan="2">Event<br>(DD/MM/AAAA UTC)</th><th colspan="2">This study</th><th colspan="2">Mereu et al. (2023)</th></tr><tr><th>MER (kg s<sup>-1</sup>)</th><th>Mean wind speed (m s<sup>-1</sup>)</th><th>MER (kg s<sup>-1</sup>)</th><th>Mean wind speed (m s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>04/12/2015 09:15</td><td>3.46 E + 06</td><td>11.4</td><td>1.46 E + 06</td><td>5</td></tr><tr><td>21/12/2020 09:15</td><td>1.07 E + 06</td><td>8.7</td><td>/</td><td>/</td></tr><tr><td>24/02/2021 19:45</td><td>1.65 E + 05</td><td>6.3</td><td>2.19 E + 05</td><td>5.8</td></tr><tr><td>28/02/2021 08:00</td><td>2.00 E + 06</td><td>15.8</td><td>2.81 E + 06</td><td>12.5</td></tr><tr><td>07/03/2021 06:30*</td><td>3.17 E + 06</td><td>23.3</td><td>1.55 E + 06</td><td>19.2</td></tr><tr><td>12/03/2021 08:00</td><td>4.00 E + 05</td><td>13</td><td>4.50 E + 05<sup>a</sup></td><td>/</td></tr><tr><td>09/08/2021 02:00</td><td>4.58 E + 05</td><td>14.2</td><td>/</td><td>/</td></tr><tr><td>23/10/2021 09:00*</td><td>4.44 E + 06</td><td>14.3</td><td>1.77 E + 06</td><td>15.8</td></tr><tr><td>10/02/2022 20:30</td><td>1.42 E + 06</td><td>9.1</td><td>8.20 E + 05</td><td>7.9</td></tr><tr><td>21/02/2022 11:30</td><td>1.63 E + 06</td><td>9</td><td>8.74 E + 05</td><td>8.8</td></tr></tbody></table>

Note. *Events for which we used the theoretical predictions for a linear wind profile with a max. wind speed of 20 m s<sup>-1</sup> to estimate the MER. <sup>a</sup>Radar estimate from Scollo et al. (2025).
<figure><img src="image_12.png" /></figure>

**Figure 9.** Evolution of the plume width, *L* (m) as a function of the downwind distance, *x* (m) measured on satellite images (MSG-SEVIRI) for four explosive eruptions of Mt Etna (orange square, 4 December 2015 09:15 UTC; lilac triangle, 24 February 2021 19:45 UTC; gray circle, 7 March 2021 06:30 UTC; purple diamond, 21 February 2022 11:30 UTC) compared with PPM's theoretical predictions when considering a linear wind profile and a MER of 10⁴ (blue solid line for a maximum wind speed of 10 m s⁻¹, gray solid line for a maximum wind speed of 20 m s⁻¹), 10⁵ (turquoise solid line for a max wind speed of 10 m s⁻¹, dashed gray line for a max wind speed of 20 m s⁻¹), or 10⁶ kg s⁻¹ (pink solid line for a max wind speed of 10 m s⁻¹, and dotted gray line for a max wind speed of 20 m s⁻¹). Calculations were made considering a fine TGSD (i.e., power-law exponent *D* = 3.3), a mid-latitude atmosphere, a gas content of 4 wt%, and an empirical parameter *μ* set at 1.

For the case of Mt Etna, the slight overestimation of the MER compared to radar data could be attributed to the lava fountaining activity usually occurring at this volcano and inducing a partitioning of the MER between the one feeding the lava fountain and the one feeding the volcanic column (Calvari et al., 2022; Snee et al., 2021). This example highlights the transient nature of explosive eruptions, which can vary rapidly in style during the same event, making it challenging to retrieve critical ESPs for real-time hazard assessment. While the improvement of the method (including a better estimate of the empirical parameter μ, and an automatized integration within the HOTVOLC system) will be the object of a future work, the results suggest that we can already use PPM and the current Web-GIS version of HOTVOLC to provide a first estimate of the source MER during an eruption.

## 6. Conclusions

We presented in this study a refined version of the 1-D volcanic column model PPM now including a simplified downwind gravity current model. We perform a sensitivity study of the model by varying the wind profile type and maximum wind speed, the entrainment at the top of the volcanic column (i.e., the empirical parameter μ introduced by Ishii et al. (2022)), the total grain-size distribution and the mass eruption rate. The main parameters affecting the theoretical predictions of the plume geometry far from the eruptive source are the maximum wind speed, the parameter μ and the MER. We then used GOES-16 images of the well observed and studied 2021 eruption of La Soufrière, St Vincent to measure the evolution of the plume width and downwind distance for 33 explosions between 9 and 22 April 2022, and compare them with our theoretical predictions to estimate the source MER of

each explosion. We found a good agreement between our estimates and those of Sparks et al. (2023) for the first two stages of the eruptions (despite greater deviations from the MER estimates of Sparks et al. (2023) for stage 1), while they remained in the same order of magnitude for the third stage. To explain this discrepancy, we used the same satellite measurements to estimate the source MER using the method of Pouget et al. (2013), and found the same difference for the stage 3 suggesting we lack some understanding about this waning stage of the eruption. Overall, our model provided a good refinement compared to the previous method of Pouget et al. (2013). We finally tested whether we could use our model to estimate the source MER in near real-time by using the Web-GIS HOTVOLC system, and determined the source MER of 10 paroxysms of Mt Etna between 2015 and 2022. Using a new GIS-based tool of measurement on the HOTVOLC web interface allows a strong reduction of the uncertainty in the measurements, and thus on the MER estimates. Our results compares well with other MER estimates found in the literature, giving us confidence in a future integration of our method into the HOTVOLC system. This method only requiring operational measurements of the plume width and length largely reduces the

<figure><img src="image_13.png" /></figure>

**Figure 10.** (a) Comparison between the modeled MER estimated in this study and those measured using radar data, in kg s⁻¹ (Mereu et al., 2023; Scollo et al., 2025). Error bars for the modeled MER are comprised within the symbols; error bars for the measured MER indicate an uncertainty of 20% based on the error propagation described by Mereu et al. (2023).
(b) Comparison between the wind speeds estimated from our downwind distance measurements with time and the wind speed estimated by (Mereu et al., 2023) from the ERA5 reanalysis, in m s⁻¹. Dashed lines indicate over- or under-estimations by a factor of 2.
uncertainties in the MER estimates, compared to the cloud-top temperature method currently used. This will add another high-value product to the ones already provided by HOTVOLC to the Toulouse VAAC, and allow more precise mapping of tephra dispersal while allowing the creation of a database including several eruptive source parameters determined by the OPGC.

Notation

<table><tr><td>A</td><td>buoyancy parameter, dimensionless</td></tr><tr><td>c<sub>a</sub></td><td>specific heat at constant pressure of the atmosphere, 998 JK<sup>-1</sup></td></tr><tr><td>c<sub>p</sub></td><td>specific heat at constant pressure of the particles, 1100 JK<sup>-1</sup></td></tr><tr><td>D</td><td>exponent of the power-law distribution of particles, dimensionless</td></tr><tr><td>d<sub>φ</sub></td><td>particle diameter, m</td></tr><tr><td>F<sub>0</sub></td><td>source buoyancy flux, m<sup>4</sup> s<sup>-3</sup></td></tr><tr><td>g</td><td>gravitational acceleration, 9.81 m s<sup>-2</sup></td></tr><tr><td>H</td><td>downwind turbulent current thickness, m</td></tr><tr><td>L</td><td>plume width, m</td></tr><tr><td>N</td><td>Brunt-Väisälä frequency, 0.01 s<sup>-1</sup></td></tr><tr><td>N<sub>φ</sub></td><td>number of classes of particle size, 20</td></tr><tr><td>P<sub>a</sub></td><td>atmospheric pressure, Pa</td></tr><tr><td>p<sub>s</sub></td><td>probability of sedimentation, 0.27 ± 0.01, dimensionless</td></tr><tr><td>Q<sub>φ</sub></td><td>mass discharge rate of particles, kg s<sup>-1</sup></td></tr><tr><td>r</td><td>column radius, m</td></tr><tr><td>R<sub>a</sub></td><td>gas constant of the atmospheric air, 287 J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>R<sub>g</sub></td><td>bulk column of the mixture in the column, R<sub>g</sub> = 461 J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>R<sub>p</sub>, r<sub>p</sub></td><td>particle sizes used to describe the power-law fragmentation</td></tr><tr><td>Re<sub>p</sub></td><td>particle Reynolds number, dimensionless</td></tr><tr><td>Ri</td><td>local Richardson number, dimensionless</td></tr><tr><td>s</td><td>curvilinear centerline of the plume, m</td></tr><tr><td>T</td><td>temperature of the volcanic mixture, T<sub>0</sub> = 1200 K</td></tr><tr><td>T<sub>a</sub></td><td>ambient atmospheric temperature for a mid-latitude atmosphere, 273 K</td></tr><tr><td>u</td><td>column vertical velocity, m s<sup>-1</sup></td></tr><tr><td>u<sub>e</sub></td><td>entrainment rate at the edge of the plume, m s<sup>-1</sup></td></tr><tr><td>U</td><td>characteristic velocity of the downwind turbulent current, m s<sup>-1</sup></td></tr><tr><td>V<sub>φ</sub></td><td>particle fallout velocity, m s<sup>-1</sup></td></tr><tr><td>V<sub>1</sub></td><td>volumetric flux of the downwind turbulent current, m<sup>3</sup> s<sup>-1</sup></td></tr><tr><td>w</td><td>wind speed, m s<sup>-1</sup></td></tr><tr><td>x</td><td>downwind distance from the eruptive vent, m</td></tr><tr><td>x<sub>g</sub></td><td>effective gas mass fraction, dimensionless</td></tr><tr><td>x<sub>φ</sub></td><td>mass proportion of φ-sized particles in the grain-size distribution, dimensionless</td></tr></table>
<table><tr><td>x<sub>tot</sub></td><td>total gas content, dimensionless</td></tr><tr><td>z</td><td>vertical axis/altitude, m</td></tr><tr><td>α</td><td>variable entrainment coefficient, dimensionless</td></tr><tr><td>α<sub>p</sub></td><td>“top-hat” entrainment coefficient for a pure plume, 0.125, dimensionless</td></tr><tr><td>β</td><td>wind entrainment coefficient, 0.1</td></tr><tr><td>θ</td><td>inclination of the plume centerline relative to the horizontal, rad</td></tr><tr><td>λ</td><td>normalization constant of the power-law distribution, dimensionless</td></tr><tr><td>λ<sub>it</sub></td><td>spreading constant depending on the flow geometry, dimensionless</td></tr><tr><td>μ<sub>air</sub></td><td>dynamic viscosity of air, μ<sub>0</sub> = 1.832 × 10<sup>-5</sup> Pa s</td></tr><tr><td>μ</td><td>ratio of volumetric flux, μ = V<sub>1</sub>/V<sub>0</sub>, dimensionless</td></tr><tr><td>φ</td><td>particle size notation, d<sub>φ</sub> = 2<sup>-φ</sup></td></tr><tr><td>ρ</td><td>bulk density of the volcanic mixture, kg m<sup>-3</sup></td></tr><tr><td>ρ<sub>a</sub></td><td>atmospheric density, kg m<sup>-3</sup></td></tr><tr><td>ρ<sub>p</sub></td><td>particle density, 2000 kg m<sup>-3</sup></td></tr><tr><td>o</td><td>values at the eruptive vent</td></tr></table>

## Data Availability Statement

The Supporting Information, available in Figshare via Michaud-Dubuy and Gouhier (2025), contains all mass eruption rates estimated in this study. The measurements of the plume geometry for each explosion of La Soufrière, St Vincent were made on GOES-16 satellite imagery combining the infrared anomaly with ash detected in visible light (“ABI RadF ash”) available as a supplementary file in Cole et al. (2023). The measurements for Mt Etna eruptions were made on the freely accessible HOTVOLC platform (Gouhier et al., 2016).

## References

* Aubry, T., & Jellinek, A. (2018). New insights on entrainment and condensation in volcanic plumes: Constraints from independent observations of explosive eruptions and implications for assessing their impacts. *Earth and Planetary Science Letters, 490*, 132–142. https://doi.org/10.1016/j.epsl.2018.03.028

* Aubry, T., Jellinek, A., Carazzo, G., Gallo, R., Hatcher, K., & Dunning, J. (2017). A new analytical scaling for turbulent wind-bent plumes: Comparison of scaling laws with analog experiments and a new database of eruptive conditions for predicting the height of volcanic plumes. *Journal of Volcanology and Geothermal Research*, 343, 233–251. https://doi.org/10.1016/j.jvolgeores.2017.07.006

* Bonadonna, C., Ernst, G., & Sparks, R. (1998). Thickness variations and volume estimates of tephra fall deposits: The importance of particle Reynolds number. *Journal of Volcanology and Geothermal Research*, 81(3–4), 173–187. https://doi.org/10.1016/S0377-0273(98)00007-9

* Bonadonna, C., & Phillips, J. (2003). Sedimentation from strong volcanic plumes. *Journal of Geophysical Research: Solid Earth*, 108(B7), 1–28. https://doi.org/10.1029/2002JB002034

* Bursik, M. (2001). Effect of wind on the rise height of volcanic plumes. *Geophysical Research Letters*, 28(18), 3621–3624. https://doi.org/10.1029/2001GL013393

* Bursik, M., Sparks, R., Gilbert, J., & Carey, S. (1992). Sedimentation of tephra by volcanic plumes: I. Theory and its comparison with a study of the Fogo A plinian deposit, Sao Miguel (Azores). *Bulletin of Volcanology*, 54(4), 329–344. https://doi.org/10.1007/bf00301486

* Calvari, S., Biale, E., Bonaccorso, A., Cannata, A., Carleo, L., Currenti, G., et al. (2022). Explosive paroxysmal events at etna volcano of different magnitude and intensity explored through a multidisciplinary monitoring system. *Remote Sensing*, 14(16), 4006. https://doi.org/10.3390/rs14164006

* Carazzo, G., Girault, F., Aubry, T., Bouquerel, H., & Kaminski, E. (2014). Laboratory experiments of forced plumes in a density-stratified crossflow and implications for volcanic plumes. *Geophysical Research Letters*, 41(24), 8759–8766. https://doi.org/10.1002/2014GL061887

* Carazzo, G., Kaminski, E., & Tait, S. (2006). The route to self-similarity in turbulent jets and plumes. *Journal of Fluid Mechanics*, 547, 137–148. https://doi.org/10.1017/S002211200500683X

* Carazzo, G., Kaminski, E., & Tait, S. (2008). On the rise of turbulent plumes: Quantitative effects of variable entrainment for submarine hydrothermal vents, terrestrial and extra terrestrial explosive volcanism. *Journal of Geophysical Research: Solid Earth*, 113(B9), 1–19. https://doi.org/10.1029/2007JB005458

* Chazette, P., Bocquet, M., Royer, P., Winiarek, V., Raut, J.-C., Labazuy, P., et al. (2012). Eyjafjallajökull ash concentrations derived from both lidar and modeling. *Journal of Geophysical Research: Atmospheres*, 117, D00U14. https://doi.org/10.1029/2011JD015755

## Acknowledgments

We thank two anonymous referees and the editor M. Edmonds for their insightful comments that greatly improved the manuscript. This work and recent developments of the HOTVOLC system were supported by the Centre National d'Etudes Spatiales (CNES-France), and a CNES Postdoctoral Grant. This study contributes to the Service National d'Observations en Volcanologie (SNOV, CNRS-INSU). This is contribution number 721 of the ClerVolc program of the International Research Center for Disaster Sciences and Sustainable Development of the University of Clermont Auvergne.
Cole, P., Barclay, J., Robertson, R., Mitchell, S., Davies, B., Constantinescu, R., et al. (2023). Explosive sequence of La Soufrière, St Vincent, April 2021: Insights into drivers and consequences via eruptive products. In R. Robertson, E. Joseph, J. Barclay, & R. Sparks (Eds.), *The 2020-21 Eruption of La Soufrière Volcano, St Vincent* (Vol. 539(1), pp. 81–106). Geological Society, London, Special Publications. https://doi.org/10.1144/SP539-2022-292

Constantinescu, R., White, J., Connor, C., Cole, P., Fontijn, K., Barclay, J., & Robertson, R. (2023). Estimation of eruption source parameters for the 2021 La Soufrière eruption (St Vincent): Implications for quantification of eruption magnitude on volcanic islands. In R. Robertson, E. Joseph, J. Barclay, & R. Sparks (Eds.), *The 2020-21 Eruption of La Soufrière Volcano, St Vincent* (Vol. 539(1), pp. 107–120). Geological Society, London, Special Publications. https://doi.org/10.1144/SP539-2023-38

Costa, A., Suzuki, Y. J., Cerminara, M., Devenish, B. J., Esposti Ongaro, T., Herzog, M., et al. (2016). Results of the eruptive column model inter-comparison study. *Journal of Volcanology and Geothermal Research*, **326**, 2–25. https://doi.org/10.1016/j.jvolgeores.2016.01.017

Degruyter, W., & Bonadonna, C. (2013). Impact of wind on the condition for column collapse of volcanic plumes. *Earth and Planetary Science Letters*, **377–378**, 218–226. https://doi.org/10.1016/j.epsl.2013.06.041

Devenish, B., Rooney, G., Webster, H., & Thomson, D. (2010). The entrainment rate for buoyant plumes in a crossflow. *Boundary-Layer Meteorology*, **134**(3), 411–439. https://doi.org/10.1007/s10546-009-9464-5

Ernst, G. J., Stephen, R., Sparks, J., Carey, N., & Bursik, M. I. (1996). Sedimentation from turbulent jets and plumes. *Journal of Geophysical Research*, **101**(95), 5575–5589. https://doi.org/10.1029/95jb01900

Girault, F., Carazzo, G., Tait, S., Ferrucci, F., & Kaminski, E. (2014). The effect of total grain-size distribution on the dynamics of turbulent volcanic plumes. *Earth and Planetary Science Letters*, **394**, 124–134. https://doi.org/10.1016/j.epsl.2014.03.021

Girault, F., Carazzo, G., Tait, S., & Kaminski, E. (2016). Combined effects of total grain-size distribution and crosswind on the rise of eruptive volcanic columns. *Journal of Volcanology and Geothermal Research*, **326**, 103–113. https://doi.org/10.1016/j.jvolgeores.2015.11.007

Gouhier, M., Deslandes, M., Guéhenneux, Y., Hereil, P., Cacault, P., & Josse, B. (2020). Operational response to volcanic ash risks using HOTVOLC satellite-based system and MOCAGE-accident model at the toulouse VAAC. *Atmosphere*, **864**, https://doi.org/10.3390/atmos11080864

Gouhier, M., Eychenne, J., Azzaoui, N., Guillin, A., Deslandes, M., Poret, M., et al. (2019). Low efficiency of large volcanic eruptions in transporting very fine ash into the atmosphere. *Scientific Reports*, **9**(1), 1449. https://doi.org/10.1038/s41598-019-38595-7

Gouhier, M., Guéhenneux, Y., Labazuy, P., Cacault, P., Decriem, J., & Rivet, S. (2016). *HOTVOLC: A web-based monitoring system for volcanic hot spots* (Vol. 426, pp. 223–241). Geological Society of London Special Publications. https://doi.org/10.1144/SP426.31

Guéhenneux, Y., & Gouhier, M. (2024). Hotvolc: The official French satellite-based service for operational monitoring and early warning of volcanic ash plumes. *Bulletin of Volcanology*, **86**(4), 29. https://doi.org/10.1007/s00445-024-01716-w

Hedin, A., Fleming, E., Manson, A., Schmidlin, F., Avery, S., Clark, R., et al. (1996). Empirical wind model for the upper, middle and lower atmosphere. *Journal of Atmospheric and Terrestrial Physics*, **58**(13), 1421–1447. https://doi.org/10.1016/0021-9169(95)00122-0

Hewett, T., Fay, J., & Hoult, D. (1971). Laboratory experiments of smokestack plumes in a stable atmosphere. *Atmospheric Environment*, **5**(9), 459–461. https://doi.org/10.1016/0004-6981(71)90028-x

Hoult, D. P., Fay, J. A., & Forney, L. J. (1969). A theory of plume rise compared with field observations. *Journal of the Air Pollution Control Association*, **19**(8), 585–590. https://doi.org/10.1080/00022470.1969.10466526

Ishii, K., Nishijo, A., Koyaguchi, T., & Suzuki, Y. (2022). A physics-based source model for real-time tephra dispersal forecasting for weak eruption plumes. *Journal of Applied Volcanology*, **11**(15), 15. https://doi.org/10.1186/s13617-022-00127-w

Joseph, E., Camejo-Harry, M., Christopher, T., Contreras-Arratia, R., Edwards, S., Graham, O., et al. (2022). Responding to eruptive transitions during the 2020-2021 eruption of La Soufrière volcano, St. Vincent. *Nature Communications*, **13**(1), 4129. https://doi.org/10.1038/s41467-022-31901-4

Kaminski, E., & Jaupart, C. (1998). The size distribution of pyroclasts and the fragmentation sequence in explosive volcanic eruptions. *Journal of Geophysical Research*, **103**(B12), 29759–29779. https://doi.org/10.1029/98jb02795

Kaminski, E., Tait, S., & Carazzo, G. (2005). Turbulent entrainment in jets with arbitrary buoyancy. *Journal of Fluid Mechanics*, **526**, 361–376. https://doi.org/10.1017/s0022112004003209

Mastin, L. (2014). Testing the accuracy of a 1-D volcanic plume model in estimating mass eruption rate. *Journal of Geophysical Research*, **119**(5), 2474–2495. https://doi.org/10.1002/2013JD020604

Mastin, L., Guffanti, M., Servranckx, R., Webley, P., Barsotti, S., Dean, K., et al. (2009). A multidisciplinary effort to assign realistic source parameters to models of volcanic ash-cloud transport and dispersion during eruptions. *Journal of Volcanology and Geothermal Research*, **186**(1–2), 10–21. https://doi.org/10.1016/j.jvolgeores.2009.01.008

Mereu, L., Scollo, S., Garcia, A., Sandri, L., Bonadonna, C., & Marzano, F. (2023). A new radar-based statistical model to quantify mass eruption rate of volcanic plumes. *Geophysical Research Letters*, **50**(7), 1–10. https://doi.org/10.1029/2022GL100596

Michaud-Dubuy, A., Carazzo, G., & Kaminski, E. (2020). Wind entrainment in jets with reversing buoyancy: Implications for volcanic plumes. *Journal of Geophysical Research: Solid Earth*, **125**(10), e2020JB020136. https://doi.org/10.1029/2020JB020136

Michaud-Dubuy, A., Carazzo, G., Kaminski, E., & Girault, F. (2018). A revisit of the role of gas entrapment on the stability conditions of explosive volcanic columns. *Journal of Volcanology and Geothermal Research*, **357**, 349–361. https://doi.org/10.1016/j.jvolgeores.2018.05.005

Michaud-Dubuy, A., & Gouhier, M. (2025). Additional file of "Validation of near real-time retrieval of plume mass eruption rates": The case of the 2021 eruption of La Soufrière, St Vincent. Figshare. https://doi.org/10.6084/m9.fishshare.29209784.v2

Morton, B., Taylor, G., & Turner, J. (1956). Turbulent gravitational convection from maintained and instantaneous sources. *Proceedings of the Royal Society of London. Series A*, **234**, 1–23.

Picquout, A., Lavigne, F., Mei, E., Grancher, D., Noer, C., Vidal, C., & Hadmoko, D. (2013). Air traffic disturbance due to the 2010 merapi volcano eruption. *Journal of Volcanology and Geothermal Research*, **261**, 366–375. https://doi.org/10.1016/j.jvolgeores.2013.04.005

Pouget, S., Bursik, M., Webley, P., Dehn, J., & Pavolonis, M. (2013). Estimation of eruption source parameters from umbrella cloud or downwind plume growth rate. *Journal of Volcanology and Geothermal Research*, **258**, 100–112. https://doi.org/10.1016/j.jvolgeores.2013.04.002

Prata, A. (1989). Infrared radiative transfer calculations for volcanic ash clouds. *Geophysical Research Letters*, **16**(11), 1293–1296. https://doi.org/10.1029/gl016i011p01293

Prata, A., & Prata, A. (2012). Eyjafjallajökull volcanic ash concentrations determined using spin enhanced visible and infrared imager measurements. *Journal of Geophysical Research*, **117**, D00U23.

Prata, F., Prata, A., Tanner, R., Grainger, R., Borgas, M., & Aubry, T. (2024). The radial spreading of volcanic umbrella clouds deduced from satellite measurements. *Volcanica*, **8**(1), 1–29. https://doi.org/10.3090/vol08.01.0129

Robertson, R. (2005). St. Vincent. In J. Lindsay, R. Robertson, J. Shepherd, & S. Ali (Eds.), *Volcanic hazard atlas of the lesser antilles* (pp. 240–261). Seismic Research Unit, The University of the West Indies, Trinidad and Tobago, WI.
Robertson, R., Barclay, J., Joseph, E., & Sparks, R. (2023). An overview of the eruption of La Soufrière Volcano, St Vincent 2020-21. In R. Robertson, E. Joseph, J. Barclay, & R. Sparks (Eds.), *The 2020-21 Eruption of La Soufrière Volcano, St Vincent* (Vol. 539, pp. 1–24). Geological Society, London, Special Publications. https://doi.org/10.1144/SP539-2023-95

Schneider, C., Rasband, W., & Eliceiri, K. (2012). NIH image to ImageJ: 25 years of image analysis. *Nature Methods*, 9(7), 671–675. https://doi.org/10.1038/nmeth.2089

Scollo, S., Bonadonna, C., Garcia, A., Mereu, L., Prestifilippo, M., Romeo, F., et al. (2025). New developments in the estimation of tephra fallout hazard at Mt. Etna, in Italy, during the PANACEA project. *Annals of Geophysics*, 68(1), V102. https://doi.org/10.4401/ag-9174

Shepherd, J., Aspinall, W., Rowley, K., Pereira, J., Sigurdsson, H., Fiske, R., & Tomblin, J. (1979). The eruption of Soufriere volcano, St. Vincent, April-June 1979. *Nature*, 282(5734), 24–28. https://doi.org/10.1038/282024a0

Sigl, M., Winstrup, M., McConnell, J., Welten, K., Plunkett, G., Ludlow, F., et al. (2015). Timing and climate forcing of volcanic eruptions for the past 2,500 years. *Nature*, 523(7562), 543–549. https://doi.org/10.1038/nature14565

Snee, E., Degruyter, W., Bonadonna, C., Scollo, S., Rossi, E., & Freret-Lorgeril, V. (2021). A model for buoyant tephra plumes coupled to lava fountains with an application to the 29th of August 2011 paroxysmal eruption at Mount Etna, Italy. *Journal of Geophysical Research: Solid Earth*, 126(12), e2020JB021360. https://doi.org/10.1029/2020jb021360

Sparks, R., Aspinall, W., Barclay, J., Renfrew, I., Contreras-Arratia, R., & Stewart, R. (2023). Analysis of magma flux and eruption intensity during the 2021 explosive activity at La Soufrière, St Vincent, West Indies. In R. Robertson, E. Joseph, J. Barclay, & R. Sparks (Eds.), *The 2020-21 Eruption of La Soufrière Volcano, St Vincent* (Vol. 539(1), pp. 63–79). Geological Society, London, Special Publications. https://doi.org/10.1144/SP539-2022-286

Sparks, R., Bursik, M., Carey, S., Gilbert, J. S., Glaze, L. S., Sigurdsson, H., & Woods, A. W. (1997). *Volcanic plumes*. John Wiley.

Stewart, C., Damby, D., Horwell, C., Elias, T., Ilyinskaya, E., Tomašek, I., et al. (2022). Volcanic air pollution and human health: Recent advances and future directions. *Bulletin of Volcanology*, 84(1), 11. https://doi.org/10.1007/s00445-021-01513-9

Stinton, A. (2023). Growth and evolution of the lava dome and coulée during the 2020-21 eruption of La Soufrière, St Vincent. In R. Robertson, E. Joseph, J. Barclay, & R. Sparks (Eds.), *The 2020-21 Eruption of La Soufrière Volcano, St Vincent* (Vol. 539, pp. 25–39). Geological Society, London, Special Publications. https://doi.org/10.1144/SP539-2023-304

Suzuki, Y., & Koyaguchi, T. (2009). A three-dimensional numerical simulation of spreading umbrella clouds. *Journal of Geophysical Research*, 114(B03209). https://doi.org/10.1029/2007JB005369

Suzuki, Y., & Koyaguchi, T. (2015). Effects of wind on entrainment efficiency in volcanic plumes. *Journal of Geophysical Research: Solid Earth*, 120(9), 6122–6140. https://doi.org/10.1002/2015JB012208

Suzuki, Y., Koyaguchi, T., Ogawa, M., & Hachisu, I. (2005). A numerical study of turbulent mixing in eruption clouds using a three-dimensional fluid dynamics model. *Journal of Geophysical Research: Solid Earth*, 110(B08201). https://doi.org/10.1029/2004JB003460

Taylor, I., Grainger, R., Prata, A., Proud, S., Mather, T., & Pyle, D. (2023). A satellite chronology of plumes from the April 2021 eruption of La Soufrière, St Vincent. *Atmospheric Chemistry and Physics*, 23(24), 15209–15234. https://doi.org/10.5194/acp-23-15209-2023

Wang, H., & Law, A.-K. (2002). Second-order integral model for a round turbulent buoyant jet. *Journal of Fluid Mechanics*, 459, 397–428. https://doi.org/10.1017/S0022112002008157

Wilson, L. (1976). Explosive volcanic eruptions: III. Plinian eruption columns. *Geophysical Journal International*, 45(3), 543–556. https://doi.org/10.1111/j.1365-246x.1976.tb06909.x

Woods, A. (1988). The fluid dynamics and thermodynamics of eruption columns. *Bulletin of Volcanology*, 50(3), 169–193. https://doi.org/10.1007/bf01079681

Woods, A., & Bower, S. (1995). The decompression of volcanic jets in a crater during explosive eruptions. *Earth and Planetary Science Letters*, 131(3–4), 189–205. https://doi.org/10.1016/0012-821x(95)00012-2

Woods, A., & Bursik, M. (1991). Particle fallout, thermal disequilibrium and volcanic plumes. *Bulletin of Volcanology*, 53(7), 559–570. https://doi.org/10.1007/bf00298156

Woods, A., & Kienle, J. (1994). The dynamics and thermodynamics of volcanic clouds: Theory and observations from the April 15 and April 21, 1990 eruptions of Redoubt volcano, Alaska. *Journal of Volcanology and Geothermal Research*, 62(1–4), 273–299. https://doi.org/10.1016/0377-0273(94)90037-x

Woods, A., & Self, S. (1992). Thermal disequilibrium at the top of volcanic clouds and its effect on estimates of the column height. *Nature*, 355(6361), 628–630. https://doi.org/10.1038/355628a0