Tephra characterization and multi-disciplinary determination of Eruptive Source Parameters of a weak paroxysm at Mount Etna (Italy)

The determination of EruptiveSourceParameters(ESPs)isa majorchallenge especially for weakvolcanicexplo- sionsassociatedwith poorlyexposed tephra-fallout deposits. Insuch a case, the combination of deposit analyses and remote sensing observationscan provide fundamental insights. We use the 29August 2011weak paroxysm at Mount Etna (Italy) as a case study to discuss some of the challenges associated with multi-disciplinary determination of ESPs of poorly exposed tephra-fallout deposits. First, we have determined the erupted mass from a combinationof ﬁ eldandsyntheticdatato ﬁ llasigni ﬁ cantgapindatasampling;syntheticdatahavebeenderived based on extrapolation of ﬁ eld observations and validated based on comparisons with other tephra deposits at Etna and TEPHRA2 modelling. Second, we have combinedthe estimatesof erupted mass and grain-sizedistribu-tion asderived both from deposit observations and satellite retrievals. Analytical modelling was appliedto char- acterize the size fractions most likely represented in satellite retrievals and tephra deposits, respectively. In addition, the Rosin-Rammler distribution ﬁ tting is shown to inform on missing parts of the grain-size distribu- tions and reproduce a tail of very ﬁ ne ash (1 – 20 μ m) whose mass proportion is close to the satellite estimates (1.3 – 1.6% versus 1.9%, respectively). Finally, it was found that this very- ﬁ ne-ash fraction increases as a function of satellite-derived Mass Eruption Rate for a set of eruptions for which independent estimates are available. This critical combination of ﬁ eld observations, analytical modelling and satellite retrievals demonstrates the po- tential and importance of multidisciplinary strategies for the derivation of ESPs even for small-size explosive events and poorly exposed deposits such as that of the 29 August 2011 paroxysm of Mt. Etna. © 2021 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://


Tephra characterization and multi-disciplinary determination of Eruptive Source Parameters of a weak paroxysm at Mount Etna (Italy)
Valentin

Introduction
The determination of Eruptive Source Parameters (ESPs) (e.g., the vent location, Total Erupted Mass -TEM, Mass Eruption Rate -MER, Plume Height -H t , Total Grain-Size Distribution -TGSD) is crucial to understand, assess, and forecast the hazards associated with volcanic eruptions (e.g., Costa et al., 2006;Scollo et al., 2008;Folch, 2012;Bonadonna et al., 2015;Beckett et al., 2020;Aubry et al., 2021). The radius of umbrella clouds can also represent a key ESP of large eruptions (Constantinescu et al., 2021). Typically, ESPs are derived based on the application of dedicated models either to deposit data or to real-time geophysical observations. As a result, the exposure of tephra deposits and/or the availability of accurate geophysical data are of crucial importance to the determination of ESPs. As a matter of fact, more information is available for well-exposed deposits which are typically associated with high-intensity and large-magnitude events than for small events, which are in fact more frequent but more difficult to study (e.g., Connor et al., 2001). In addition, these small events are typically not even preserved in the stratigraphic records as they are easily eroded and mixed with soils or larger events (Gurioli et al., 2013;Kiyosugi et al., 2015). As a result, hazard assessments tend to be biased towards the largest events (e.g., Connor et al., 2001;Bonadonna, 2006). Given the low potential of preservation in the geological record, an effort is required to characterize small events during or soon after the emplacement of tephra deposits in order to have a comprehensive characterization of eruption dynamics at active volcanoes and accurate hazard assessments. Nonetheless, even in case of rapid sampling, the characterization of small explosive events is associated with a variety of challenges (e.g., small accumulations, mix with previous tephra layers, rapid erosion in case of tephra sampling, dilute plumes and clouds in case of remote sensing observations). In addition, a poor sampling coverage may also introduce large biases in the determination of the TEM and TGSD from tephra-fallout deposits, especially when plumes are dispersed off-shore and/or over other inaccessible areas (e.g., Andronico et al., 2014a;Bonadonna et al., 2015). In such a case, the combination of multi-disciplinary strategies might provide important insights.
Multi-disciplinary strategies have been recently developed to compute key ESPs (e.g., TEM, MER, TGSD) based on various methodologies and sensors including the combination of deposit data and geophysical observations Bonadonna et al., 2011;Corradini et al., 2016;Calvari et al., 2018;Poret et al., 2018aPoret et al., , 2018bPioli and Harris, 2019) and the combination of various sensors (Marzano et al., 2018Mereu et al., 2020;Freret-Lorgeril et al., 2021). In particular, the combination of satellite-based observations with tephra-fallout deposit analyses has been suggested to provide additional insights into the TGSD Corradini et al., 2016;Poret et al., 2018aPoret et al., , 2018b. Indeed, satellite data have the potential to provide information on the airborne fraction of tephra, i.e., the so-called very fine ash <30 μm , that is not usually sampled on ground. In fact, even though part of the very fine ash fraction detected by satellite retrievals might fall on the ground at proximal to medial distances from the vent, part of it might remain suspended in the atmosphere for very long times (Durant, 2015).
Here we use the 29 of August 2011 weak paroxysm of Mount Etna (Italy) as a case study to illustrate various multi-disciplinary strategies for the determination of ESPs even in case of poorly exposed deposits and dilute plumes and clouds. At Mount Etna, which is one of the most active volcanoes in Europe, frequent paroxysms have occurred between 2011 and 2015 (Behncke et al., 2014;Scollo et al., 2014;De Beni et al., 2015;Calvari et al., 2018;Corradini et al., 2018;Freret-Lorgeril et al., 2018;Poret et al., 2018a;Andronico et al., 2021). Mount Etna's paroxysms are typically characterized by the emission of lava fountain-fed tephra plumes whose height can reach up to 15 km above sea level (a.s.l.) Calvari et al., 2018;Corradini et al., 2018). Due to the complexity of combining remote sensing and field data and to the fact that remote sensing data are not always available for the studied events, such a combination has been performed only for a few case studies (i.e., Corradini et al., 2016;Poret et al., 2018aPoret et al., , 2018bFreret-Lorgeril et al., 2021).
Using the 29 August 2011 paroxysm as a case study, our work mostly aims at (i) addressing the challenges associated with the determination of key ESPs (TEM and TGSD) of poorly exposed tephra deposits based on dedicated multi-disciplinary strategies and (ii) exploring the assumptions associated with the combination of observations of tephrafallout deposits and satellite data.

Deposit sampling and determination of erupted mass and plume height
A total of 13 samples were collected between 0.4 and 21.7 km from the vent (i.e., New Southeast Crater; NSEC, renamed SEC since 2020) ( Fig. 2A and Table 1). It is important to note that a sampling gap of about 10 km between the proximal samples and the first medial sample exists mostly due to the presence of la Valle del Bove that does not allow access to the deposit (Fig. 2).
The mass of the deposit was computed based on the observed ground accumulation data as a function of the square root of the area inside the isomass contours shown in Fig. 2A using exponential fitting (Pyle, 1989), Power-law  and Weibull integration methods (Bonadonna and Costa, 2012). We also computed an isoMd Φ map (Fig. 2B), with Md Φ calculated based on Inman (1952), in order to correlate the grain-size data with the Plume height (Bonadonna andCosta, 2013a, 2013b). Finally, erupted mass and plume height were also determined based on the inversion of observed values of mass/area with the analytical model TEPHRA2 that relies on the solution of the advection-diffusion-sedimentation equation Connor and Connor, 2006; https://github. com/geoscience-community-codes/tephra2). In particular, TEPHRA2 allows to distinguish the accumulation of several particle size classes thanks to the use of a size selective diffusion law and different particle densities. The inversion employs the downhill simplex method which is a useful optimization technique that allows to find the minimum of a function of N independent variables. During inversion, we sequentially run the forward model along the downhill simplex algorithm with different sets of input parameters and end the simulation whenever the goodness-of-fit falls within the user tolerance.

Whole Deposit Grain-Size Distribution
Individual grain-size distributions, referred to as empirical GSDs hereafter (see Table 1), were obtained by mechanical sieving of all   1.4-2 0.14 10 1-1.4 0.11 11 0.5-1 0.11 samples using ½ Φ intervals from −3.5 Φ (11,300 μm) down to 4 Φ (63 μm). The mass percentage of all sieved fractions was obtained with a 10 −4 g resolution weighing scale. The cumulative ash fraction >4 Φ was analyzed with a BETTERSIZER laser diffraction machine at 0.5 Φ intervals. In fact, the >4 Φ fraction of individual samples was too small to be analyzed individually. It is important to note that only the 7 medial samples contained particles >4 Φ.
Here we distinguish between field-based Whole Deposit Grain-Size Distribution (WDGSD) and satellite-based Grain-Size Distribution (GSD SAT ; see next section). It is important to note that the GSD SAT combines information of both the plume and the horizontal cloud depending on plume inclination and location of plume corner, while the WDGSD only contains information of the sampled deposit. A critical combination of these two distributions is later discussed to provide the best estimate of the Total Grain-Size Distribution, i.e., the size distribution of all particles ejected during the explosive event.
One efficient way to determine the WDGSD of a non-uniform tephra-fallout deposit data set is to use the Voronoi Tessellation method Bonadonna et al., 2015). However, the Voronoi Tessellation strategy cannot reproduce missing parts of the deposit associated with significant sampling gaps, as it is the case of the 29 August 2011 paroxysm (Fig. 2) (e.g., Alfano et al., 2016). To address this issue, we built synthetic GSDs following the strategy developed by Alfano et al. (2016). Such GSDs were obtained by fitting a Weibull function (Bonadonna andCosta, 2012, 2013a) on the variation of observed Md Φ (in Φ scale; Inman, 1952) as a function of the distance from the vent. Hereafter, we refer to Md Φ ⁎ as Md Φ values in metric scale. By extrapolating this fit and using the mean sorting coefficient σ Φ (Inman, 1952) of all samples, we can determine synthetic GSDs at various distances from the vent that cover the sampling gap. Values of Md Φ and σ Φ of synthetic samples were validated based on comparison with other paroxysms at Etna and with the forward modelling of TEPHRA2 .
Finally, we used the following Rosin-Rammler equation to fit grainsize data: where wt% is the weight percentage of particle with a diameter d (μm), and x 0 and l are the length scale and the shape of the distribution, respectively; in particular, x 0 typically corresponds to the diameter at which~63% of the particles are smaller (Vesilind, 1980;. In fact, the Rosin-Rammler fitting has been shown to best reproduce the distribution tail associated with very fine material that is often missing from the deposit (Murow et al., 1980;. To determine the best Rosin-Rammler parameters x 0 and l to fit our data, we used the Curve fitting application of Matlab.

Particle densities
Density of particles of different size was determined between 63 μm and 36 mm. The density of 100 scoria clasts between 16 and 36 mm selected in the proximal samples was determined using the water immersion technique of Houghton and Wilson (1989) at the Laboratory Magmas and Volcanoes (Clermont-Ferrand) (see Bonny, 2012 for more details). Densities of particles between 63 and 2400 μm selected from the most proximal of the medial samples (S3 in Table 1) were determined using water pycnometry at the University of Geneva (Eychenne and Le Pennec, 2012;Freret-Lorgeril et al., 2019). Three pycnometry measurements were made per analyzed fraction to assess the uncertainty. Our measurements of particle density were used to model the mass decay of particles at the base of the volcanic cloud (see section 2.5).

Cloud mass and Grain-Size Distribution from satellite retrievals
We used observations made by the Spinning Enhanced Visible and Infrared Imager (SEVIRI) on board the Meteosat Second Generation (MSG) geostationary satellite to complement the determination of total erupted mass and WDGSD from tephra-fallout deposit analyses. SEVIRI is a multispectral instrument with 12 channels from visible to thermal infrared, a temporal resolution that ranges from 5 min (Rapid Scan Mode) to 15 min (Earth full disk) and a nadir spatial resolution of 3 km at sub satellite point. The ash particle retrievals have been performed using the Volcanic Plume Retrieval procedure Guerrieri et al., 2015) considering the SEVIRI thermal infrared channels centered at 10.8 and 12 μm and applied to SEVIRI images from 04:15 to 11:00 UTC (resolution of 15 min). For each volcanic cloud pixel, the Volcanic Plume Retrieval procedure computes simultaneously the particles effective radius (micron), Aerosol Optical Depth (at 550 nm) and total column abundance-M a (kg/m 2 ). An ash mass flux F(t) (kg/s) was computed from each 15 min SEVIRI image considering only pixels at a distance of 15 ± 1.5 km from the summit craters (nearest pixels to the crater usually have large opacity that may lead to a large uncertainty in the retrievals). The flux was then computed using the following equation: where l trans is the transect width (in this case 3 km), v w is the wind speed (m/s). Finally, by means of the wind speed, the time t of the flux at 15 km was then reported to the crater (0 km). The average mass flux multiplied by the whole eruption time gives the total erupted mass (kg). Finally, the GSD SAT and mass accumulation of particles belonging to the 1-20 μm fraction (10-5.5 Φ) were obtained considering the whole cloud extent and combining all the images together (Wen and Rose, 1994;Prata and Grant, 2001;Stevenson et al., 2015;Corradini et al., 2016Corradini et al., , 2018Poret et al., 2018aPoret et al., , 2018bGouhier et al., 2019). It is important to keep in mind that particles coarser than 5.5 Φ, and potentially present in the volcanic cloud, are not discriminable in the thermal infrared spectral range. The reason is that the radiative effect of these particles is approximately the same as particles of 5.5 Φ. This could lead to an underestimation of the volcanic cloud total ash mass retrieved if only satellite measurements are used (see Stevenson et al., 2015, and references therein).
Following the methodology presented in Freret-Lorgeril et al. (2021), the GSD SAT was computed in wt% using the ratio between the total mass of pixels containing ash particles at each time step, within a given Φ range, and the total mass of the detected volcanic plume/ cloud. Finally, we averaged all individual distributions from images containing a minimum of 100 pixels with ash signal to derive the mean GSD SAT of the whole event.

Mass/area decay of individual size fractions at the base of the volcanic cloud
With the aim of better understanding the relation between the tephra fraction deposited on the ground and that detected by satellite retrievals, we used the 1D integral model of Bonadonna and Phillips (2003) to estimate the mass decay of each size fraction in the umbrella cloud with distance from vent. The variation of mass per unit area (kg/ m 2 ) for each particle size at the base of the umbrella cloud is calculated along the main dispersal axis starting from the plume corner (calculated based on the eq. 12b of Bonadonna and Phillips (2003)). Particles were considered spherical, their density was determined in the lab (section 2.3) and their terminal velocity v t (m/s) was computed based on the drag equation of Bonadonna (2016a, 2016b). We used meteorological data from the hydrometeorological service of ARPA in Emilia Romagna Scollo et al., 2019) to derive wind velocity. It is important to note that results associated with 1D modelling are associated with some uncertainties even though they provide a first approximation estimate of tephra accumulation both at the base of umbrella clouds and on the ground (Bursik et al. 1992a(Bursik et al. , 1992bSparks et al., 1992;Bonadonna and Phillips, 2003). In addition to the description of mass/area decay at the base of the volcanic cloud, we also inverted the tephra-fallout deposit to investigate the most characteristic position of release of individual size fractions. The height of the cloud base H cb of 5.4 km (a.s.l.) was determined from the Plume Top Height H t of 6.7 km above vent (i.e., 9.6 km a.s.l.; Corradini et al., 2018) based on the Eqs. (1) and (2) of Bonadonna and Phillips (2003).

Total erupted mass
Fig . 3A shows the ash mass as a function of time as estimated from the different SEVIRI images. The TEM SAT visible by the satellite images for the entire period (i.e., satellite-based erupted mass) has been estimated to 2.5 × 10 6 kg. Interestingly, the ash mass of the imaged volcanic cloud does not vary significantly after the end of the paroxysm (04:53 UTC; vertical dashed line in Fig. 3A). The cloud is tracked up tõ 400 km from the vent until 11:00 UTC.
The TEM DEP (i.e., erupted mass as derived from the tephra-fallout deposit) based on the integration of the exponential, Weibull and Power-Law functions fitted on the variation of ground accumulation with the square root of the area of isomass contours (Fig. 3B;Pyle, 1989;Bonadonna and Costa, 2012) is 1.3 × 10 8 kg, 1.4 × 10 8 kg, and 1.4 ± 0.1 × 10 8 kg, respectively. The average TEM DEP considering all aforementioned values is 1.4 ± 0.1 × 10 8 kg. The proximal integration limit for the power-law was set as established by the strategy proposed by . The distal limit is varied between 100 km and 300 km, i.e., the distances at which the ground accumulation becomes negligible (10 −3 -10 −4 kg/m 2 ; Fig. 3B). In any case, given that the power-law exponent is >2, the volume estimates are not sensitive to the distal integration limit (e.g., Bonadonna et al., 2015). It is important to note that the exponential and Weibull best-fits are integrated to infinity, and, therefore, Fig. 3. A) Variation in time of plume/cloud mass (kg) as retrieved from satellite images. Vertical dashed line at 04:53 UTC represents the end of the paroxysm based on radar retrievals (Freret-Lorgeril et al., 2018). The arrow indicates the time at which the total mass starts to decrease. B) Variation of ground accumulation as a function of the square root of the isomass contour areas and associated Exponential (blue), Power-Law (red line; y = 19.85 × x -2.11 ) and Weibull (Purple dashed line) best-fits. C) Variation of ground accumulation as a function of the square root of the isomass contour areas for tephra-fallout deposits produced by small and large Etna paroxysms, Basaltic Plinian and Subplinian eruptions and other Plinian eruptions. References: Cotopaxi L3 and L5 (Biass and Bonadonna, 2011); Pululagua (Volentik et al., 2010); Taupo (Walker, 1980); Fontana Lapilli (Costantini et al., 2009);Fuego 1974(Rose et al., 2008; Cerro Negro 1992 (Connor and Connor, 2006); Tarawera 1886 (Walker et al., 1984); Etna 122 BC (Coltelli et al., 1998) they include all deposited material but not what remained suspended in the atmosphere. Moreover, the proportion of material that sedimented on land was about 7.1 × 10 7 kg based on the integration of the powerlaw fitting with a distal integration limit of 12 km, i.e., the approximate distance between the vent and the coastline along the main plume dispersal axis converted in square root of area (Figs. 2A and 3B). Using an event duration of 60-70 min (Calvari et al., 2018;Freret-Lorgeril et al., 2018), the average MER based on the TEM DEP is equal to 3.5 ± 0.3 × 10 4 kg/s. Regardless of the poor deposit exposure, the thinning trend of the 29 August 2011 paroxysm (red squares in Fig. 3C) is in good agreement with the trend of other small events at Etna; associated values of MER are also in agreement (between 10 2 and 10 4 kg/s) (see also Appendix A). Note that values of MER of strong paroxysms have been estimated between 10 5 and 10 6 kg/s and such events clearly display higher ground accumulation than those obtained for small paroxysms (Fig. 3C). The inversion of tephra ground accumulation using the model TEPHRA2 Connor and Connor, 2006) gives an erupted mass of 1.8 ± 0.1 × 10 8 kg (see Table B1 in Appendix B). Based on the inversion of the tephra deposit, we used TEPHRA2 to compute an isomass map (i.e., forward modelling; Fig. 4A) and the ground accumulation at each sampling site (Fig. 4B) which provide very similar results than the observed data (e.g., Table 1). In particular, on 13 computed datapoints, 12 of them lie within ±50% from the perfect match (1:1 line) with the observed samples.

Determination of plume height
Based on the ground accumulation data, we can find an inversion solution with TEPHRA2 that provides a good agreement with satellite observations of plume height, i.e., 11.1 ± 1.4 km a.s.l. against 9.0-9.6 km a. s.l., respectively. Moreover, plume height can also be derived based on Md Φ ⁎ values converted in cm (Bonadonna andCosta, 2013a, 2013b). As an example,   In "Other Plinian eruptions" Cotopaxi L3 is andesitic, Cotopaxi L5 is basaltic-andesite and Pululagua 2450BP is dacitic; figure modified from Bonadonna andCosta (2013a, 2013b). B) Variation of plume heights above vent level (a.v.l. km) as a function of the Weibull parameter λ MdΦ (Bonadonna andCosta, 2013a, 2013b). References are in Bonadonna andCosta (2013a, 2013b)  values between 23 and 50 km. The data we obtained for the 23 November 2013, 29 August 2011, 10 April 2011 and 12-13 January 2011 paroxysms of Etna extend this trend to lower heights (7.8 km, 6.4 km, 4.3 km and 7.1 km above vent, respectively) with a new equation that presents a very good coefficient R 2 of 0.94 (Fig. 5B):

Whole Deposit Grain-Size Distribution
The application of the Voronoi Tessellation to the 13 collected samples (Table 1) results in a bimodal WDGSD with a gap of sizes between −2 and 1 Φ centered on −1 Φ (i.e., white histogram in Fig. 6A). We followed the strategy of Alfano et al. (2016) to interpolate the proximal samples with the medial to distal samples (S1 to S7) based on the variation of Md Φ values with distance from vent fitted by a Weibull distribution (Fig. 6B). Five equally spaced synthetic samples were created assigning the corresponding Md Φ ⁎ value based on the Weibull fit (blue dots in Fig. 6B) and an average sorting σ Φ value (i.e., 0.7; Table 2). The resulting WDGSD obtained applying the Voronoi Tessellation on both the collected samples and synthetic data is well sorted with a Md Φ value of −0.9 (1866 μm) and a σ Φ of 2.1 (blue histogram in Fig. 6A).
Proximal samples display Md Φ values between −3.3 and − 3.9 with σ Φ between 0.7 and 0.9, in agreement with values observed for 8 Etna paroxysms (Table 1, Fig. 6C-D). Moreover, values of Md Φ and σ Φ are also in agreement with modelled values derived from the forward simulations of TEPHRA2 using the input parameters obtained from the inversion (Fig. B1 in Appendix B).  Table 2 Grain-size parameters based on mechanical sieving. Md Φ and σ Φ are from Inman (1952) and x 0 indicates the Rosin-Rammler parameter below which 63% of the grain-size distribution is smaller (see Eq. (1)). name Md

Particle density
Scoria clasts of size between 16 and 36 mm taken from proximal samples (Fig. 2) have densities between 279.3 and 1537.7 kg/m 3 with an average value of 555.7 ± 149.4 kg/m 3 (Fig. 7). The densities of ash and small lapilli clasts from 63 to 2400 μm are between 703.1 ± 224.6 kg/m 3 and 2650 ± 472.8 kg/m 3 , respectively. The variation of particle densities from fine to coarse-grained samples has a sigmoidal shape, similar to observations made on samples from the 2006 subplinian eruption of Tungurahua volcano (Ecuador) by Eychenne and Le Pennec (2012). Herein, we determine a relationship between the particle size (d in m) and densities ρ p (in kg/m 3 ) in the following form: with an excellent R 2 of 0.96. This law is used hereafter to determine the densities of all particles from 10 to 6 Φ in the 1D model of Bonadonna and Phillips (2003) (see Sections 3.5 and 3.6).

WDGSD, GSD SAT and TGSD
As already mentioned, the TEM SAT is derived from the mass flux computed at 15 km from the vent. This value is associated with particles theoretically comprised between 1 and 20 μm (i.e., > 5.5 Φ). Our 1D simulations based on the model of Bonadonna and Phillips (2003) and using the particle density calculated in this work, show that, unless falling as part of aggregates and/or gravitational instabilities, particles <20 μm would sediment starting from about 350 km from vent (red arrow in Fig. 8). In addition, the main tephra fraction deposited and sampled on the ground was likely released from the base of the umbrella cloud at distances <15 km (green arrows in Fig. 8). This shows that, for the largest part, the ash fraction detected by satellite retrievals at the plume/cloud top is not included in the analysis of WDGSD (see sampling distances in Table 1). The TGSD can, therefore, be obtained by combining WDGSD and GSD SAT .
The empirical WDGSD discussed in previous section was fitted by a Rosin-Rammler distribution (red line in Fig. 9A) with best parameters x 0 = 3474 μm and l = 0.837 (see Eq. (1)), a correlation coefficient R 2 of 0.99 and a Root Mean Square Error (RMSE) of 0.03. This Rosin-Rammler best fit has a Md Φ of −1.2 (2235 μm) and is well sorted with a σ Φ of 2.0.
Over the whole duration of the event, the average GSD SAT derived from satellite that only considers particle between 10 and 5.5 Φ shows a Md Φ value of 7.2 Φ and a σ Φ of 0.8 (Fig. 9B). While mass computation methods from deposit also contain particles that sedimented beyond the coastline, the determination of the WDGSD is limited to the analyzed samples and does not extrapolate to nonsampled areas, e.g., distal zones covered by satellite images. Hence, as mentioned above, a sum of both the WDGSD and GSD SAT , considering the ratio between the field and satellite total masses (i.e., 1.9%; see Section 3.1), needs to be done in order to best approximate the TGSD (black histogram in Fig. 9A).
A significant result is the similarity that exists between the wt% of the 1-20 μm fraction (i.e., the fraction detected by satellite) observed in the TGSD and in the Rosin-Rammler fitted sieve-based WDGSD and TGSD (Fig. 9). Indeed, the 1-20 μm fraction represents 1.9 wt% of the TGSD, 1.6 wt% of the Rosin-Rammler fitted on the TGSD (Fig. 9A). Moreover, this very fine fraction represents 1.3 wt% of the sieve-based Rosin-Rammler WDGSD. In contrast, the 1-20 μm fraction of the sieve-based WDGSD not fitted by Rosin-Rammler is 0.03 wt%.
The combined TGSD shows a clear gap in between 3 and 6 Φ. This gap might represent a lack of fine material sampled on ground because of the limited deposit exposure, i.e., presence of the coastline (Figs. 1B and 3), and the physical limitation of the retrievals in the thermal infrared spectral range to discriminate particles coarser than 20 μm. The Rosin-Rammler best fits suggest that this gap represents 4.6-5.2% of missing material in our dataset.
Finally, 1.5% of the total 5.5-10 Φ fraction (combined ground and satellite data) sedimented within 25 km from vent even though the characteristic sedimentation distance of this fraction is beyond 350 km (Fig. 8). Even though 1.5% is significantly lower with respect to the 46% observed for the same fraction for the Eyjafjallajökull 2010 eruption (up to 56 km from the vent; Bonadonna et al., 2011), we conclude that this fraction has been probably sedimented as part of aggregates and/or gravitational instabilities. Aggregates were not observed in the deposit, but they could have been broken with impact on the ground.

Mass/area decay at the base of the volcanic cloud
The mass/area decay of various particle sizes at the base of the volcanic cloud is determined along the dispersal axis based on both satellite and 1D model estimates (Fig. 10A). The mass accumulation (kg/m 2 ) from satellite was obtained considering the maximum value of each transect perpendicular to the dispersal cloud axis and by considering all the SEVIRI images from the start of the eruption until 11:00 UTC. The modelled decay is computed with the model of Bonadonna and Fig. 7. Variation of clast density ρ p (black dots) in kg/m 3 as a function of particle size (Φ unit). Red dots are the average values obtained for each analyzed fraction (± standard deviation). The blue line corresponds to the best fit between density and particle size (see Eq. (4)) with a 95% confidence interval (dashed blue line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Phillips (2003) based on the assumption that the TGSD as calculated in previous sections is associated with the average TEM as derived from tephra-fallout deposit analysis (1.4 × 10 8 kg) combined with the satellite-derived mass (2.5 × 10 6 kg). It is interesting to notice that the modelled accumulation of the 5.5-10 Φ fraction is twice to three times higher than the same fraction as retrieved from satellite data but presents a similar decay. The overestimation of the modelled fraction could be related to the fact that the 5.5-10 Φ fraction sedimented closer to the vent than expected due to aggregation and/or gravitational instabilities. Such processes are not taken into account in our 1D model. In addition, our 1D model also shows that at 15 km from vent the 3-6 Φ fraction is associated with significant values of mass/area (>10 −1 kg/ m 2 ) even though it is not detected by the satellite retrieval algorithm (Fig. 10B). This might explain the gap between 3 and 6 Φ shown in Fig. 9A. Finally, it is also interesting to notice that around 450 km from vent, when the cloud is very diluted, all the size fractions detectable by the satellite retrievals (5.5-10 Φ) are below 10 −2 kg/m 2 (Fig. 10C).

Plume height determination
Even though the distribution of the largest clasts around a volcano is typically used to determine the maximum plume height (e.g., Carey & Sparks, 1986;Rossi et al., 2019), such strategy is scarcely used at Etna due to the typical lack of proximal data and the small size of distal tephra clasts (e.g., Coltelli et al., 1998). Here we have shown how the variation of Md Φ as a function of the square root of isoMd Φ areas (Figs. 2B and 5) can also be used at Etna to retrieve an average plume height of the event based on the value of λ MdΦ (using a Weibull function; Bonadonna and Costa, 2012, 2013a, 2013b. In particular, with this study we have extended the empirical relation proposed by Bonadonna andCosta (2013a, 2013b) to average H t values from 23 km down to about 4 km above vent (i.e., height of 10 April 2011 paroxysm; Freret-Lorgeril et al., 2021). We expect such a trend to help characterize H t from tephra-fallout deposits of future eruptions, especially in case of missing remote sensing data. We also found that, regardless of the poor exposure of the deposit, the inversion with TEPHRA2 of ground accumulation data provides consistent results for plume height with satellite observations (i.e., between 9 and 11 km a.s.l.; see Appendix B). The height differences between TEPHRA2 and satellite estimates shown in Appendix B are not significant in comparison to the typical uncertainties observed by satellite and visible imagery, i.e., of about ±0.5 km Scollo et al., 2019 and references therein). This agreement also confirms that the mass/area values used for TEPHRA2 inversion are related to the spreading of the cloud that was observed in the atmosphere for about 6-7 h (based on satellite images) with a duration of the associated paroxysm of a little more than one hour. In fact, most samples were collected beyond the plume corner (around 1-2 km), being the height above the vent around 6-8 km and the plume being strong (Calvari et al., 2018).

Insights into TGSDs from the Rosin-Rammler fitting
Based on a large data set,  have found the Rosin-Rammler formulation (Eq. (1)) as the best strategy to reconstruct the grain-size distribution tails of tephra-fallout deposits. This is particularly interesting given that WDGSDs often lack proximal and/or distal data due to erosion and/or difficult accessibility. Accordingly, using the best x 0 and l parameters (e.g., with Matlab), we were able to reproduce the 29 August 2011 empirical WDGSD with an excellent R 2 >0.99 (with Md Φ value of −1.0; Fig. 6A). We found that the wt% of the 5.5-10 Φ very fine ash fraction well correlates with the WDGSD and the TGSD (Fig. 9A). This suggests the capacity of the Rosin-Rammler equation to reproduce the tails of WDGSDs  that are rarely sampled on ground and to provide a first order estimate of the airborne fraction potentially detected by satellite.
Both Rosin-Rammler parameters, x 0 and l, associated with the empirical WDGSD and the TGSD (resulting from the combination of the WDGSD and the GSD SAT ) are in the same range with respect to those found for 49 tephra deposits at other volcanoes analyzed by   (Fig. 11A). They also show similar trends of x 0 versus Md Φ ⁎ (m) and l versus σ Φ (Fig. 11B,C). In particular, whereas x 0 intuitively increases as a function of Md Φ ⁎ (m) (Fig. 11B), the better sorted the WDGSDs, the smaller the parameter l (Fig. 11C). In particular, the Rosin-Rammler WDGSD of the 29 August 2011 event presents close l value (0.84) with respect to the 12-13 January SEC paroxysm (l = 0.94) and the 27 October 2002 South-East Crater strong eruption (l = 0.98), for which proximal data (i.e., <5 km from the vent) were also acquired (Andronico et al., 2008(Andronico et al., , 2014a. Even if other data obtained at Etna present highly variable values of l (blue squares in Fig. 11A-C), it is important to bear in mind that extreme l values are not necessarily associated with poor deposit exposure. As an example, Fig. 9. A) Empirical WDGSD (blue histograms, equivalent to the blue histograms in Fig. 7A) and associated Rosin-Rammler best fit (red line). The empirical WDGSD was combined with GSD SAT from satellite retrievals to obtain the best approximation of TGSD (black histogram); the TGSD Rosin-Rammler best fit is also shown with the purple line. The gray area corresponds to the 10-5.5 Φ area detected by satellite retrievals. B) Satellite-based (computed) average GSD SAT used to derive the TGSD in Fig. 9A (black histograms). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) the highest l of 2.47 (Fig. 11A-C) was derived from the 24 November 2006 deposit that was sampled from proximal (1 km) to distal (80 km) areas at 27 sampling sites (Andronico et al., 2014b).
Various attempts to correlate grain-size parameters (i.e., fine ash content, Md Φ , σ Φ ) with eruptive style and ESPs exist (e.g., Rust and Cashman, 2011;Costa et al., 2016;. The percentage of tephra <1 mm, measured along the dispersal axis where the isopach is 10% of the maximum deposit thickness, was also proposed by Walker (1973) as an indicator of explosiveness of an eruption (i.e., the fragmentation index) used to classify explosive eruptions together with the area of pyroclastic dispersal. Nonetheless, all these correlations and classifications of eruptive styles based on grain-size parameters are affected by the large uncertainties associated with both the determination of GSDs and WDGSDs and with the determination of ESPs (e.g., MER, plume height); therefore, they often only provide limited and misleading insights. Similarly, the overall 1-20 μm fractions as derived from the Rosin-Rammler best fits of the aforementioned WDGSDs/TGSDs present a poor correlation with MER ( Fig. 12A; Pioli et al., 2019 and references therein). Nonetheless, a better correlation is shown if data are grouped in 3 categories: 1) tephra-fallout deposits older than 1000 years (i.e., most likely depleted in very fine ash due to erosion and/or poor deposit exposure; open squares), 2) recent events rich in very fine ash due to specific eruptive styles and dynamics (e.g., events associated with PDCs), collective sedimentation processes (e.g., aggregation and/or gravitational instabilities) and/or efficient sampling strategy (i.e., sampling carried out during or just after deposit emplacement) and/or good deposit exposure (i.e., large part of the deposit land and accessible for sampling) (green triangles) and 3) all other events of  dataset including Etna's eruptions (black and light blue circles, respectively) (Fig. 12A). The outlier with very small content in very fine ash (black dot) corresponds to the 1986 coarse-grained basaltic-andesite eruption of Izu-Oshima volcano that forms a small island in Sagami Bay (Japan), and, for which, only the proximal deposit could be analyzed (Mannen, 2006). Moreover, the content of very fine ash increases as a function of satellite-based MERs ( Fig. 12B; data from Gouhier et al., 2019). These results suggest that, in case of good deposit exposure and sampling, the content of very fine ash should increase as a function of eruption intensities, with some variation due to specific eruptive dynamics (e.g., presence of co-PDC plumes). It is important to also note that some of the data scatter in Fig. 12A is partly due to the various strategies used in literature to determine MER either from plume height (e.g., Wilson and Walker, 1987;Mastin et al., 2009;Degruyter and Bonadonna, 2012) or from the ratio between deposit mass and eruption duration.  Bonadonna and Phillips (2003) for the 5.5-10 Φ fraction (blue line). Modelled mass accumulation at the cloud base of particle sizes between −3 and 10 Φ up to 25 (B) and 1000 km (C) from the vent. The vertical dashed line in B indicates the distance of 15 km from the vent at which satellite mass estimates are computed. The black vertical line in B indicates the plume corner location (=0.2H t ; Bonadonna and Phillips, 2003). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Multidisciplinary determination of ESPs
The application of the three empirical fitting methods for the integration of EM from tephra-fallout deposits (i.e., Exponential, Power-Law and Weibull) provides similar results due to the relatively good exposure of the deposit in relation to a relatively small paroxysm. In fact, the exponent of the Power-Law fit is >2 (Fig. 3B), suggesting both a small magnitude and limited dispersal of the associated eruptive event and a low impact of distal extreme of integration on the erupted mass estimates (Bonadonna et al., 2015). As a result, we consider the average deposit TEM of 1.4 ± 0.2 × 10 8 kg derived from the Exponential, the Power-Law (with integration limits of 100 and 300 km) and Weibull empirical methods (Fig. 3B) as the best estimate. It should be noted that all these values of TEM are associated with various levels of uncertainty which originate from sampling of mass/area (up to 30%), data contouring, availability of samples and the deposit exposure (up to 70%; e.g., Bonadonna et al., 2015). Nonetheless, this estimate is still in agreement with results from the inversion of ground accumulation data with TEPHRA2 (i.e., 1.8 ± 0.1 × 10 8 kg, see Appendix B) and compares favorably to 4.4 × 10 8 kg obtained from L-band radar retrievals (Freret-Lorgeril et al., 2018). It is important to bear in mind that, the Exponential and Weibull fits are integrated to infinity, whereas the Power-Law is ideally integrated to the most distal extension of the tephra-fallout deposit. Hence, as stated earlier, our TEM DEP accounts both for the material that sedimented on land (i.e.,~50% of the TEM DEP ) and for the material that sedimented beyond the coastline (i.e., >25 km), including the material detected by satellite that eventually has fallen either on the ground or in the sea (Fig. 1B).
TEM SAT and GSD SAT are likely to represent both airborne and deposited very fine ash (either on land or out at sea), even though it is difficult for the two fractions to be quantified. The mass accumulation of the 5.5-10 Φ fraction detected by SEVIRI is shown to decrease as a function of the distance from vent (Fig. 10A). Our estimates show that this very fine ash, if falling as individual particles, would eventually reach the ground at distances beyond 350 km from the vent (Fig. 8) and is, in any case, considered in the deposit-based mass estimates (i.e., TEM DEP ) (Fig. 3B). The ratio between TEM SAT and TEM DEP is 1.8%, 1.9% and 2.0% depending on the integration based on the Weibull, Power-law, and Exponential fit, respectively, with a mean of 1.9 ± 0.1%. This means that at 15 km, 1.9 ± 0.1% of the 5.5-10 Φ fraction is still in the volcanic cloud. Interestingly, this ratio with the TEM DEP computed with TEPHRA2 is similar, with a value of 1.4%. Nonetheless, we cannot exclude that some of the ash retrieved in the satellite images might never sediment. In such a case, instead of using the ratio between both TEM SAT and TEM DEP to determine the proportion of very fine ash that is detected by satellite (see section 3.2), TEM SAT should be divided by TEM SAT + TEM DEP . In the case of the 29 August 2011 paroxysm, the two strategies result in a similar proportion of 1.8-1.9%, respectively. In addition, it is important to note that particles coarser than 5 Φ are likely to remain in the volcanic cloud beyond 15 km as shown by our 1D model (accumulation >10 −1 kg/m 2 at the base of the volcanic cloud up to at least 25 km from vent; Fig. 10B and C). If these coarse-ash particles contribute to the cloud thermal signature, they could not be discriminated from 20 μm particles due to Mie Scattering effects (Prata and Grant, 2001;Stevenson et al., 2015;Gouhier et al., 2019). In such a case, TEM SAT should be considered as minimal values.
Differently with respect to the TEM estimates, the very fine ash fraction detected by satellite retrievals needs to be added to the WDGSD derived from tephra-fallout deposits. In fact, all current strategies used to determine the WDGSD (including the Voronoi Tessellation) only consider the sampled tephra-fallout deposit and do not extrapolate grainsize data to infinity or to the distal extent of the deposit. Our 1D model simulations show that the >5.5 Φ fraction would deposit in the sea (Fig. 8), and, therefore, if it did not fall closer to the vent as aggregates and/or gravitational instabilities, it is not accounted for in the WDGSD. Hence, WDGSD and GSD SAT can be combined based on the relative mass proportion.
It is important to note that our combined satellite-field TEM estimate characterizes the fraction of tephra that is transported in the atmosphere by the buoyant plume and umbrella cloud and deposited away from the eruptive vent (i.e., the SEC). Nonetheless, in the case of lava fountain-fed tephra plume events such as Mt. Etna's paroxysms, lava fountains generate an additional very proximal deposit that builds the cone (Behncke et al., 2014;Andronico et al., 2014a;De Beni et al., 2015;Spanu et al., 2016;Freret-Lorgeril et al., 2018). Such a very proximal tephra-fallout deposit associated with the lava fountain should be further investigated to better understand the dynamics of lava fountainfed tephra plumes generated during paroxysmal events at Etna and other volcanoes (e.g., Snee et al., 2021).
Finally, based on the comparison with observations from other Etna paroxysms and modelled Md Φ values (Figs. 3, 6 and Appendix B), we have shown that using synthetic data to fill sampling gaps in order to determine the WDGSD not only is applicable to large deposits such as that of the 2008-2013 Chaitén eruption (Alfano et al., 2016) but also to small deposits such as that of the 29 August 2011 paroxysm of Etna.

Conclusions
Quantifying ESPs at very active volcanoes such as Mount Etna is crucial for both real-time ash dispersal forecasting and long-term hazard assessment. The determination of ESPs of weak explosive events associated with poorly exposed tephra-fallout deposits has been proven to be especially challenging. Nonetheless, the characterization of weak explosive events is crucial to an accurate assessment of the whole range of eruption intensity associated with a given volcano. In this work, we combined particle size analyses performed on a sparsely sampled tephra-fallout deposit with numerical modelling and satellite-based observations of a weak paroxysm of Mount Etna (i.e., the 29 of August 2011 event) to investigate the potential of multidisciplinary strategies for the determination of ESPs.
In terms of characterization of small, poorly exposed tephra-fallout deposits we have shown how: 1) all empirical integrations used and the inversion of ground accumulation data with TEPHRA2 provide similar estimates of TEM DEP (1.4 ± 0.0 × 10 8 kg vs 1.8 ± 0.1 × 10 8 ), which could be due to the limited size of the tephra-fallout deposit; 2) plume height of small paroxysms that are not associated with large clasts (especially large lithic clasts) can also be derived from tephrafallout deposits based on inversion of ground accumulation data as well as on the variation of Md Φ values with distance from vent; 3) poorly-exposed tephra-fallout deposits associated with large sampling gaps related to critical sedimentation regimes should be carefully treated for the determination of WDGSD. Here we show that synthetic data, validated by both the comparison with observed data from other Etna paroxysms and simulations with TEPHRA2, can be used to reconstruct WDGSD that would otherwise result in a bimodal distribution. The resulting WDGSD obtained applying the Voronoi Tessellation on both the collected samples and the synthetic data is well sorted with a Md Φ value of −0.9 and a σ Φ of 2.1. The Rosin-Rammler fitting of such a WDGSD displays similar distribution parameters with respect to other well-sampled WDSGDs at Etna; In terms of multi-disciplinary determination of ESPs, we have shown how: 4) the variation of particle density with particle size shows a sigmoidal distribution that can be used to model particle density in analytical and numerical models; 5) the ash fraction detected by satellite retrievals is not included in the analysis of WDGSD. In fact, unless falling as part of aggregates and/or gravitational instabilities, particles <20 μm would sediment starting from about 350 km from vent, while the main tephra fraction deposited and sampled on the ground was released from the base of the umbrella cloud at distances <15 km where satellite data are retrieved for grain-size. As a result, in order to obtain a TGSD, WDGSD and GSD SAT should be combined considering their relative mass proportion; 6) the Rosin-Rammler fitting already shown to best reproduce the tails of WDGSDs , is here shown to also best reproduce the tail of the TGSD. In fact, the very fine ash fraction detected by satellite retrievals (1-20 μm) observed in the TGSD (1.9-1.2 wt%) is similar to that observed in the Rosin-Rammler fitted sieve-based WDGSD. In contrast, the 1-20 μm fraction of the sieve-based WDGSD not fitted by Rosin-Rammler is 0.03 wt%; 7) the Rosin-Rammler fitting was also used to reconstruct the gap of 3-6 Φ observed in the TGSD suggesting that this fraction is associated with 3.8-5.3 wt% even though it is not represented either in the GSD SAT or in the WDGSD. In fact, our 1D modelling shows that at 15 km where satellite data for mass flux are retrieved, the 3-6 Φ fraction is associated with significant values of mass/area (>10 −1 kg/m2) even though it is not detected by the satellite retrieval algorithm; Additionally, such a fraction is not represented in the WDGSD as it is lost at sea; 8) the ratio between TEM SAT and TEM DEP is especially important for the determination of TGSD; in case of the 29 August 2011 paroxysm, TEM SAT represents 1.9% ± 0.1% of the total deposited mass (i.e., TEM DEP ) and 1.4% of the computed TEM DEP with TEPHRA2. It is important to note that TEM DEP accounts both for the material that sedimented on land (i.e.,~50% of the TEM DEP in case of the 29 August 2011 event) and for the material that sedimented beyond the coastline (i.e., >25 km), including the material detected by satellite that eventually has fallen either on the ground or in the sea. So, differently with respect to the determination of TGSD (point 4), the TEM SAT should be considered as mostly contained in TEM DEP even though we cannot exclude that some of the airborne material never sediments and, therefore, is not included in TEM DEP ; 9) 1.5 wt% of the 5.5-10 Φ fraction of the whole TGSD fell on land (within 25 km from vent) and was observed in the deposit even though the characteristic sedimentation distance of this fraction is beyond 350 km; interestingly, our modelled accumulation of the 5.5-10 Φ fraction is twice to three times higher than the same fraction as retrieved from satellite data but presents a similar decay supporting the possibility that this small fraction was sedimented as part of aggregates and/or gravitational instabilities; 10) at 400 km from vent, when the cloud is very diluted, the modelled mass/area at the base of the volcanic cloud of all the size fractions detectable by the satellite retrievals (5.5-10 Φ) are below 10 −2 kg/m 2 indicating that most of the mass sedimented at distances <400 km as suggested by SEVIRI. 11) the content in very fine ash as derived from the Rosin-Rammler fits of various WDGSDs/TGSDs tends to increase as a function of both plume-height-derived MER and satellitederived MER. This suggests an important correlation between fine-ash content and ESPs often attempted but complicated by the large associated uncertainties and that should be further investigated based on more accurate estimates of WDGSDs/TGSDs (i.e., associated with recent and rapidly sampled tephra-fallout deposits and multi-sensor strategies such as the one presented in this paper).

Author contribution statement
FD, LuG and SS were in charge of collecting the samples and of their preliminary analyses. SC and LoG have provided the results of the SEVIRI retrievals. JL has performed the inversion and forward modelling with TEPHRA2. ER has performed the 1D model simulations. VFL has carried out all WDGSD, TGSD, TEM and particles density estimates. VFL and CB have provided a first interpretation of the results and written the first draft of the paper. All authors have contributed to the editing and finalization of the paper.

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