Reconstructing volcanic plume evolution integrating satellite and ground-based data : Application to the 23 rd November 2013 Etna eruption

Recent explosive eruptions recorded from different volcanoes worldwide (e.g. Hekla in 2000, Eyjafjallajökull in 15 2010, Cordón-Caulle in 2011) demonstrated the necessity of a better assessment of the eruption source parameters (e.g. column height, mass eruption rate and especially the Total Grain-Size Distribution – TGSD) to reduce the uncertainties associated with the far-travelling airborne ash mass. To do so, volcanological studies started to integrate observations in order to use more realistic numerical inputs, crucial for taking robust volcanic risk mitigation actions. On 23 November 2013, Etna volcano (Italy) erupted producing a 10-km height plume, from which two volcanic clouds were observed at two different altitudes from 20 satellite (MSG-SEVIRI, MODIS). One was described as mainly composed by very fine ash (i.e. PM20), whereas the second one as made of ice/SO2 droplets (i.e. not measurable in terms of ash mass). Atypical north-easterly winds transported the tephra from Etna towards the Puglia region (southern Italy), permitting tephra sampling in proximal (i.e. ~5-25 km from source) and medial areas (i.e. Calabria region, ~160km). Based on the field data analysis, we estimated the TGSD but the paucity of data (especially related to the fine ash fraction) prevented it from being entirely representative of the initial magma fragmentation. 25 To better estimate the TGSD covering the entire grain-size spectrum, we integrated the available field data with X-band weather radar and satellite retrievals. The resulting TGSD is used as input for the FALL3D tephra dispersal numerical model to reconstruct the tephra loading and the far-travelling airborne ash mass. The optimal TGSD is selected by solving an inverse problem through a best-fit with the field, ground-based and satellite-based measurements. The results suggest a total erupted mass of 1.2 × 10 kg, which is very similar to the field-derived value of 1.3 × 10 kg, and a TGSD with a PM20 fraction between 30

worldwide (e.g.Hekla in 2000, Eyjafjallajökull in 2010, Cordón-Caulle in 2011) demonstrated the necessity for a better assessment of the eruption source parameters (ESPs; e.g.column height, mass eruption rate, eruption duration, and total grain-size distribution -TGSD) to reduce the uncertainties associated with the far-travelling airborne ash mass.Volcanological studies started to integrate observations to use more realistic numerical inputs, crucial for taking robust volcanic risk mitigation actions.On 23 November 2013, Etna (Italy) erupted, producing a 10 km height plume, from which two volcanic clouds were observed at different altitudes from satellites (SEVIRI, MODIS).One was retrieved as mainly composed of very fine ash (i.e.PM 20 ), and the second one as made of ice/SO 2 droplets (i.e.not measurable in terms of ash mass).An atypical north-easterly wind direction transported the tephra from Etna towards the Calabria and Apulia regions (southern Italy), permitting tephra sampling in proximal (i.e.∼ 5-25 km from the source) and medial areas (i.e. the Calabria region, ∼ 160 km).A primary TGSD was derived from the field measurement analysis, but the paucity of data (especially related to the fine ash fraction) prevented it from being entirely representative of the initial magma fragmentation.To better constrain the TGSD assessment, we also estimated the distribution from the X-band weather radar data.We integrated the field and radar-derived TGSDs by inverting the relative weighting averages to best fit the tephra loading measurements.The resulting TGSD is used as input for the FALL3D tephra dispersal model to reconstruct the whole tephra loading.Furthermore, we empirically modified the integrated TGSD by enriching the PM 20 classes until the numerical results were able to reproduce the airborne ash mass retrieved from satellite data.The resulting TGSD is inverted by best-fitting the field, ground-based, and satellite-based measurements.The results indicate a total erupted mass of 1.2 × 10 9 kg, being similar to the fieldderived value of 1.3 × 10 9 kg, and an initial PM 20 fraction between 3.6 and 9.0 wt %, constituting the tail of the TGSD.

Introduction
Volcanic explosive eruptions pose hazards related to the release of large quantities of material into the atmosphere.The observation of the eruption features, such as the eruptive column, the tephra loading, or the far-travelling volcanic plume, aims at characterizing the eruption source parameters (ESPs).Hazard assessment related to tephra dispersal, and its implications for aviation safety and public health, is one of the Table 1.Field measurements (locations, loadings, and modes) with the computed tephra loadings obtained with the ARPAE database for the explored TGSDs (Fig. 5).
To mitigate the risk to aviation traffic, nine VAACs (volcanic ash advisory centres) were created worldwide for volcanic cloud monitoring purposes.By making use of operational volcanic ash transport and dispersion models, VAACs aim at alerting for the presence of volcanic ash in the atmosphere.
Besides other ESPs (e.g.eruption start and duration, column height, and mass eruption rate -MER), such models require the total grain-size distribution (TGSD) as input (e.g.Folch, 2012), being one of the most critical ESPs, significantly affecting tephra dispersal model outputs (e.g.Scollo et al., 2008;Beckett et al., 2015).Typically, the TGSD is derived from the field sample analysis through the Voronoi tessellation method (Bonadonna and Houghton, 2005).However, collecting field data on tephra deposit greatly depends on the atmospheric conditions, land/sea deposition, site accessibility, etc.As a consequence, for an inadequate sample dataset in terms of sampling distance from the source (Andronico et al., 2014;Costa et al., 2016a), spatial distribution, and density of samples (Bonadonna et al., 2015;Spanu et al., 2016), the field-derived TGSD is uncertain and cannot be assumed as representative of the whole tephra loading and dispersal.Additionally, the atmospheric residence time of the very fine ash (i.e.hereinafter in this work PM 20 ), ranging from hours to weeks (Rose and Durant, 2009), prevents any rapid deposition, implying their substantial under-estimation within the TGSD (Bonadonna et al., 2011).This raises the necessity to integrate field data with measurements from other sensors (e.g.ground-based radar and satellite) capable of retrieving the missing information in terms of airborne ash.Moreover, the recent eruptions (e.g.Hekla in February 2000, Eyjafjallajökull in April 2010, Cordón-Caulle in June 2011) have shown the impact of the very fine ash on air traffic (e.g.Guffanti et al., 2010;Folch et al., 2012;Sulpizio et al., 2012), but also on public health (e.g.respiratory diseases; Andronico and Del Carlo, 2016;Tomašek et al., 2016;Horwell et al., 2017).
The non-existence of a single instrument capable of covering entirely the grain-size spectrum motivated this study in proposing a method based on the synergic use of field, ground-based, and satellite data for better constraining the TGSD, and therefore the numerical simulations (here FALL3D; Costa et al., 2006;Folch et al., 2009) to reconstruct the tephra loading and the far-travelling airborne ash dispersal.Actually, excluding a few studies (Bonadonna et al., 2011;Folch et al., 2012), simulations are commonly run by using the field-based TGSD or adopting subjective parameterizations (e.g.assuming a constant mass fraction for fine ash).Here, we expanded the reconstruction of the tail of the field-derived TGSD by using radar and satellite retrievals.
We applied this methodology to the 23 November 2013 Etna paroxysm, which occurred from the New South-East Crater (hereinafter NSEC), being the most active crater in the last 20 years (Behncke et al., 2014;De Beni et al., 2015).Atypical winds dispersed the plume north-eastwards, driving the tephra towards the Calabria and Apulia regions (∼ 400 km from the source), where ash fallout was reported (Bonaccorso et al., 2014;Andronico et al., 2015;Montopoli, 2016).A few hours after the eruption, tephra was sampled along the plume axis from Etna (i.e.5-25 km from the NSEC) to Calabria (i.e.∼ 160 km; Fig. 1 and Table 1).Meanwhile, the eruption benefited from being observed through ground-based (i.e.X-band weather radar -X-Radar and L-band Doppler radar -VOLDORAD 2B) and satellitebased (i.e.infrared satellite radiometer) remote sensing instruments.Although they operate in different parts of the electromagnetic spectrum, their integration aims at providing a more complete view of the eruption, especially of the plume dynamic.
In Sect.2, the paper presents the 23 November 2013 Etna eruption, the field, and remote sensing data.Section 3 reports the TGSD estimation, the modelling approach, and the methodology used to reproduce the eruption features.Section 4 is devoted to the results together with their discussions.Section 5 presents the main concluding remarks.1.

The 23 November 2013 lava fountain
In 2013, the 17th lava fountain episode took place on 23 November from the NSEC (De Beni et al., 2015).Mild Strombolian explosions initiated on 22 November afternoon and increased after 07:00 of the following day (all times are in UTC).The transition between Strombolian and lava fountaining activity (i.e. between resumption and the paroxysmal phase; Alparone et al., 2003) started at 09:30, producing intense lava fountains which increased rapidly in height and intensity.During the 50 min of duration of the paroxysmal phase, a sustained 10 km height eruptive column was observed (Bonaccorso et al., 2014;Andronico et al., 2015).Moreover, a peculiar feature was recorded from INGV-OE, showing a greyish volcanic plume that rose above a denser brownish one, from which tephra fallout was visible (Fig. 2).Such observation is attributed to the release of a large amount of water vapour/gas rising higher than tephra (Corradini et al., 2016).This is relevant for characterizing the far-travelling airborne ash, which becomes more complex with the presence of two distinct volcanic clouds.In this case, volcanic ash in the far-field region was testified by a A319 pilot flying over the Albanian coasts at 13:50 and 10.3 km a.s.l.(above sea level), i.e.FL 339, reporting ash between 10.9 and 11.5 km a.s.l., i.e.FL 360-380 (Crompton and Husson, 2015).

Field data
Samples were collected and tephra loading per unit area measured at seven locations (Fig. 1 and Table 1).They were oven-dried at 110 • C for 12 h and analysed in the sedimentology laboratory at INGV-OE, in Catania (Italy).The individual grain-size distributions (GSDs; available in the Supplement in Fig. S1) were measured optically at a 1 -interval through the CAMSIZER ® (Retsch Technology), covering the range from −5 to 5 (where d = 2 − , with the diameter d in millimetres).Although field measurements are commonly used for determining the total erupted mass (TEM) by integrating the isomass lines (Bonadonna andCosta, 2012, 2013), the paucity of samples with their wide dispersion (Fig. 1) limits the reliability of the estimation based on field observations only.However, on the basis of the field data analysis, Andronico et al. (2015) estimated a TEM of 1.3 ± 1.1 × 10 9 kg making use of the Weibull distribution method (Bonadonna andCosta, 2012, 2013).Then, combining the field-derived TEM with the paroxysmal duration (∼ 50 min), they calculated an average MER of 4.5 ± 3.6 × 10 5 kg s −1 .Furthermore, considering the climax phase only (i.e. from 09:55 to 10:14), the MER reached 10 6 kg s −1 , ejecting more than 80 wt % of the erupted mass (Donnadieu et al., 2017).It is worth noting that such MER estimations represent average (or peak) values for the entire duration of the paroxysmal phase without considering its time evolution (i.e. the variation of eruption intensity).Indeed, the time-series MER can be assessed from the relationships between MER and the column height (e.g.Mastin et al., 2009;Degruyter and Bonadonna, 2012;Woodhouse et al., 2013;Folch et al., 2016) and from velocity variations at the vent recorded by VOLDORAD 2B.

Satellite and ground-based remote sensing data
The simultaneous record of the eruption from both satellites and ground-based instruments permits retrieval, on the one hand, of the plume spreading and airborne ash mass dispersal (see Animation A1 in the Supplement), collected by the Spinning Enhanced Visible and Infrared Imager (SE-VIRI) on board the geostationary Meteosat Second Gen- eration (MSG) satellite.The Moderate Resolution Imaging Spectro-radiometer (MODIS) aboard the NASA-Aqua polarorbit satellite was also used to describe the eruption features (Corradini et al., 2016).On the other hand, concerning ground-based instruments, the X-Radar (Montopoli, 2016;Vulpiani et al., 2016) and the visible/thermal cameras (Corradini et al., 2016) provided time-series data of the plume height and the erupted mass.
The available data mentioned above were integrated through a multi-disciplinary approach in Corradini et al. (2016) to improve the volcanic cloud retrievals and the source characterization, and to generate new products.In particular, the satellite observations (Fig. 3) showed the formation of the two distinct volcanic clouds described in Sect. 2. Although both spread north-eastwards, one reached ∼ 6 km a.s.l., being mainly made of ash (ash cloud -AC), and therefore retrieved in terms of airborne ash mass and cloud altitude.The second cloud was higher (∼ 11 km a.s.l.) with enough ice/gas droplets (ice/gas cloud -IC) to significantly alter the cloud characteristics, blinding the satellite from any ash mass measurement (Prata and Kerkmann, 2007).Initially, the clouds were united and split out over the Calabria region (around 11:00).In a final stage, the AC reached the Apulia region, whereas the IC moved over the Ionian Sea towards Albania (around 14:00).In terms of mass, Fig. 4 shows ash was dominant from the onset of the eruption until 11:30, and then ice replaced ash.In fact, from SEVIRI retrievals, ash was likely released between 10:00 and 12:00 prior the emitted water vapour being transformed into ice (i.e.11:00-12:45).This is also shown in Fig. 4, where ice formation starts later than SO 2 and ash emission.SO 2 was released all along the eruption (i.e.10:00-12:30), although with a lower contribution than ash and ice.
The data integration presented in Corradini et al. ( 2016) permits us to reduce the uncertainties associated with the volcanic cloud top height, the ash/ice/SO 2 masses (Fig. 4), and the aerosol optical depth (AOD) retrievals.On the basis of the satellite and X-Radar data, Corradini et al. ( 2016) improved the mass estimation of 30 % and reported a X-Radar-derived TEM of ∼ 3.0 × 10 9 kg with a PM 20 fraction between 1 and 2 wt %, that is ∼ 30-60 t.The source characterization can also be better described by means of the ESP and the eruptive phases.The plume height time series was recorded from the visible cameras at INGV-OE, indicating values from the NSEC (∼ 3300 m a.s.l.) to ∼ 11 km a.s.l., with a rapid increase around 09:30 followed by a decay at 10:20.
The VOLDORAD 2B radar is a pulsed Doppler radar operating at 23.5 cm wavelength (L-band) allowing tephra from block-sized to lapilli-sized to be detected.VOLDO-RAD 2B has continuously monitored Etna's summit craters since 2009 (Donnadieu et al., 2015(Donnadieu et al., , 2016) ) at 3 km from the NSEC (La Montagnola station).Inferred radar parameters (e.g.backscattered echo power) are proportional to the quantity of tephra detected through the radar beam.In addition, the along-beam radial velocities permit lava fountains from being observed at high time resolution (i.e.0.2 s), inferring near-source detection of the ejection velocities by means of the following equation (Freret-Lorgeril et al., 2016;Donnadieu et al., 2017): where V e is the ejection velocities (in m s −1 ), v r+ is the radial velocity (in m s −1 ), and θ is the elevation angle of the radar beam (here θ = 14.9 • ).Such an approach is relevant for integrating the time-dependent ejection velocities with the corresponding observed eruptive column heights.In particular, we used the VOLDORAD 2B data associated with the 23 November 2013 eruption to better constrain the eruption phases' characterization.

Methodology
Simulating the tephra loading and airborne ash dispersal of the 23 November 2013 Etna eruption requires us to assess the related ESPs, and in particular the TGSD.Their use as input parameters into the FPlume model (Folch et al., 2016) aims at describing the eruption column, representing the source term required by the FALL3D tephra dispersal model (Costa et al., 2016b).In the following methodology, we present the TGSD reconstruction and modelling approach.Then, the simulations are analysed in terms of tephra loadings and airborne ash mass dispersal to best fit the field and satellite measurements.

TGSD estimation
The seven field samples are not sufficient for assuming the field-derived TGSD as the full spectrum TGSD (Andronico et al., 2014;Beckett et al., 2015;Bonadonna et al., 2015;Costa et al., 2016a;Spanu et al., 2016).Although such a field-based TGSD is being biased toward coarse ash, we first estimated the TGSD (hereinafter Field TGSD; Fig. 5) from the individual GSDs using the Voronoi tessellation method (Bonadonna and Houghton, 2005).However, the Field TGSD needs to be better characterized prior to being used within atmospheric ash dispersal models.Considering the Field TGSD representativeness on the grain-size spectrum (i.e.−5 to 5 ; Sect.2.1), we used the X-Radar retrievals to constrain the mass relative to coarse and fine ash (i.e.−1 to 5 ; Corradini et al., 2016).The X-Radar-derived TGSD is inverted from the particle-size distribution (PSD), given as ash number density distribution (Corradini et al., 2016).It is worth noting that we considered a spatial and temporal average of the X-Radar-based PSD for the whole event.The average takes in input each PSD estimated from each single radar resolution volume delineated by horizontal angle, vertical angle, and range distance at each available time step for the airborne ash mass seen by the radar.We converted the PSD into number of particles per unit of volume with the particle-size bins.Then, by means of the volume and density associated with the size bins, we calculated the mass density distribution (hereinafter Radar TGSD; Fig. 5).However, we would like to highlight that retrieval of Radar data is done assuming a Gamma distribution for the number particles per unit of volume for each particle size interval.Then this distribution is converted to express the mass fraction as a function of .In particular, since a single Gamma distribution is not able to adequately describe large size spectra, a Gamma distribution, with different parameters, is assumed in each particle size range of fine ash, coarse ash, small lapilli, and large lapilli, so the final total distribution is a combination of several gamma distributions.However, such an empirical derived distribution can be approximated using other distributions, such as a lognormal or a Weibull distribution.The latter point will be investigated in future studies.
It is worth noting that the Field and Radar TGSDs are distributions observed through their own grain-size window, which explains the substantial difference in shape (Fig. 5).It follows that assessing accurately the TGSD covering both windows can be done by integrating the Field and Radar TGSDs only.Although, in principle, their integration is possible, the grain-size windows' discrepancy prevents merging of the Field and Radar TGSDs without knowing their relative weighting averages.We determined empirically the weight combination by integrating the distributions at regular intervals (i.e. from full Field TGSD to full Radar TGSD).The resulting distribution (i.e.−5 to 5 ; hereinafter Integrated TGSD; Fig. 5) is obtained, best fitting the tephra loading at the sampled sites.
However, due to the instrument/method grain-size limit, none of the three TGSDs (Field, Radar, or Integrated TGSD; Fig. 5) contains enough PM 20 to reproduce the far-travelling airborne ash mass retrieved by satellite.We assessed the tail of the Integrated TGSD (i.e.≥ 6) by modifying empirically the PM 20 fraction, adding mass into the corresponding classes.We calculated the fractions based on an empirical power-law dependence of the classes with through the following parameterization: where X( i ) is the fraction (in wt %) allocated to the ith bin, X( 5) is the fraction obtained for = 5, and γ is the empirical factor (γ < 1).The explored γ values span from 0.1 to 0.7, giving, respectively, PM 20 fractions between ∼ 0.6 and 10.7 wt % of the TEM.The best fraction to use within the TGSD (hereinafter Whole TGSD; Fig. 5) is chosen as best fitting the satellite retrievals.

Modelling approach
To furnish the ESPs required by the FALL3D tephra dispersal model, we used the FPlume integral plume model (Folch et al., 2016) describing the eruptive column based on the buoyant plume theory (Morton et al., 1956).FPlume solves a set of 1-D cross-section-averaged equations for mass, momentum, and energy conservation in the eruption column, accounting for wind coupling, air moisture, particle re-entrainment, and ash aggregation effects (Folch et al., 2016).Among the source conditions, FPlume feeds into FALL3D by describing the mass flow rate for each particle bin and the vertical distribution within the column.As inputs, FPlume uses the TGSD, initial magma temperature, and water content (Table 2) to calculate the mass released per unit of time within the column.Indeed, FPlume uses the TGSD to solve the mass conservation equation for each class distributing along the column.Then, the mass for each particle class at each level is transported laterally using FALL3D.
In our case, Etna's magmas have a temperature of 1300 K with ∼ 2.5 wt % of water (Carbone et al., 2015;Spilliaert et al., 2006).FPlume calculates MER from the column height (or vice versa) for a given wind profile (Folch et al., 2016) by describing the air mixing within the plume through two turbulent air entrainment coefficients (i.e.radial -α and crossflow -β; Bursik, 2001;Kaminski et al., 2005;Suzuki and Koyaguchi, 2015;Folch et al., 2016;Costa et al., 2016b).Here, α and β are obtained empirically through the solution of an inverse problem best-fitting the erupted mass derived from the field measurements (Poret et al., 2017).Ash aggregation can be considered negligible during Etna eruptions, with less than 2 wt % of the fine ash removed by aggregation.For this reason, we did not consider such a process in this study.The effects of the typical uncertainties associated with the input parameters of FPlume on the source term characterization are described in Macedonio et al. (2016).
FALL3D is used for simulating tephra dispersal and is a 3-D time-dependent Eulerian model based on the advectiondiffusion-sedimentation equation computed over a terrainfollowing domain (Costa et al., 2006;Folch et al., 2009).Besides the ESPs, FALL3D needs the time-dependent meteorological fields over the computational domain for the corresponding period (i.e. from 00:00 on 23 November up to 00:00 on 29 November 2013).The first series of simulations are run by means of a local high-resolution meteorological database (ARPAE from INGV-OE) to better constrain the computed tephra loadings against the field measurements in proximal and medial areas (Fig. 1 and Table 1).Indeed, ARPAE provides a 7 km × 7 km spatial and 15 min temporal resolution over the domain highlighted in Fig. 1.Then, FALL3D internally interpolates the meteorological data over a grid set at 1 km × 1 km resolution.The parameterizations used for the simulations with the ARPAE database are summarized in the Appendix.The related main atmospheric profiles (e.g.temperature, air moisture, and wind speed) over the NSEC are displayed in Fig. 6.
The second series of simulations aims at reproducing the satellite retrievals, expanding the computational domain to Albania.The ARPAE data do not cover such a domain, for which we use the meteorological fields from the European Center for Medium-Range Weather Forecasts (ECMWF, ERA-Interim-Reanalysis; hereinafter ERA-Interim).They provide a 6 h interval for 37 pressure levels of data at 0.75 • horizontal resolution.For computational cost reasons, the internal grid resolution into FALL3D is set at 5 km × 5 km, which is still consistent with the satellite data resolution (3 km × 3 km at nadir).The parameterization used with the ERA-Interim database is summarized in the Appendix.
The consistency between the two databases is checked by adding the profiles retrieved over the NSEC with ERA-Interim in Fig. 6.Although ARPAE and ERA-Interim tend to have the same temperature and wind speed patterns, the air moisture from ERA-Interim is slightly lower than ARPAE for 3-6 km a.s.l. and higher for 7-11 km a.s.l.These observations are not significant to produce a substantial effect on the simulations.Moreover, Fig. 6 also shows the conditions over the Albanian capital (Tirana).With such meteorological conditions, airborne tephra needs 4.5 h to be transported from Etna to Albania (Fig. 6), being consistent with the pilot report mentioning ash.Wind speed is moderate to strong, with higher velocities near the volcano than at Tirana.As indicative values at 9 km a.s.l., we report ∼ 48 and ∼ 45 m s −1 over the NSEC (at 09:30) for ERA-Interim and ARPAE, respectively, and ∼ 34 m s −1 over Tirana at 14:00.Besides the velocities, the wind direction (Fig. 6) shows a strong northeasterly orientation over the NSEC, which is consistent with the tephra dispersion towards Calabria.The profiles indicate a substantial variation between middle (5-6 km a.s.l.) and high altitudes (> 7 km a.s.l.), which probably resulted in the different spreading orientations for the two volcanic clouds (AC and IC) at their own altitudes (Fig. 3).Besides the profiles, the consistency for using alternatively the two meteorological databases is checked by constraining the simulations with ERA-Interim to converge the TEM towards the same value as for the Integrated TGSD and the ARPAE database.
Tephra dispersal simulations are commonly carried out using the field-based TGSD and assuming a constant average column height (or MER) for the entire duration of the paroxysmal phase (Fig. 7a).However, it is evident that eruption intensity varies substantially with time and consequently the column height (e.g.Scollo et al., 2014Scollo et al., , 2015)).To account for such variability, we discretized the eruption into a set of phases consistent with (i) the plume height observations from the remote sensing measurements (Corradini et al., 2016) and (ii) the exit velocities retrieved by VOLDORAD-2B (Donnadieu et al., 2015(Donnadieu et al., , 2016(Donnadieu et al., , 2017)).The improved simulation scheme (Fig. 7b and c) is achieved by coupling this discretization with the ARPAE or ERA-Interim databases and the Integrated TGSD or Whole TGSD, respectively, depending on the inversion purpose.

Inversion modelling strategy
Simulation optimization is carried out to assess the ESP, and among them the TGSD, leading to the numerical reconstruction of the tephra loading and airborne ash mass dispersal.Input parameters in Table 2 were varied at constant steps within their ranges facing the inherent non-uniqueness solution for assessment purposes (e.g.Anderson and Segall, 2013).Starting by inverting the Integrated TGSD, we tested each weighting average combination of the Field and Radar TGSDs, ranging from 100 wt % Field TGSD to 100 wt % Radar TGSD, with a step of 5 wt %.To select the best combination, we compared the tephra loadings computed at the sampled sites until we best fit the field measurements.
Considering the simulations, we used the scheme described in Sect.3.2 (Fig. 7b and c), which implies a set of column height values (and hence the corresponding MERs) with the average exit velocity.Therefore, neither the column height, the MER, nor the exit velocity were changed in each simulation.However, we inverted the plume parameters (i.e.α and β) from 0.05 to 0.15 and 0.05 to 1.0, respectively (Costa et al., 2016b), by means of the following goodnessof-fit procedure.
The goodness of fit between simulations and field observations was evaluated through different statistical metrics (see Poret et al., 2017).In particular, we used the normalized root mean square error (i.e.RMSE) assuming three different error distributions (i.e.RMSE 1 , RMSE 2 , and RMSE 3 ) described in Folch et al. (2010).We also used the Aida (1978) indexes K (i.e.geometric average of the distribution) and k (i.e.geometric standard deviation of the distribution).
where i refers to the ith sample over N , and Sim and Obs are the simulated and observed tephra loadings, respectively.For a given set of ESPs, K gives the gap between the theoretical optimal tephra loading samples and the simulated ones.The reliability of the simulation is obtained for K between 0.95 and 1.05, which means a threshold of ±5 wt % from the derived theoretical optimal TEM.It follows that the best simulations are selected for K close to 1 with k and the three RM-SEs minimized.Additionally, we estimated the bias, the correlation, and the Student T test (Folch et al., 2010).
After the Integrated TGSD, the Whole TGSD is inverted by quantitatively analysing the effect of different PM 20 fractions (i.e.0.6-10.7 wt %; Sect.3.1) on the computed airborne ash dispersal.The best fraction is selected by means of the following three statistical metrics.The mass difference (i.e.Mass) between the satellite measurements and the FALL3D estimates.We compared the masses over the number of pixels given by the plume mask (obtained for the threshold of 0.1 t km −2 ) retrieved from SEVIRI: where M Obs and M Sim are the observed and simulated masses integrated over the whole event (i.e. from t 0 = 09:30 to t f = 14:30, with T = t f − t 0 ).This index gives the discrepancy (in t) for each γ factor (i.e.PM 20 fractions).Additionally, we also calculated for each γ factor the absolute average difference of mass per unit area (Sum( ) in t km −2 ) for the entire volcanic cloud by the following: where N is the number of pixels (i.e.plume mask), and M Obs (N ) and M Sim (N ) are the observed and modelled masses associated with the Nth pixel for SEVIRI and FALL3D, respectively.Area P refers to the area covered for the related time interval, which is calculated by means of N and the pixel resolution (i.e. 9 km 2 ).This index indicates the uncertainty of the simulated airborne ash mass per unit area with respect to the satellite retrieval.
Considering that Mass and Sum( ) are discrepancy estimates, the selection is done on the basis of their minimization.Nonetheless, Sum( ) gives absolute values preventing any over-or under-estimation characterization.It follows that we also evaluated the following index: where ε refers to an over-estimation per pixel when ε < 0 and an under-estimation per pixel for ε > 0, with a best fit for ε = 0.Moreover, the index indicates the average mass difference per unit area (i.e.t km −2 ) between the satellite measurements and the simulation.The synergic use of these metrics aims at providing a simple way of comparing spatially and temporally the simulation outputs with the field and remote system measurements.

Results and discussion
This section describes the results of the inversion of (i) the ESPs, and among them, (ii) the Integrated TGSD reproducing the tephra loading.Then, (iii) we report the results for assessing the PM 20 fraction needed within the Whole TGSD to capture the airborne ash transported in distal area.

ESP inversion
Regarding the Integrated TGSD inversion (Sect.3.3), Table 3 shows the statistical analysis for the best simulation (i.e.75,25), being selected as the best weighting average combination for composing the Integrated TGSD (Table 3 and Fig. 8).It is worth noting that RMSE 2 and k indicate relatively high values yielding a mean error factor nearby 3, which is comparable to uncertainties associated with other classical methods (Bonadonna andCosta, 2012, 2013;Bonadonna et al., 2015).Figure 9 illustrates the statistical analysis of the Whole TGSD inversion (Sect.3.3) for the best simulation for each PM 20 fraction.Considering the whole airborne ash mass, the results yield a best value for Mass at γ = 0.65 (i.e.PM 20 = 9.0 wt %), indicating an overall underestimation of ∼ 76 t of ash by FALL3D for the entire eruption.Then, Sum( ) shows a minimum for γ = 0.40 (i.e.PM 20 = 3.6 wt %), giving an absolute average difference of mass per unit area of ∼ 0.37 t km −2 for the whole sequence.The third index returns a best value of ε = −0.03t km −2 for γ = 0.65 (i.e.PM 20 = 9.0 wt %), being consistent with Mass.ε likely reflects that FALL3D slightly over-estimates the average mass per pixel of 0.03 t km −2 .By integrating the results (Fig. 9), the Whole TGSD required the minimum PM 20 fraction of 3.6 wt % to best reproduce in absolute terms the average ash mass per unit area.However, such a fraction is not sufficient for best simulating the whole airborne ash mass released during the eruption, and minimizing the over-or underestimation, which tends to be satisfied with higher PM 20 fractions (i.e.9.0 wt %).The corresponding input TGSD is www.atmos-chem-phys.net/18/4695/2018/Atmos.Chem.Phys., 18, 4695-4714, 2018   2).
displayed in Fig. 5.Moreover, Mass and ε in Fig. 9 both indicate that FALL3D under-estimates substantially the airborne mass for PM 20 fractions lower than ∼ 7 wt % and overestimates it above ∼ 10 wt %.
Regarding the other ESPs, although the column height values were not changed throughout the simulations (Fig. 7b  and c), we report here the MER inverted by FPlume for the climax phase only, which is ∼ 7.0 × 10 5 kg s −1 .The calibration of α and β returns values ranging from 0.06 to 0.15 and from 0.21 to 1.00, respectively, depending on the weighting average combination (Table 3).The latter ranges are con-sistent with the literature (Devenish et al., 2010;Suzuki and Koyaguchi, 2015).

Tephra loading
During the Integrated TGSD inversion, the six proximal samples were relatively stable when varying the weighting average combination, whereas the farthest sample (i.e.TER) was substantially affected.Figure 8 shows the comparison between the computed and measured tephra loadings with the Integrated TGSD (details in Table 1).It is worth noting that making use of the Field TGSD prevents FALL3D from cap- turing the TER sample, while the Radar TGSD fails on most of the samples as indicated in Table 1.These observations argue the necessity of combining the two different distributions through the Integrated TGSD, especially when field measurements are few.Figure 8 shows the seven samples lying within the 1/5-5-times threshold of the measured tephra loadings, especially the unique medial sample (i.e.TER).As indicative values from Table 1, the six proximal samples indicate tephra loadings ranging from 1 to 17 kg m −2 .In contrast, FALL3D computed them between 3 and 7 kg m −2 for the Integrated TGSD.Such narrower ranges compared to the field data can be attributed to the complexity for modelling in proximal areas (< 20 km from the source) and the field samples' location with respect to the main plume axis.
Besides the tephra loadings, we also compared the fieldderived GSD at the sampled sites with the numerical results for the Integrated TGSD (see Fig. S1).Although FALL3D reproduces accurately three of the seven samples by peaking at the same modes, four proximal samples (i.e.CRT, PDM, FFD, and GDN) are shifted by 1 , indicating the field measurements being slightly finer than the computed ones.This discrepancy argues for the difficulty in computing accurately at such proximal areas due to plume dynamic complexities (e.g.Cerminara et al., 2016).Nonetheless, the mode shift can also be attributed to the sampling distance from the source as explained in Spanu et al. (2016).Indeed, in proximal areas the coarse tephra (−4 ≥ ≥ −2) deposits rapidly, increasing the difficulty in estimating accurately this part of the TGSD with the Voronoi tessellation method together with a paucity of field measurements (Andronico et al., 2014).Moreover, we cannot exclude partial breakages of a few coarse-grained clasts when impacting the ground (Andronico et al., 2015), which also may result in grain sizes slightly finer than expected.
Although we used the improved simulation scheme (Sect.3.2; Fig. 7b), we run a simulation through the simplified procedure (Fig. 7a) to highlight the effect on the tephra loading, and therefore the statistical analysis.The results show that making use of a constant plume height (here ∼ 11.3 km a.s.l.) for the entire paroxysmal phase gives K = 1.01 and k = 5.76 with RMSE 1 = 0.80, RMSE 2 = 3.36, and RMSE 3 = 1.33, which are significantly higher than for the improved procedure (details in Table 3).Regarding the TEM, the simplified scheme returns 1.5 × 10 9 kg, which is ∼ 34 % higher than for the integrated approach with 1.2 × 10 9 kg.The latter TEM is in good agreement with the estimation of 1.3 × 10 9 kg reported in Andronico et al. (2015).It is worth noting that varying the weighting av- erage from 100 wt % Field TGSD towards 100 wt % Radar TGSD yields an increasing TEM going from 1 to 6 × 10 9 kg, respectively (Table 3).This observation of TEM is consistent with the results described in Corradini et al. (2016), which indicates a X-Radar-derived total mass of 3.0 × 10 9 kg compared to the field-derived TEM of 1.3 × 10 9 kg from Andronico et al. (2015).Such a difference between X-Radar and field-based TEM estimates can be explained by considering the following aspects: (i) X-Radar samples airborne particles during their fallout, whereas the field measurements are based on deposited tephra; (ii) the operative window focuses the X-Radar retrievals on detecting the ash particles (−1 to 5 ), while the field sampling method expands the measurements to block-sized (−5 to 5 ); (iii) the Radar TGSD refers to the average over the duration observed from the radar at the sampled grid points, which does not necessarily coincide with the duration and location characterized by the Field TGSD; (iv) as explained in Sect.3.1, the X-Radar measurements are made with assumptions using a regression model of radar simulations, which can add a further degree of uncertainty.The assumptions mainly affecting the final radar retrieval involve the radar forward model used to set up the radar retrieval scheme.It follows that assumptions made about particle shape, density, orientation, and PSD play the key role.However, the presented integrated approach by weighting the distributions issued from different methods aims at preventing the resulting Integrated TGSD from being associated with the full uncertainty of a single source.
The use of the different distributions (i.e.Field, Radar, Integrated, and Whole TGSDs) presented in this study permits comparison of the resulting tephra loading maps (Fig. 10).The tephra loading scale reported in Fig. 10 refers to the use of the ERA-Interim database, indicating slightly different tephra loadings than the values in Table 1 (ARPAE).Here, Fig. 10 is used as indicative tephra loading maps to display the effect of the input TGSD on the resulting tephra dispersal, showing the affected areas (e.g. the Calabria and Apulia regions).In particular, the use of the Field TGSD (Fig. 10a) permits FALL3D to compute the tephra loadings at the sampled sites up to Calabria, but not in the Apulia region where ash was reported.The Radar TGSD (Fig. 10b) operates in the ash window, preventing its use from reproducing any tephra loading and airborne ash data.In contrast, the Integrated and Whole TGSDs (Fig. 10c and d) capture all the tephra loading samples, but only the Whole TGSD succeeds in simulating the far-travelling airborne ash mass retrieved from the satellite.The corresponding time-series animation of the tephra loading associated with the Whole TGSD is available as the Supplement (Animation A2).
Figure 11.Illustration of the comparative study between the SEVIRI and FALL3D airborne ash masses for a given time (i.e.12:00, 13:00 and 14:00) to best reproduce the satellite retrievals.

Airborne ash dispersal
As mentioned in Sect.2, large quantities of ash, water vapour (transformed into ice), and SO 2 gas (Fig. 4) were released from Etna, preventing the remote systems from quantifying the whole event easily.The formation of two volcanic clouds (AC and IC) following their own trajectory at different altitudes (Fig. 3) increased substantially the complexity of comparing quantitatively the far-travelling airborne ash masses (i.e.SEVIRI and FALL3D).Indeed, the columnar satellite measurements and FALL3D results prevent the two clouds from being isolated, which motivated this study to focus on the plume mask retrieved by SEVIRI for each time (Fig. 11).Figure 11 illustrates the comparison between the retrieved and computed airborne ash mass.By means of the inverted PM 20 range (i.e.3.6-9.0wt %), we displayed the airborne ash mass maps.The left column refers to the minimum PM 20 fraction (i.e.3.6 wt %) required to capture accurately the absolute average difference of mass per unit area (i.e.Sum( )), whereas the right column corresponds to the fraction (i.e.9.0 wt %) best reproducing the whole airborne ash mass (i.e.Mass and ε).Each panel in Fig. 11 (i.e. a,b,and c) shows the overlapping between the SE-VIRI retrievals and the FALL3D outputs for a given time.Although the overlap tends to decrease with time (i.e.12:00, 13:00, and 14:00, respectively), the results for γ = 0.65 (i.e.PM 20 = 9.0 wt %) indicate a better performance than for γ = 0.40 (i.e.PM 20 = 3.6 wt %).The entire time-series an-imations are available in the Supplement (Animations A3 and A4 for γ = 0.40 and γ = 0.65, respectively).
The PM 20 range obtained for the 23 November 2013 Etna paroxysm tends to be relatively high with respect to the literature (1-2 wt %; Corradini et al., 2016), eventually attributed to the observational data used and the instrument properties.However, in terms of mass to the TEM, the estimated PM 20 fractions indicate consistent values.Indeed, 1-2 wt % of the X-Radar TEM (3.0 × 10 9 kg) refers to 30-60 t, while 3.6-9.0wt % of the integrated TEM (1.2 × 10 9 kg) gives 43-108 t.In fact, Corradini et al. ( 2016) integrated X-Radar data with satellite retrievals to assess the PM 20 fraction.However, the satellite cannot quantify any ash mass from pixels mainly filled by ice or gas (e.g.SO 2 ).In other words, although the volcanic ice/gas clouds (i.e.IC) are assumed to be produced from ash nucleus (Corradini et al., 2016), the probable presence of ash within such clouds will be missed from SEVIRI.
Being the airborne ash mass spreading downwind towards the far field, the very fine ash fraction (i.e.here 3.6-9.0wt % of the erupted mass) is a critical input into operational tephra dispersal models (e.g.HYSPLIT, Stunder et al., 2007;NAME, Witham et al., 2007;FALL3D, Folch et al., 2012), which are widely used for aviation safety.Although few studies have attempted to better constrain the fraction estimation, eruptions from different volcanoes are not comparable, as such a fraction is very different from one case to the other, ranging from 50 wt % to a few wt % (Rose and Durant, 2009).As discussed by Costa et al. (2016aCosta et al. ( , 2017)), the very fine ash fraction varies with eruption intensity, magma composition, and eruption style.In particular, at the Spurr 1992 eruption, Wen and Rose (1994) estimated ∼ 2 wt % dispersed into the distal area.At the 2010 Eyjafjallajökull eruption, the estimated range spans from ∼ 0.9 to 11 wt % (Bonadonna et al., 2011;Dacre et al., 2011;Devenish et al., 2012).However, some operational models assume a fraction of ∼ 5 wt %, which is not related to our estimate for the Etna eruption.In fact, assuming a constant fraction (e.g. 5 wt %) would represent the very fine ash fraction that escapes to aggregation processes and travels in the far field.In the case of basaltic eruptions, like at Etna, the eruption intensity and the very fine ash content are low, and hence aggregation less efficient (Costa et al., 2010), implying that most of the fraction can be transported distally.These observations yield the necessity for better considering such a fraction as input, suggesting further investigations on both basaltic and silicic volcanoes.
Regarding the FALL3D results in Fig. 11, the airborne ash maps show the two volcanic clouds (AC and IC) observed from the satellite (Corradini et al., 2016), although they are still connected to each other.Dispersing simultaneously from the source, the FALL3D simulations yield the presence of volcanic ash following the trajectory of AC below FL 250.In addition, FALL3D also indicates a major contribution of the airborne mass associated with the IC trajectory spreading over FL 250.The results in terms of temporal dispersal (Animation A3) are corroborated by the SEVIRI retrievals (An-imation A1) and the pilot report, which mentioned volcanic ash and probably gas near Albania at FL 360-380 (Crompton and Husson, 2015).
As a consequence of being blind to any ash within the IC, the comparative study results represent partially the whole airborne ash.This raises questions related to volcanic hazards, such as the air traffic safety.In fact, on the basis of the FALL3D results, the IC appears to have a significant amount of erupted material (i.e.PM 20 , ice, and gas).This observation highlights the necessity for quantifying entirely the far-travelling airborne tephra, perhaps benefitting from other sensors capable of characterizing such aerosol clouds.In particular, this study inferred from quantitative analysis based on the observations in terms of tephra loading and airborne ash mass the interest in integrating retrievals from diverse instruments to assess accurately the initial magma fragmentation (i.e.TGSD of the whole erupted tephra).

Conclusions
Recent studies have shown the need to improve the assessment of the eruption source parameters to reduce the uncertainties and present more realistic numerical outputs, which can be used for hazard mitigation.Here, we worked on better estimating the initial magma fragmentation (i.e. total grainsize distribution -TGSD) by integrating measurements from field samples and ground-based (X-band weather radar) and satellite-based (SEVIRI) systems.We applied the methodology to the 23 November 2013 Etna paroxysm, which benefited from a north-easterly wind direction that dispersed the tephra over Calabria towards the Apulia (Italy) and Albania regions.The available observations in terms of tephra loadings and airborne ash dispersal were used to reconstruct numerically (through the FALL3D model) the eruption features from the source to distal areas.In fact, the field-based TGSD reproduces only the sampled tephra loadings, whereas the Radar TGSD refers to a limited range of ash classes preventing its use within FALL3D as the initial TGSD.We produced an Integrated TGSD (i.e.weighting average of field + radar distributions) to best fit the tephra loadings.The inversion results yield a TGSD made of 75 wt % of the Field TGSD and 25 wt % of the Radar TGSD.However, the Integrated TGSD does not account for the far-travelling airborne ash mass retrieved from satellites (i.e.PM 20 ).We empirically modified the Integrated TGSD to implement the SEVIRI retrievals by investigating diverse PM 20 fractions (i.e.0.6-10.7 wt %) until we best fit the measurements.The inverted PM 20 fraction best matching the SEVIRI data ranges from 3.6 to 9.0 wt %, depending on capturing the whole airborne ash mass or the mass per unit area.This study highlighted the need to improve the integration of data from different instruments to better quantify tephra loading and airborne mass (i.e.PM 20 , ice, and gas), especially when aerosol clouds are produced during the eruption.From a computational point of view, the assessment of the initial TGSD would benefit from such integration, being widely used for modelling purposes such as for air traffic safety.This work aims at being of interest for developing new methods or tools capable of assessing the full size-spectrum TGSD.
Data availability.The ERA-Interim Reanalysis meteorological database was retrieved from the European Center for Medium-Range Weather Forecasts (ECMWF) and the ARPAE database from the INGV-OE archives.The L-band Doppler radar (VOLDO-RAD 2B) data were provided by the open-access database on the OPGC website: http://voldorad.opgc.fr/(Donnadieu et al., 2015).

Information about the Supplement
The supplement associated with this paper serves to illustrate the results in terms of individual grain-size distributions with the Integrated TGSD, which is validated on the basis of the tephra samples (Fig. S1).The time-series animations aim at highlighting the main eruption features (i.e.whole tephra loading and airborne ash dispersal).
Comparison of the 7-individual field-derived GSDs with the computed ones through the FALL3D model.The figure indicates the reproducibility of the local GSD by peaking at the same mode.The shifted GSDs are discussed in Sect.4.2.

Animation A1
The time-series animation refers to the dynamic evolution of the volcanic ash cloud travelling from the source retrieved from SEVIRI (i.e.09:30-14:30 UTC).

Animation A2
The time-series animation corresponds to the simulation of the tephra loading obtained for the Whole TGSD with γ = 0.65.The animation shows the temporal expansion of the tephra fallout indicating the affected areas (i.e.09:30-14:30 UTC).

Animation A3
The time-series animation shows the simulation of the airborne ash dispersal associated with the Whole TGSD produced with γ = 0.40 (i.e.09:30-14:30 UTC).This animation indicates the temporal dispersal obtained with the initial injection of 3.6 wt % of PM 20 into the atmosphere.The major lobe goes towards Albania, which corresponds to the ice/gas volcanic cloud, whereas the minor lobe (i.e.tail) spreads towards the Apulia region (southern Italy) and is related to the volcanic ash cloud.

Animation A4
The time-series animation refers to the simulation of the fartravelling airborne ash dispersal computed with the Whole TGSD for γ = 0.65 (i.e.09:30-14:30 UTC).This animation shows a similar dispersal to Animation A3.However, using γ = 0.65 means the initial injection of 9.0 wt % of PM 20 into the atmosphere, which results in higher ash mass values, especially for the major lobe spreading towards Albania.

Figure 1 .
Figure 1.Tephra sample locations (Sicily and Calabria regions, Italy).(a) shows the local to medial areas (up to ∼ 160 km from the NSEC) affected by the fallout.(b) is a zoom indicating the proximal zone (up to ∼ 25 km from the NSEC) and the dispersion of the samples.Details in Table1.

Figure 2 .
Figure 2. Photograph of the eruption showing the formation of the two volcanic clouds rising at different altitudes (greyish above the brownish).Source: courtesy of Boris Behncke (INGV-OE).

Figure 3 .
Figure 3. Satellite image (SEVIRI) showing the trajectories of the two volcanic clouds (modified from Fig. 17 in Corradini et al., 2016).The ash cloud dispersed towards the Apulia region (southern Italy) at ∼ 6 km a.s.l., whereas the ice/gas cloud moved over Albania at ∼ 11 km a.s.l.

Figure 4 .
Figure 4. Ash, ice, and SO 2 mass time series retrieved from SEVIRI for the 23 November 2013 Etna eruption.

Figure 5 .
Figure 5. Input TGSDs estimated from either field or X-Radar data.The Integrated TGSD emerges from a weighting average combination of the Field and Radar TGSDs.The Whole TGSD derives from the Integrated TGSD modified to implement the satellite measurements.

Figure 6 .
Figure 6.The main meteorological profiles over the NSEC from ARPAE (INGV-OE) and ERA-Interim (ECMWF), and over Tirana for ERA-Interim.

Figure 7 .
Figure 7. Simulation schemes.(a) Simplified procedure.(b) Discretization of the eruption into a set of phases to account for the temporal variation of the intensity (i.e.column height, hence MER, and exit velocity).The improved scheme is accompanied by the Integrated TGSD and ARPAE database.(c) Same procedure as (b) with the Whole TGSD and ERA-Interim.

Figure 8 .
Figure 8.(a) Comparative study between the measured and computed tephra loadings for inverting the Integrated TGSD.(b) Graphic of the k index showing the optimization for assessing the best weighting average combination to apply to the Field and Radar TGSDs (details in Table2).

Figure 9 .
Figure 9. Quantitative analysis of the airborne ash mass measured from SEVIRI and computed by FALL3D to invert the PM 20 fraction to use within the Whole TGSD for best reproducing the SEVIRI retrievals.The upper part compares the whole airborne ash masses for the entire eruption, whereas the middle part gives the difference of the absolute average difference of mass per unit area.The lower part quantifies the difference in terms of mass per unit area (details in Sect.3.3).

Figure 10 .
Figure 10.Tephra loading maps computed with the (a) Field, (b) Radar, (c) Integrated, and (d) Whole TGSDs, respectively.They indicate the relevance of the integrated approach reproducing the affected areas.

Table 2 .
Input parameters used within the FPlume and FALL3D models.Multiple TGSDs are tested as input for the simulations (see Sect. 3.1).The column height, MER, and exit velocity are set as multiple values (see Sect. 3.2).The simulation scheme is presented in Sect.3.2 and Fig. 7.

Table 3 .
(Folch et al., 2016)or the best simulations (i.e.calibration of α and β) for each weighting average combination tested during the inversion of the Integrated TGSD.α is described through α 1 and α 2 within the calculation(Folch et al., 2016).TEM indicates the associated theoretical value for each combination.