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

Matthieu Poret, Stefano Corradini, Luca Merucci, Antonio Costa, Daniele Andronico, Mario Montopoli, Gianfranco Vulpiani, Valentin Freret-Lorgeril

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

## To cite this version:

Matthieu Poret, Stefano Corradini, Luca Merucci, Antonio Costa, Daniele Andronico, et al. Re-constructing volcanic plume evolution integrating satellite and ground-based data: application to the 23 November 2013 Etna eruption. Atmospheric Chemistry and Physics, 2018, 18, pp.4695-4714. 10.5194/acp-18-4695-2018. insu-03708974

HAL Id: insu-03708974

https://insu.hal.science/insu-03708974

Submitted on 30 Jun 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

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

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

Matthieu Poret<fnref n="1" />,<fnref n="2" />, Stefano Corradini<fnref n="3" />, Luca Merucci<fnref n="3" />, Antonio Costa<fnref n="1" />, Daniele Andronico<fnref n="4" />, Mario Montopoli<fnref n="5" />, Gianfranco Vulpiani<fnref n="6" />, and Valentin Freret-Lorgeril<fnref n="7" />

**Correspondence:** Matthieu Poret (matthieu.poret@gmail.com)

Received: 7 December 2017 – Discussion started: 2 January 2018

Revised: 1 March 2018 – Accepted: 14 March 2018 – Published: 6 April 2018

**Abstract.** Recent explosive volcanic eruptions recorded 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₂₀), and the second one as made of ice/SO₂ 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₂₀ 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⁹ kg, being similar to the field-derived value of 1.3 × 10⁹ kg, and an initial PM₂₀ fraction between 3.6 and 9.0 wt %, constituting the tail of the TGSD.

## 1 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

<footnote n="1">¹Istituto Nazionale di Geofisica e Vulcanologia, Bologna, Italy</footnote>

<footnote n="2">²University of Bologna, Geophysics department, Bologna, Italy</footnote>

<footnote n="3">³Istituto Nazionale di Geofisica e Vulcanologia, CNT, Rome, Italy</footnote>

<footnote n="4">⁴Istituto Nazionale di Geofisica e Vulcanologia, Osservatorio Etneo, Catania, Italy</footnote>

<footnote n="5">⁵Institute of Atmospheric Sciences and Climate, National Research Council of Italy, Rome, Italy</footnote>

<footnote n="6">⁶Department of Civil Protection, Presidency of the Councils of Ministers, Rome, Italy</footnote>

<footnote n="7">⁷Université Clermont Auvergne, CNRS, IRD, OPGC, Laboratoire Magmas et Volcans, 63000 Clermont-Ferrand, France</footnote>
<table><caption>**Table 1.** Field measurements (locations, loadings, and modes) with the computed tephra loadings obtained with the ARPAE database for the explored TGSDs (Fig. 5).</caption><thead><tr><th rowspan="2">Sample</th><th rowspan="2">Location</th><th colspan="4">Field observations</th><th colspan="4">Computed loading (kg m<sup>-2</sup>)</th></tr><tr><th>Longitude</th><th>Latitude</th><th>Mode (Φ)</th><th>Loading (kg m<sup>-2</sup>)</th><th>Field TGSD</th><th>Radar TGSD</th><th>Integrated TGSD</th><th>Whole TGSD</th></tr></thead><tbody><tr><td>CTL</td><td>Citelli</td><td>15.060</td><td>37.765</td><td>-3</td><td>1.7 x 10<sup>1</sup></td><td>7.2 x 10<sup>0</sup></td><td>3.4 x 10<sup>-1</sup></td><td>4.1 x 10<sup>0</sup></td><td>2.0 x 10<sup>0</sup></td></tr><tr><td>CRT</td><td>Cerrita</td><td>15.092</td><td>37.774</td><td>-2</td><td>1.4 x 10<sup>1</sup></td><td>5.2 x 10<sup>0</sup></td><td>3.5 x 10<sup>-1</sup></td><td>2.8 x 10<sup>0</sup></td><td>2.0 x 10<sup>0</sup></td></tr><tr><td>PDM</td><td>Piedimonte</td><td>15.177</td><td>37.810</td><td>-2</td><td>6.1 x 10<sup>0</sup></td><td>1.3 x 10<sup>1</sup></td><td>1.3 x 10<sup>-1</sup></td><td>6.6 x 10<sup>0</sup></td><td>1.8 x 10<sup>0</sup></td></tr><tr><td>FFD</td><td>Fiumefreddo</td><td>15.215</td><td>37.799</td><td>-1</td><td>1.6 x 10<sup>0</sup></td><td>9.6 x 10<sup>0</sup></td><td>2.9 x 10<sup>-1</sup></td><td>4.9 x 10<sup>0</sup></td><td>1.5 x 10<sup>0</sup></td></tr><tr><td>CPV</td><td>Campovolo</td><td>15.228</td><td>37.801</td><td>-2</td><td>9.5 x 10<sup>-1</sup></td><td>8.6 x 10<sup>0</sup></td><td>3.2 x 10<sup>-1</sup></td><td>4.4 x 10<sup>0</sup></td><td>1.4 x 10<sup>0</sup></td></tr><tr><td>GDN</td><td>Giardini</td><td>15.250</td><td>37.819</td><td>-1</td><td>4.0 x 10<sup>0</sup></td><td>9.8 x 10<sup>0</sup></td><td>3.8 x 10<sup>-1</sup></td><td>5.0 x 10<sup>0</sup></td><td>1.4 x 10<sup>0</sup></td></tr><tr><td>TER</td><td>T. Ellera</td><td>16.548</td><td>38.417</td><td>3</td><td>1.6 x 10<sup>-2</sup></td><td>4.0 x 10<sup>-4</sup></td><td>3.5 x 10<sup>-1</sup></td><td>1.5 x 10<sup>-2</sup></td><td>2.4 x 10<sup>-2</sup></td></tr></tbody></table>

major motivations for developing robust automated tools to forecast tephra loading and airborne ash dispersal (e.g. Costa et al., 2006; Barsotti et al., 2008; Folch et al., 2008, 2009). 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 tesselation 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, Eyjafjalajö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 satellite-based (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.
<figure><img src="image_5.png" /><figcaption>**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 Table 1.</figcaption></figure>

## 2 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).

<figure><img src="image_6.png" /><figcaption>**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).</figcaption></figure>

### 2.1 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),
<figure><img src="image_7.png" /><figcaption>**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.</figcaption></figure>

covering the range from −5 to 5Φ (where $d = 2^{-\Phi}$, 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 and Costa, 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 \pm 1.1 \times 10^9$ kg making use of the Weibull distribution method (Bonadonna and Costa, 2012, 2013). Then, combining the field-derived TEM with the paroxysmal duration (~50 min), they calculated an average MER of $4.5 \pm 3.6 \times 10^5$ kg s⁻¹. Furthermore, considering the climax phase only (i.e. from 09:55 to 10:14), the MER reached $10^6$ kg s⁻¹, 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.

## 2.2 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 (SEVIRI) on board the geostationary Meteosat Second Generation (MSG) satellite.

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Ash, ice, and SO₂ mass time series retrieved from SEVIRI for the 23 November 2013 Etna eruption.</figcaption></figure>

The Moderate Resolution Imaging Spectro-radiometer (MODIS) aboard the NASA-Aqua polar-orbit 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₂ and ash emission. SO₂ 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₂ 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
<figure><img src="image_9.png" /><figcaption>**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.</figcaption></figure>

TEM of ~3.0 × 10⁹ kg with a PM₂₀ 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 erup- tive 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 in- crease 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, 2016) at 3 km from the NSEC (La Montagnola station). Inferred radar parameters (e.g. backscattered echo power) are proportional to the quan- tity 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; Don- nadieu et al., 2017):

$$V_e = \frac{v_{r+}}{\sin\theta} \approx 3.89v_{r+}, \qquad (1)$$

where $V_e$ is the ejection velocities (in m s⁻¹), $v_{r+}$ is the ra- dial velocity (in m s⁻¹), and $\theta$ is the elevation angle of the radar beam (here $\theta = 14.9^{\circ}$). Such an approach is relevant for integrating the time-dependent ejection velocities with the corresponding observed eruptive column heights. In partic- ular, we used the VOLDORAD 2B data associated with the 23 November 2013 eruption to better constrain the eruption phases' characterization.

## 3 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 simula- tions are analysed in terms of tephra loadings and airborne ash mass dispersal to best fit the field and satellite measure- ments.

### 3.1 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 spec- trum (i.e. −5 to 5Φ; Sect. 2.1), we used the X-Radar re- trievals 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 $\Phi$. 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$\Phi$; 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 PM20 to reproduce the far-travelling airborne ash mass retrieved by satellite. We assessed the tail of the Integrated TGSD (i.e. $\Phi \ge 6$) by modifying empirically the PM20 fraction, adding mass into the corresponding classes. We calculated the fractions based on an empirical power-law dependence of the classes with $\Phi$ through the following parameterization:

$$X(\Phi_i) = X(\Phi_5) \times \gamma^{(\Phi_i - \Phi_5)}, \quad (2)$$

where $X(\Phi_i)$ is the fraction (in wt %) allocated to the ith bin, $X(\Phi_5)$ is the fraction obtained for $\Phi = 5$, and $\gamma$ is the empirical factor ($\gamma < 1$). The explored $\gamma$ values span from 0.1 to 0.7, giving, respectively, PM20 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.

## 3.2 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 – $\alpha$ and cross-flow – $\beta$; Bursik, 2001; Kaminski et al., 2005; Suzuki and Koyaguchi, 2015; Folch et al., 2016; Costa et al., 2016b). Here, $\alpha$ and $\beta$ 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 advection-diffusion-sedimentation equation computed over a terrain-following 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
<figure><img src="image_10.png" /><figcaption>**Figure 6.** The main meteorological profiles over the NSEC from ARPAE (INGV-OE) and ERA-Interim (ECMWF), and over Tirana for ERA-Interim.</figcaption></figure>

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⁻¹ over the NSEC (at 09:30) for ERA-Interim and ARPAE, respectively, and ~34 m s⁻¹ over Tirana at 14:00. Besides the velocities, the wind direction (Fig. 6) shows a strong north-easterly 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., 2014, 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 (Donadieu et al., 2015, 2016, 2017). The improved simulation scheme (Fig. 7b and c) is achieved by coupling this discretization with the ARPAE or ERA-Interim databases and
<figure><img src="image_11.png" /><figcaption>**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.</figcaption></figure>

the Integrated TGSD or Whole TGSD, respectively, depending on the inversion purpose.

## 3.3 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 goodness-of-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₁, RMSE₂, and RMSE₃) described

<table><caption>**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.</caption><thead><tr><th>Parameter</th><th colspan="2">Explored range</th></tr></thead><tbody><tr><td>TGSD</td><td colspan="2">Multiple</td></tr><tr><td>Column height</td><td colspan="2">Multiple</td></tr><tr><td>MER</td><td colspan="2">Multiple</td></tr><tr><td>Exit velocity</td><td colspan="2">Multiple</td></tr><tr><td>Initial magma temperature (°K)</td><td colspan="2">1300</td></tr><tr><td>Exit water fraction (wt %)</td><td colspan="2">2.5</td></tr><tr><td>Radial entrainment coefficient (α)</td><td>0.05</td><td>0.15</td></tr><tr><td>Cross-flow entrainment coefficient (β)</td><td>0.05</td><td>1.00</td></tr></tbody></table>

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

$$K = \exp \left[ \frac{1}{N} \sum_{i}^{N} \log \left( \frac{\text{Obs}_i}{\text{Sim}_i} \right) \right], \quad (3)$$

$$k = \exp \left[ \sqrt{ \frac{1}{N} \sum_{i}^{N} \log \left( \frac{\text{Obs}_i}{\text{Sim}_i} \right)^2 - \left( \frac{1}{N} \sum_{i}^{N} \log \left( \frac{\text{Obs}_i}{\text{Sim}_i} \right) \right)^2 } \right], \quad (4)$$

where *i* refers to the *i*th 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 PM20 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⁻²) retrieved from SEVIRI:

$$ \Delta Mass = \frac{1}{\Delta T} \int_{t_0}^{t_f} (M_{\text{Obs}} - M_{\text{Sim}}) dt, \qquad (5) $$

where $M_{\text{Obs}}$ and $M_{\text{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 $\Delta T = t_f - t_0$). This index gives the discrepancy (in t) for each $\gamma$ factor (i.e. PM20 fractions). Additionally, we also calculated for each $\gamma$ factor the absolute average difference of mass per unit area (Sum(Δ) in t km⁻²) for the entire volcanic cloud by the following:

$$ \overline{\text{Sum}(\Delta)} = \frac{1}{\Delta T} \int_{t_0}^{t_f} \frac{\sum_{N} |M_{\text{Obs}}(N) - M_{\text{Sim}}(N)|}{\text{Area}_\text{p}} dt, \qquad (6) $$

where N is the number of pixels (i.e. plume mask), and $M_{\text{Obs}}(N)$ and $M_{\text{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²). This index indicates the uncertainty of the simulated airborne ash mass per unit area with respect to the satellite retrieval.

Considering that ΔMass and $\overline{\text{Sum}(\Delta)}$ are discrepancy estimates, the selection is done on the basis of their minimization. Nonetheless, $\overline{\text{Sum}(\Delta)}$ gives absolute values preventing any over- or under-estimation characterization. It follows that we also evaluated the following index:

$$ \varepsilon = \frac{1}{\Delta T} \int_{t_0}^{t_f} \frac{\left[ \sum_{N} M_{\text{Obs}}(N) - M_{\text{Sim}}(N) < 0 \right] + \left[ \sum_{N} M_{\text{Obs}}(N) - M_{\text{Sim}}(N) > 0 \right]}{\text{Area}_\text{p}} dt, \qquad (7) $$

where $\varepsilon$ refers to an over-estimation per pixel when $\varepsilon < 0$ and an under-estimation per pixel for $\varepsilon > 0$, with a best fit for

$\varepsilon = 0$. Moreover, the index indicates the average mass difference per unit area (i.e. t km⁻²) 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.

# 4 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 PM20 fraction needed within the Whole TGSD to capture the airborne ash transported in distal area.

## 4.1 ESP inversion

Regarding the Integrated TGSD inversion (Sect. 3.3), Table 3 shows the statistical analysis for the best simulation (i.e. $K \approx 1$, RMSE₁, RMSE₂, RMSE₃, and k minimized) for each weighting average combination. Regardless of the weights, RMSE₁ and RMSE₃ have flat patterns, motivating us to rely on the RMSE₂ and k. They show relevant combinations from (65, 35; i.e. 65 and 35 in wt% for the Field and Radar TGSDs, respectively) to (85, 15). Although RMSE₂ ranges between 1.56 and 1.85 from (65, 35) to (85, 15), k is minimized at 2.95 for (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₂ and k indicate relatively high values yielding a mean error factor nearby 3, which is comparable to uncertainties associated with other classical methods (Bonadonna and Costa, 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 PM20 fraction. Considering the whole airborne ash mass, the results yield a best value for ΔMass at $\gamma = 0.65$ (i.e. PM20 = 9.0 wt%), indicating an overall underestimation of ~76t of ash by FALL3D for the entire eruption. Then, $\overline{\text{Sum}(\Delta)}$ shows a minimum for $\gamma = 0.40$ (i.e. PM20 = 3.6 wt%), giving an absolute average difference of mass per unit area of ~0.37 t km⁻² for the whole sequence. The third index returns a best value of $\varepsilon = -0.03$ t km⁻² for $\gamma = 0.65$ (i.e. PM20 = 9.0 wt%), being consistent with ΔMass. $\varepsilon$ likely reflects that FALL3D slightly over-estimates the average mass per pixel of 0.03 t km⁻². By integrating the results (Fig. 9), the Whole TGSD required the minimum PM20 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 under-estimation, which tends to be satisfied with higher PM20 fractions (i.e. 9.0 wt%). The corresponding input TGSD is
<table><caption>**Table 3.** Statistical metric for the best simulations (i.e. calibration of α and β) for each weighting average combination tested during the inversion of the Integrated TGSD. α is described through α₁ and α₂ within the calculation (Folch et al., 2016). TEM indicates the associated theoretical value for each combination.</caption><thead><tr><th rowspan="2">Integrated TGSD<br>Combination<br>(in wt%)</th><th colspan="2">Input</th><th colspan="8">Statistical metric</th><th>Output</th></tr><tr><th>α<br>(α₁-α₂)</th><th>β</th><th>K</th><th>k</th><th>RMSE₁</th><th>RMSE₂</th><th>RMSE₃</th><th>Correlation</th><th>Bias</th><th>T test</th><th>TEM<br>(× 10⁹ in kg)</th></tr></thead><tbody><tr><td>Radar TGSD</td><td>0.15-0.15</td><td>1.00</td><td>6.97</td><td>9.82</td><td>0.97</td><td>7.71</td><td>0.87</td><td>-0.2</td><td>-0.3</td><td>0.1</td><td>5.73</td></tr><tr><td>20 Field/80 Radar</td><td>0.06-0.09</td><td>0.72</td><td>1.00</td><td>4.35</td><td>0.84</td><td>2.95</td><td>0.74</td><td>0.8</td><td>0.0</td><td>1.0</td><td>2.8</td></tr><tr><td>40 Field/60 Radar</td><td>0.06-0.09</td><td>0.40</td><td>1.02</td><td>3.48</td><td>0.81</td><td>1.61</td><td>0.74</td><td>0.8</td><td>0.0</td><td>1.0</td><td>1.66</td></tr><tr><td>60 Field/40 Radar</td><td>0.06-0.09</td><td>0.28</td><td>1.01</td><td>3.08</td><td>0.78</td><td>1.53</td><td>0.77</td><td>0.9</td><td>0.0</td><td>1.0</td><td>1.28</td></tr><tr><td>65 Field/35 Radar</td><td>0.06-0.09</td><td>0.26</td><td>1.01</td><td>3.02</td><td>0.77</td><td>1.56</td><td>0.77</td><td>0.9</td><td>0.0</td><td>1.0</td><td>1.22</td></tr><tr><td>70 Field/30 Radar</td><td>0.06-0.09</td><td>0.25</td><td>0.98</td><td>2.98</td><td>0.77</td><td>1.67</td><td>0.80</td><td>0.9</td><td>0.0</td><td>1.0</td><td>1.21</td></tr><tr><td>75 Field/25 Radar</td><td>0.06-0.09</td><td>0.22</td><td>1.02</td><td>2.95</td><td>0.76</td><td>1.64</td><td>0.79</td><td>0.9</td><td>0.0</td><td>1.0</td><td>1.13</td></tr><tr><td>80 Field/20 Radar</td><td>0.06-0.09</td><td>0.22</td><td>0.99</td><td>2.96</td><td>0.75</td><td>1.77</td><td>0.82</td><td>0.9</td><td>0.0</td><td>1.0</td><td>1.13</td></tr><tr><td>85 Field/15 Radar</td><td>0.06-0.09</td><td>0.21</td><td>1.01</td><td>3.00</td><td>0.75</td><td>1.85</td><td>0.84</td><td>0.9</td><td>0.0</td><td>1.0</td><td>1.10</td></tr><tr><td>90 Field/10 Radar</td><td>0.06-0.09</td><td>0.21</td><td>1.00</td><td>3.13</td><td>0.74</td><td>2.02</td><td>0.88</td><td>0.9</td><td>0.0</td><td>1.0</td><td>1.12</td></tr><tr><td>Field TGSD</td><td>0.06-0.09</td><td>0.35</td><td>0.99</td><td>6.56</td><td>0.83</td><td>3.65</td><td>1.44</td><td>0.9</td><td>0.0</td><td>1.0</td><td>1.60</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**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 Table 2).</figcaption></figure>

displayed in Fig. 5. Moreover, ΔMass and ε in Fig. 9 both indicate that FALL3D under-estimates substantially the airborne mass for PM₂₀ fractions lower than ~7 wt% and over-estimates 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⁵ kg s⁻¹. 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).

## 4.2 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-
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Quantitative analysis of the airborne ash mass measured from SEVIRI and computed by FALL3D to invert the PM<sub>20</sub> 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).</figcaption></figure>

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<sup>-2</sup>. In contrast, FALL3D computed them between 3 and 7 kg m<sup>-2</sup> 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 field-derived 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<sub>1</sub> = 0.80, RMSE<sub>2</sub> = 3.36, and RMSE<sub>3</sub> = 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<sup>9</sup> kg, which is ~34% higher than for the integrated approach with 1.2 × 10<sup>9</sup> kg. The latter TEM is in good agreement with the estimation of 1.3 × 10<sup>9</sup> kg reported in Andronico et al. (2015). It is worth noting that varying the weighting av-
<figure><img src="image_14.png" /><figcaption>**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.</figcaption></figure>

erage from 100 wt % Field TGSD towards 100 wt % Radar TGSD yields an increasing TEM going from 1 to $6 \times 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 \times 10^9$ kg compared to the field-derived TEM of $1.3 \times 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\Phi$), while the field sampling method expands the measurements to block-sized ($-5$ to $5\Phi$); (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><img src="image_15.png" /><figcaption>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.</figcaption></figure>

## 4.3 Airborne ash dispersal

As mentioned in Sect. 2, large quantities of ash, water vapour (transformed into ice), and SO₂ 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₂₀ range (i.e. 3.6–9.0 wt %), we displayed the airborne ash mass maps. The left column refers to the minimum PM₂₀ 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 SEVIRI 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₂₀ = 9.0 wt %) indicate a better performance than for γ = 0.40 (i.e. PM₂₀ = 3.6 wt %). The entire time-series an-
imations are available in the Supplement (Animations A3 and A4 for $\gamma = 0.40$ and $\gamma = 0.65$, respectively).

The PM20 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 PM20 fractions indicate consistent values. Indeed, 1–2 wt% of the X-Radar TEM (3.0 × 10⁹ kg) refers to 30–60 t, while 3.6–9.0 wt% of the integrated TEM (1.2 × 10⁹ kg) gives 43–108 t. In fact, Corradini et al. (2016) integrated X-Radar data with satellite retrievals to assess the PM20 fraction. However, the satellite cannot quantify any ash mass from pixels mainly filled by ice or gas (e.g. SO₂). 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.0 wt% 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. (2016a, 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. PM20, 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).

## 5 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 grain-size 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. PM20). We empirically modified the Integrated TGSD to implement the SEVIRI retrievals by investigating diverse PM20 fractions (i.e. 0.6–10.7 wt%) until we best fit the measurements. The inverted PM20 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. PM20, 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).
Appendix A

Appendix completes Table 2 in terms of parameterizations (i.e. parameters and models) used to run the simulations under the ARPAE and ERA-Interim meteorological databases.

<table><caption>**Table A1.** Additional parameterizations used in the simulations.</caption><thead><tr><th>Parameterization</th><th>ARPAE</th><th>ERA-Interim</th></tr></thead><tbody><tr><td>Vent elevation (m a.s.l.)</td><td>3300</td><td>3300</td></tr><tr><td>Vent longitude (°)</td><td>15.002012</td><td>15.002012</td></tr><tr><td>Vent latitude (°)</td><td>37.746548</td><td>37.746548</td></tr><tr><td>Time step meteo data (min)</td><td>30</td><td>30</td></tr><tr><td>Longitude nodes</td><td>160</td><td>115</td></tr><tr><td>Latitude nodes</td><td>100</td><td>100</td></tr><tr><td>Grid resolution (km²)</td><td>1</td><td>5</td></tr><tr><td>Altitude layers (from 0 m a.s.l., 500 m step)</td><td>12 000</td><td>12 000</td></tr><tr><td>Eruption column model</td><td>FPlume<sup>a</sup></td><td>FPlume<sup>a</sup></td></tr><tr><td>Terminal velocity model</td><td>Ganser<sup>b</sup></td><td>Ganser<sup>b</sup></td></tr><tr><td>Vertical turbulence model</td><td>Similarity<sup>c</sup></td><td>Similarity<sup>c</sup></td></tr><tr><td>Horizontal turbulence model</td><td>CMAQ<sup>d</sup></td><td>CMAQ<sup>d</sup></td></tr><tr><td>Gravity current</td><td>Yes<sup>e</sup></td><td>Yes<sup>e</sup></td></tr></tbody></table>

<footnote><sup>a</sup> The eruption column model uses the buoyant plume theory (Folch et al., 2016). <sup>b</sup> The terminal settling velocity is calculated through the Ganser model (Ganser, 1993). <sup>c</sup> The vertical component of the eddy diffusivity tensor ($K_{\perp}$) is estimated using the similarity option (Costa et al., 2006: Ulke, 2000). <sup>d</sup> The horizontal component of the eddy diffusivity tensor ($K_H$) is evaluated as in Byun and Schere (2006) by the CMAQ option. <sup>e</sup> The gravity current effects in the umbrella region are negligible in the far-field region, but were considered in the simulations (Costa et al., 2013; Suzuki and Koyaguchi, 2009).</footnote>
## 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 $\gamma = 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 $\gamma = 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 far-travelling airborne ash dispersal computed with the Whole TGSD for $\gamma = 0.65$ (i.e. 09:30–14:30 UTC). This animation shows a similar dispersal to Animation A3. However, using $\gamma = 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.

**The Supplement related to this article is available online at https://doi.org/10.5194/acp-18-4695-2018-supplement.**

*Author contributions.* MP conceived the idea and defined the project approach together with AC, SC, and LM. MP and AC com-

ported the FALL3D simulations and co-wrote the manuscript. SC and LM provided and processed the SEVIRI data. GV and MM provided and processed the X-band weather radar data. DA provided and processed the field data. VF-L provided and processed the L-band VOLDORAD 2B data. All the co-authors worked on the interpretation of the results and finalization of the manuscript.

*Competing interests.* The authors declare that they have no conflict of interest.

*Acknowledgements.* Matthieu Poret received funding from FP7 Marie Curie Actions Framework “FP7-PEOPLE-2013-ITN”, project VERTIGO grant number 607905. Antonio Costa and Daniele Andronico acknowledge European project EUROVOLC (grant number 731070) and MIUR project Premiale Ash-RESILIENCE. The radar is operated jointly by the OPGC and INGV-OE (Catania, Italy) in the framework of a collaboration agreement between INGV-OE, the French CNRS, and the OPGC-Université Clermont Auvergne in Clermont-Ferrand (France). The X-band weather radar data were provided by the Civil Protection Department (Rome, Italy) and the MSG-SEVIRI data by the INGV in Rome. We are also grateful to Boris Behncke (INGV-OE) for the imagery support of the eruption. We warmly thank Nicola and Salvatore Costa for collecting the sample at T. Ellera. We also warmly acknowledge Massimo Cantanero, Rosa Anna Corsaro, and Antonio Cristaldi, who helped to collect the tephra samples and analyse them. Finally, we are grateful to Larry Mastin and an anonymous reviewer for their valuable comments that improved the quality and clarity of the manuscript.

Edited by: Johannes Quaas

Reviewed by: Larry Mastin and one anonymous referee

## References

Aida, I.: Reliability of a tsunami source model derived from fault parameters, J. Phys. Earth, 26, 57–73, 1978.

Alparone, S., Andronico, D., Lodato, L., and Sgroi, T.: Relationship between tremor and volcanic activity during the Southeast Crater eruption on Mount Etna in early 2000, J. Geophys. Res., 108, 2241, https://doi.org/10.1029/2002JB001866, 2003.

Anderson, K. and Segall, P.: Bayesian inversion of data from effusive volcanic eruptions using physics-based models: Application to Mount St. Helens 2004–2008, J. Geophys. Res.-Sol. Ea., 118, 2017–2037, https://doi.org/10.1002/jgrb.50169, 2013.

Andronico, D. and Del Carlo, P.: PM$_{10}$ measurements in urban settlements after lava fountain episodes at Mt. Etna, Italy: Pilot test to assess volcanic ash hazard on human health, Nat. Hazards Earth Syst. Sci., 16, 29–40, https://doi.org/10.5194/nhess-16-29-2016, 2016.

Andronico, D., Scollo, S., Cristaldi, A., and Lo Castro, M. D.: Representivity of incompletely sampled fall deposits in estimating eruption source parameters: a test using the 12–13 January 2011 lava fountain deposit from Mt. Etna volcano, Italy, Bull. Volcanol., 76, 861, https://doi.org/10.1007/s00445-014-0861-3, 2014.
Andronico, D., Scollo, S., and Cristaldi, A.: Unexpected hazards from tephra fallouts at Mt. Etna: The 23 November 2013 lava fountain, J. Volcanol. Geotherm. Res., 304, 118–125, https://doi.org/10.1016/j.jvolgeores.2015.08.007, 2015.

Barsotti, S., Neri, A., and Scire, J.: The VOL-CALPUFF model for atmospheric ash dispersal. I. Approach and physical formulation, J. Geophys. Res., 113, B03208, https://doi.org/10.1029/2006JB004623, 2008.

Beckett, F. M., Witham, C. S., Hort, M. C., Stevenson, J. A., Bonadonna, C., and Millington, S. C.: Sensitivity of dispersion model forecasts of volcanic ash clouds to the physical characteristics of the particles, J. Geophys. Res.-Atmos., 120, 11636–11652, https://doi.org/10.1002/2015JD023609, 2015.

Behncke, B., Branca, S., Corsaro, R. A., De Beni, E., Miraglia, L., and Proietti, P.: The 2011–2012 summit activity of Mount Etna: birth, growth and products of the new SE crater, J. Volcanol. Geotherm. Res., 270, 10–21, https://doi.org/10.1016/j.jvolgeores.2013.11.012, 2014.

Bonaccorso, A., Calvari, S., Linde, A., and Sacks, S.: Eruptive processes leading to the most explosive lava fountain at Etna volcano: The 23 November 2013 episode, Geophys. Res. Lett., 41, 4912–4919, https://doi.org/10.1002/2014GL060623, 2014.

Bonadonna, C. and Costa, A.: Estimating the volume of tephra deposits: a new simple strategy, Geology, 40, 415–418, https://doi.org/10.1130/G32769.1, 2012.

Bonadonna, C. and Costa, A.: Plume height, volume, and classification of explosive volcanic eruptions based on the Weibull function, Bull. Volcanol., 75, 742, https://doi.org/10.1007/s00445-013-0742-1, 2013.

Bonadonna, C. and Houghton, B. F.: Total grain-size distribution and volume of tephra-fall deposits, Bull. Volcanol., 67, 441–456, https://doi.org/10.1007/s00445-004-0386-2, 2005.

Bonadonna, C., Genco, R., Gouhier, M., Pistolesi, M., Cioni, R., Alfano, F., Hoskuldsson, A., and Ripepe, M.: Tephra sedimentation during the 2010 Eyjafjallajökull eruption (Iceland) from deposit, radar, and satellite observations, J. Geophys. Res.-Sol. Ea., 116, B12202, https://doi.org/10.1029/2011JB008462, 2011.

Bonadonna, C., Biass, S., and Costa, A.: Physical characterization of explosive volcanic eruptions based on tephra deposits: Propagation of uncertainties and sensitivities analysis, J. Volcanol. Geotherm. Res., 296, 80–100, https://doi.org/10.1016/j.jvolgeores.2015.03.009, 2015.

Bursik, M.: Effect of wind on the rise height of volcanic plumes, Geophys. Res. Lett., 28, 3621–3624, https://doi.org/10.1029/2001GL013393, 2001.

Byun, D. and Schere, K. L.: Review of the governing equations, computational algorithms, and other components of the Models-3 Community Multiscale Air Quality (CMAQ) modeling system, Appl. Mech. Rev., 59, 51–77, https://doi.org/10.1115/1.2128636, 2006.

Carbone, D., Zuccarello, L., Messina, A., Scollo, S., and Rymer, H.: Balancing bulk gas accumulation and gas output before and during lava fountaining episodes at Mt. Etna, Sci. Rep., 5, 18049, https://doi.org/10.1038/srep18049, 2015.

Cerminara, M., Esposti Ongaro, T., and Neri, A.: Large eddy simulation of gos-particle kinematic decoupling and turbulent entrainment in volcanic plumes, J. Volcanol. Geotherm. Res., 326, 143–171, https://doi.org/10.1016/j.jvolgeores.2016.06.018, 2016.

Corradini, S., Montopoli, M., Guerrieri, L., Ricci, M., Scollo, S., Merucci, L., Marzano, F. S., Pugnaghi, S., Prestifilippo, M., Ventress, L., Grainger, R. G., Carboni, E., Vulpiani, G., and Coltelli, M.: A multi-sensor approach for the volcanic ash cloud retrievals and eruption characterization, Remote Sensing, 8, 58, https://doi.org/10.3390/rs8010058, 2016.

Costa, A., Macedonio, G., and Folch, A.: A three-dimensional Eulerian model for transport and deposition of volcanic ashes, Earth Planet. Sc. Lett., 241, 634–647, https://doi.org/10.1016/j.epsl.2005.11.019, 2006.

Costa, A., Folch, A., and Macedonio, G.: A model for wet aggregation of ash particles in volcanic plumes and clouds: I. Theoretical formulation, J. Geophys. Res., 115, B09201, https://doi.org/10.1029/2009JB007175, 2010.

Costa, A., Folch, A., and Macedonio, G.: Density-driven transport in the umbrella region of volcanic clouds: implications for tephra dispersion models, Geophys. Res. Lett., 40, 4823–4827, https://doi.org/10.1002/grl.50942, 2013.

Costa, A., Pioli, L., and Bonadonna, C.: Assessing tephra total grain-size distribution: Insights from field data analysis, Earth Planet. Sc. Lett., 443, 90–107, 2016a.

Costa, A., Suzuki, Y. J., Cerminara, M., Devenish, B., Esposti Ongaro, T., Herzog, M., Van Eaton, A., Denby, L. C., Bursik, M., de Michieli Vitturi, M., Engwell, S., Neri, A., Barsotti, S., Folch, A., Macedonio, G., Girault, F., Carazzo, G., Tait, S., Kaminski, E., Mastin, L., Woodhouse, M., Phillips, J. C., Hogg, A. J., Degruyter, W. J., and Bonadonna, C.: Results of the eruption column model intercomparison study, J. Volcanol. Geotherm. Res., 326, 2–25, https://doi.org/10.1016/j.jvolgeores.2016.01.017, 2016b.

Costa, A., Pioli, L., and Bonadonna, C.: Corrigendum to “Assessing tephra total grain-size distribution: Insights from field data analysis”. (Earth and Planetary Sci. Lett. 443, 90–107, 2016), Earth Planet. Sc. Lett., 465, 205–209, https://doi.org/10.1016/j.epsl.2017.03.003, 2017.

Crompton, E. and Husson, P.: Volcanic Ash Advisories: how the VAACs use the ‘Discernible Ash’ definition to draw their lines now and in the future, in: 7th International Workshop on Volcanic Ash (IWVA/7), 19–23 October 2015, Anchorage, Alaska, USA, https://www.wmo.int/aemp/sites/default/files/P-09_Crompton-Husson.pdf (last access: 4 April 2018), 2015.

Dacre, H. F., Grant, A. L. M., Hogan, R. J., Belcher, S. E., Thomson, D. J., Devenish, B., Marenco, F., Haywood, J., Ansmann, A., and Mattis, I.: The structure and magnitude of the ash plume during the initial phase of the Eyjafjallajökull eruption, evaluated using lidar observations and NAME simulations, J. Geophys. Res., 116, D00U03, https://doi.org/10.1029/2011JD015608, 2011.

De Beni, E., Behncke, B., Branca, S., Nicolosi, I., Carluccio, R., D’Ajello Caracciolo, F., and Chiappini, M.: The continuing story of Etna’s New Southeast Crater (2012–2014): Evolution and volume calculations based on field surveys and aerophotogrammetry, J. Volcanol. Geotherm. Res., 303, 175–186, https://doi.org/10.1016/j.jvolgeores.2015.07.021, 2015.

Degruyter, W. and Bonadonna, C.: Improving on mass flow rate estimates of volcanic eruptions, Geophys. Res. Lett., 39, L16308, https://doi.org/10.1029/2012GL052566, 2012.

Devenish, B. J., Rooney, G. G., Webster, H. N., and Thomson, D. J.: The entrainment rate for buoyant plumes in a crossflow, Bound-
Lay. Meteorol., 134, 411-439, https://doi.org/10.1007/s10546-009-9464-5, 2010.

Devenish, B., Francis, P. N., Johnson, B. T., Sparks, R. S. J., and Thomson, D. J.: Sensitivity analysis of dispersion modeling of volcanic ash from Eyjafjallajökull in May 2010, J. Geophys. Res., 117, D00U21, https://doi.org/10.1029/2011JD016782, 2012.

Donnadieu, F., Freville, P., Rivet, S., Hervier, C., and Cacault, P.: The Volcano Doppler radar data base of Etna (VOLDORAD 2B), Université Clermont Auvergne - CNRS, https://doi.org/10.18145/VOLDORAD.ETNA.2009, 2015.

Donnadieu, F., Freville, P., Hervier, C., Coltelli, M., Scollo, S., Prestifilipo, M., Valade, S., Rivet, S., and Cacault, P.: Near-source Doppler radar monitoring of tephra plumes at Etna, J. Volcanol. Geotherm. Res., 312, 26-39, 2016.

Donnadieu, F., Freret-Lorgeril, V., Coltelli, M., Scollo, S., Gouhier, M., Fréville, P., Hervier, C., and Prestifilippo, M.: Lava fountaining paroxysms generating tephra plumes at Mt. Etna: Remote sensing retrievals, in: IAVCEI meeting, Poster VO43A-139, 14-18 August 2017, Portland, USA, 2017.

Folch, A.: A review of tephra transport and dispersal models: Evolution, current status, and future perspectives, J. Volcanol. Geotherm. Res., 235-236, 96-115, https://doi.org/10.1016/j.jvolgeores.2012.05.020, 2012.

Folch, A., Cavazzoni, C., Costa, A., and Macedonio, G.: An automatic procedure to forecast tephra fall-out, J. Volcanol. Geotherm. Res., 177, 767-777, https://doi.org/10.1016/j.jvolgeores.2008.01.046, 2008.

Folch, A., Costa, A., and Macedonio, G.: FALL3D: A computational model for transport and deposition of volcanic ash, Comput. Geosci., 35, 1334-1342, https://doi.org/10.1016/j.cageo.2008.08.008, 2009.

Folch, A., Costa, A., Durant, A., and Macedonio, G.: A model for wet aggregation of ash particles in volcanic plumes and clouds: 2. Model application, J. Geophys. Res., 115, B09202, https://doi.org/10.1029/2009JB007176, 2010.

Folch, A., Costa, A., and Basart, S.: Validation of the FALL3D ash dispersion model using observations of the 2010 Eyjafjallajökull volcanic ash clouds, Atmos. Environ., 48, 165-183, https://doi.org/10.1016/j.atmosenv.2011.06.072, 2012.

Folch, A., Costa, A., and Macedonio, G.: FPLUME-1.0: An integral volcanic plume model accounting for ash aggregation, Geosci. Model Dev., 9, 431-450, https://doi.org/10.5194/gmd-9-431-2016, 2016.

Freret-Lorgeril, V., Donnadieu, F., Coltelli, M., Scollo, S., Fréville, P., Hervier, C., and Prestifilippo, M.: Doppler radar retrievals from lava fountaining paroxysms generating tephra plumes at Mt. Etna, EGU General Assembly, 23-28 April 2016, Vienna, Austria, 2016.

Ganser, G. H.: A rational approach to drag prediction of spherical and non-spherical particles, Powder Technol., 77, 143-152, https://doi.org/10.1016/0032-5910(93)80051-B, 1993.

Guffanti, M., Casadevall, T. J., and Budding, K.: Encounters of Aircraft with Volcanic Ash Clouds: A Compilation of Known Incidents, 1953-2009, US Geological Survey Data Series 545, ver. 1.0, p. 12, plus 4 appendixes including the compilation database, available at: http://pubs.usgs.gov/ds/545 (last access: 4 April 2018), 2010.

Horwell, C. J., Sargent, P., Andronico, D., Lo Castro, M. D., Tomatis, M., Hillman, S. E., Michnowicz, S. A. K., and Fubini, B.: The iron-catalysed surface reactivity and health-pertinent physical characteristics of explosive volcanic ash from Mt. Etna, Italy, J. Appl. Volcanol., 6, https://doi.org/10.1186/s13617-017-0063-8, 2017.

Kaminski, E., Tait, S., and Carazzo, G.: Turbulent entrainment in jets with arbitrary buoyancy, J. Fluid Mech., 526, 361-376, https://doi.org/10.1017/S0022112004003209, 2005.

Macedonio, G., Costa, A., and Folch, A.: Uncertainties in volcanic plume modelling: A parametric study using FPlume, J. Volcanol. Geotherm. Res., 326, 92-102, https://doi.org/10.1016/j.jvolgeores.2016.03.016, 2016.

Mastin, L. G., Guffanti, M., Servranckx, R., Webley, P., Barsotti, S., Dean, K., Durant, A., Ewert, J. W., Neri, A., Rose, W. I., Schneidner, D., Siebert, L., Stunder, B., Swanson, G., Tupper, A., Volentik, A., and Waythomas, C. F.: A multidisciplinary effort to assign realistic source parameters to models of volcanic ash-cloud transport and dispersion during eruptions, J. Volcanol. Geotherm. Res., 186, 10-21, 2009.

Montopoli, M.: Velocity profiles inside volcanic clouds from three-dimensional scanning microwave dual-polarization Doppler radars, J. Geophys. Res.-Atmos., 121, 7881-7900, https://doi.org/10.1002/2015JD023464, 2016.

Morton, B. R., Taylor, G. I., and Turner, J. S.: Turbulent gravitational convection from maintained and instantaneous source, Philos. T. Roy. Soc. Lond. A, 234, 1-23, 1956.

Poret, M., Costa, A., Folch, A., and Martì, A.: Modelling tephra dispersal and ash aggregation: The 26th April 1979 eruption, La Soufrière St. Vincent, J. Volcanol. Geotherm. Res., 347C, 207-220, https://doi.org/10.1016/j.jvolgeores.2017.09.012, 2017.

Prata, A. J. and Kerkmann, J.: Simultaneous retrieval of volcanic ash and SO₂ using MSG-SEVIRI measurements, Geophys. Res. Lett., 34, L05813, https://doi.org/10.1029/2006GL028691, 2007.

Rose, W. I. and Durant, A. J.: Fine ash content of explosive eruptions, J. Volcanol. Geoth. Res., 186, 32-39, https://doi.org/10.1016/j.jvolgeores.2009.01.010, 2009.

Scollo, S., Folch, A., and Costa, A.: A parametric and comparative study of different tephra fallout models, J. Volcanol. Geoth. Res., 176, 199-211, https://doi.org/10.1016/j.jvolgeores.2008.04.002, 2008.

Scollo, S., Prestifilippo, M., Pecora, E., Corradini, S., Merucci, L., Spata, G., and Coltelli, M.: Eruption column height estimation of the 2011-2013 Etna lava fountains, Ann. Geophys., 57, S0214, https://doi.org/10.4401/ag-6396, 2014.

Scollo, S., Boselli, A., Coltelli, M., Leto, G., Pisani, G., Prestifilippo, M., Spinelli, N., and Wang, X.: Volcanic ash concentration during the 12 August 2011 Etna eruption, Geophys. Res. Lett., 42, 2634-2641, https://doi.org/10.1002/2015GL063027, 2015.

Spanu, A., dé Vitturi, M. M., and Barsotti, S.: Reconstructing eruptive source parameters from tephra deposit: a numerical study of medium-sized explosive eruptions at Etna volcano, Bull. Volcanol., 78, 59, https://doi.org/10.1007/s00445-016-1051-2, 2016.

Spilliaert, N., Allard, P., Metrich, N., and Sobolev, A. V.: Melt inclusion record of the conditions of ascent, degassing, and extrusion of volatile-rich alkali basalt during the powerful 2002 flank eruption of Mount Etna (Italy), J. Geophys. Res., 111, B04203, https://doi.org/10.1029/2005JB003934, 2006.
Stunder, B., Heffter, J. L., and Draxler, R. R.: Airborne volcanic ash forecast reliability, Weather Forecast., 22, 1132–1139, 2007.

Sulpizio, R., Folch, A., Costa, A., Scaini, C., and Dellino, P.: Hazard assessment of far-range volcanic ash dispersal from a violent Strombolian eruption at Somma-Vesuvius volcano, Naples, Italy: Implications on civil aviation, Bull. Volcanol., 74, 2205–2218, https://doi.org/10.1007/s00445-012-0656-3, 2012.

Suzuki, Y. J. and Koyaguchi, T.: A three-dimensional numerical simulation of spreading umbrella clouds, J. Geophys. Res., 114, B03209, https://doi.org/10.1029/2007JB005369, 2009.

Suzuki, Y. J. and Koyaguchi, T.: Effects of wind on entrainment efficiency in volcanic plumes, J. Geophys. Res.-Sol. Ea., 61, 6122–6140, https://doi.org/10.1002/2015JB012208, 2015.

Tomašek, I., Horwell, C. J., Damby, D. E., Barošová, H., Geers, C., Petri-Fink, A., Rothen-Rutishauser, B., and Clift, M. J. D.: Combined exposure of diesel exhaust particles and respirable Soufrière Hills volcanic ash causes a (pro-) inflammatory response in an in vitro multicellular epithelial tissue barrier model, Part. Fibre Toxicol., 13, 67, https://doi.org/10.1186/s12989-016-0178-9, 2016.

Ulke, A. G.: New turbulent parameterization for a dispersion model in atmospheric boundary layer, Atmos. Environ., 34, 1029–1042, 2000.

Vulpiani, G., Ripepe, M., and Valade, S.: Mass discharge rate retrieval combining weather radar and thermal camera observations, J. Geophys. Res.-Sol. Ea., 121, 5679–5695, https://doi.org/10.1002/2016JB013191, 2016.

Wen, S. and Rose, W. I.: Retrieval of sizes and total masses of particles in volcanic clouds using AVHRR bands 4 and 5, J. Geophys. Res., 99, 5421–5431, 1994.

Witham, C. S., Hort, M. C., Potts, R., Servranckx, R., Husson, P., and Bonnardot, F.: Comparison of VAAC atmospheric dispersion models using the 1 November 2004 Grimsvötn eruption, Meteorol. Appl., 14, 27–38, https://doi.org/10.1002/met.3, 2007.

Woodhouse, M. J., Hogg, A. J., Phillips, J. C., and Sparks, R. S. J.: Interaction between volcanic plumes and wind during the 2010 Eyjafjallajökull eruption, Iceland, J. Geophys. Res.-Sol. Ea., 118, 92–109, https://doi.org/10.1029/2012JB009592, 2013.