Evaluating an average density for the NSEC cone, instead, is a major challenge due to its rather high heterogeneity, since it is constituted mostly by tephra ranging in grain-size from large clasts (up to meter-sized bombs) to lapilli and ash particles (with varying degrees of compaction and welding) and minor intercalated lava flows. For the purposes of our work, we collected 10 decimeter-sized fragments taken from a portion of the crater rim that collapsed and was rafted downslope during the effusive activity (**Figure 6a**), and we assumed them to be representative of the cone for providing density measurements. These fragments were dried, weighted, waterproofed with plastic wraps and finally immersed in water, resulting in weights and volumes ranging between 113 and 857 g, and 120 and 700 cm3 , respectively, which is equal to clasts from 8 to 18 cm (in dimension). By measuring the excess water volume, we obtained a mean value of density of 1.22 g/cm3 , providing a vesicularity of ∼55% for the cone. This value agrees well enough with the density value estimated by Mulas et al. (2016) for the 1669 Monti Rossi cone, considering it as consisting of coarse lapilli and bombs (without lava intercalations), of 1.8 g/cm3 . So far, in this paper we adopted the value of 50% in order to compare our volume estimation with the value used by Behncke et al. (2014) and De Beni et al. (2015) to evaluate the growth of the NSEC in the period 2011–13. Due to the irregular packing of volcanic particles (strongly varying in size and density) which compose the NSEC, we believe that this value can be considered reasonably realistic though conservative, unless new estimates of volcanic cone density will be published.

#### RESULTS

## Growth of the Cone During Paroxysmal Events

The cone of the NSEC grew over a pre-existing subsidence depression formed between 2007 and 2009 and significantly enlarged due to minor explosive activity in 2010 (Andronico et al., 2013); its growth has been then documented by repeated GPS surveys starting from 2011 (Behncke et al., 2014). The NSEC at the time of our study was mostly composed of large welded scoriaceous bombs and blocks; on 26 October, some sizeable blocks (3–5 m) were torn from the crater rim and carried along by the southern lava flow close to the 2002–03 cones (**Figures 6a,b**). These blocks permitted insights into the structure of the inner portion of the crater rim, consisting of scoria, spatter, and ash layers of reddish-brown color. In **Table 1** we reported the values of the cone volumes for the 2011–12 and 2013 sequences. The volume of the cone was estimated at 19 × 106 m3 by Behncke et al. (2014) after the 2011–12 lava fountaining episodes, i.e., 0.76 × 106 m3 per episode. Furthermore, the difference between the May 2012 and October 2014 DEMs (De Beni et al., 2015), when 22 paroxysmal events and a few Strombolian-effusive episodes took place, provides a volume increase for the NSEC of 27.0 ± 0.8 ×106 m3 . However, considering that the base of the cone

TABLE 2 | Volcanological parameters of the 25–26 October 2013 lava fountain episode.

We reported: the bulk volume increase of the cone, the bulk lava volume, and the dense rock equivalent (DRE) volume plus the correspondent magma eruption rate (MGER) for both cone (assuming 50 % vesicularity), tephra fallout (30 %) and lava flow (20 %) deposits. The value marked by the apex (\*) has been assumed by averaging the volume increases of the 2011-13 episodes. The average DRE MGER value from effusive and explosive deposits is 69.63 m3 /s.

TABLE 1 | Volcanological parameters of the NSEC cone during the 2011–12 and 2013 lava fountaining sequences.

We reported: the bulk values and the correspondent dense rock equivalent (DRE) values (assuming 50 % vesicularity) of the whole cone volume, the average volume increase for a single episode of each sequence, the average values of the two sequences, and the magma eruption rate (MGER) of the DRE total cone volume (single sequence and average of the sum of the two cones).

was smaller on 26 October 2013 than the base on October 2014, as deduced from field observation and GPS measurements, we calculate a volume of 23.28 ±1.57× 106 m3 for the scoria cone in October 2013. By subtracting the volume of the lava flows emplaced over the cone in that period, we obtained a total volume for the NSEC of 20.48 ± 1.69 × 106 m3 . Dividing this value per the number of events, we derived an average increase of the cone volume of 1.08 ± 0.08 x 106 m3 for each paroxysm. By averaging the previous two estimates, we assume that, from 2011 to 2013 the NSEC has grown by 0.90 × 106 m3 during each paroxysm, a value we consider representative also for the volume increase caused by the 25-26 October 2013 episode (**Table 2**).

#### Distal Tephra Fallout Deposit

During the paroxysmal phase, on 26 October, tephra fallout extended to at least 130 km distance from the vent, beyond the southern coast of Sicily (**Figure 7a**). The dispersal area of the tephra was quite symmetrical and the deposit well-visible on the ground due to almost constant direction of wind above the volcano from ENE to WSW during the eruption. The load on unit area declined very rapidly orthogonally to the dispersal axis. Adrano, Centuripe (∼10 km WSW from Adrano), Caltanissetta and, to a lesser degree, Enna were the towns most affected by tephra fallout (principally ash); small amounts of fine-grained ash reached also the urban settlements of Leonforte and Montedoro.

Along the dispersal axis, in the proximal area (at ∼1.4 km from the NSEC), the fallout deposit formed a 3-cm thick bed made of coarse ash to prevalently fine, black and iridescent, poorly porphyritic, vesicular scoriaceous lapilli (85% of particles in the range 1–8 mm) with scattered pluridecimetric scoria (**Figure 8a**), amounting to a mass load of 22.25 kg/m2 . At medial sites (∼5 km distant), we measured 4.7 kg/m2 of mostly fine lapilli forming a continuous <2 cm-thick layer on the ground. Further away, the tephra fall dropped to 0.4 kg/m2 in Adrano (16 km distant), where we found a continuous, thin layer of medium-sized ash (∼85% of clasts sized between 0.25 and 1 mm; **Figure 8b**). Finally, at distal sites the fallout consisted of fine to medium ash particles, ∼99% of which were <0.25 mm in Caltanissetta (∼90 km from the vent).

Laterally away from the fallout axis, the tephra cover became rapidly discontinuous; for example, close to Torre del Filosofo (∼1 km from the NSEC; **Figure 7b**), the deposit consisted of scattered 5–10 cm-sized scoria, with larger clots up to 40 cm, while we measured 0.16 and 0.087 kg/m2 of tephra at only 2 km N and 3.3 km ∼S (**Figure 8c**), respectively, off-axis from the Adrano sample (16 km from the vent). Furthermore, within <1 km from the NSEC, we found a low number of spindle- to fusiform-shaped bombs with diameters mostly around 20–40 cm.

The two different analyzed grain-size classes of ash show comparable componentry in both samples, the ash being almost totally composed of juvenile particles (∼97% both), mainly sideromelane (87.2% in the coarser sample and 91% in the finer) with minor percentages of tachylite (9.8% and 5.8%, respectively). The few present lithics are mostly gray-dark to gray lava particles and secondarily reddish clasts, while crystals (mainly plagioclase) are negligible (0.2 and 1.6%, respectively).

The fallout deposit was estimated having a TEM of 2.25 ± 0.76 × 108 kg, which is equivalent to a mean MER of 7.81 × 103 kg/s providing a DRE magma volume of 8.33 × 104 m3 (assuming an average density value of 2700 kg/m3 ; Harris et al., 2010) (**Table 2**).

## Lava Flows

On early 26 October, lava overflowed from the crater above the NSEC cone in three spots, two on the SW crater rim and one on the E rim (**Figure 7b**). The first overflow from the SW rim fed a narrow lava tongue that spilled southward, broadening on the flat terrain at the southern base of the cone and expanding into several lobes. Later, as lava fountaining intensified the effusion rate increased, leading to the breaching of the crater rim in the place where the lava was overflowing (**Figures 6a,b**), and to a second overflow a few tens of meters further west, which fed a rather narrow lava flow that spilled down the south flank of the cone for about 100 m before uniting with the earlier lava flow.

The resulting lava flow-field to the south of the NSEC expanded into four main lobes (**Figure 7b**): one toward SW, stopping at the base of the small 1971 cinder cone (**Figure 6c**); another circled the NW base of the 2002–2003 cone (stopping short of the "Torre del Filosofo" monitoring station that fortuitously survived, only to be destroyed by the next lava fountaining episode, on 11 November 2013; **Figure 6d**) and coming then to a halt near the N base of the prehistoric cone Monte Frumento Supino; the third advanced SSE near the NE base of the 2002–2003 cone; finally, the fourth went toward SE, reaching the site of the monitoring station "Belvedere" on the rim of the Valle del Bove. The most advanced lava front, near Monte Frumento Supino, was about 1.5 km from the NSEC. The eastern lava flow was emitted at a late stage of the paroxysm, and spilled down through a deep notch formed on 27 April 2013 in the flank of the NSEC cone, largely covering a lava flow emitted at that date, and stopping at a distance of about 1.3 km from the crater.

Based on the minimum and maximum values of thicknesses, we estimate an average total volume of lava emitted on 26 October 2013 of 1.84 ± 0.44 × 106 m3 (**Table 2**), 86% of which constitutes the main southern lava flow-field.

## DISCUSSION

## How Representative Is the 25-26 October 2013 Eruption Within the Paroxysmal Activity of Etna?

In order to ascertain the representativeness of the 25-26 October eruption among the almost 200 paroxysms that occurred in 1998–2015 from the SEC and NSEC, in the following we compare physical parameters, textural features and observations on lava and tephra deposits with previous, well-studied paroxysmal episodes at Etna and particularly with the 2011–12 and 2013 sequences (**Table 1**). Our study shows that the paroxysmal phase of the 25-26 October lava fountain lasted ∼8 h, producing 1.84 ± 0.44 × 106 m3 of lava flows, and 2.25 × 0.76 × 108 kg of distal tephra as fallout deposit (corresponding to a mean MER of 7.81 × 103 kg/s). Lava flows reached a maximum length of 1.5 km from the vent, extending for 0.59 km2 at a mean effusion rate

FIGURE 7 | Mapping of the volcanic deposits erupted during the 25-26 lava fountain 2013 episode. (a) Tephra fallout on the shaded relief of Sicily; the red lines are the isomasses in kg per square meter (values in full white circles). Full white squares are the main cities in this part of Sicily (Catania, Enna e Caltanissetta). Full yellow circles indicate some localities affected by the fallout (AD, Adrano; LF, Leonforte; MT, Montedoro). The bottom right inset shows the fallout deposit drawn over the map showing all of Sicily. (b) The NSEC cone and the associated lava flows on the shaded relief of Mount Etna volcano. TDF is the location of a former alpine shelter buried under recent tephra and lava, MFS is Monte Frumento Supino.

FIGURE 8 | The tephra deposit associated with the 25-26 October 2013 distal fallout. (a) Proximal lapilli layer and 40-cm sized scoria deposited at about 1 km of distance from the NSEC; (b) continuous layer of ash deposit along the dispersal axis at 16 km distant in Adrano; (c) lateral outcrop showing a scattered ash deposit in Adrano. Photos by D. Andronico.

of ∼65 m3 /s. During the paroxysmal activity, lava jets reached a maximum height of 500 m and the eruption column is estimated to have risen to up to ∼8 km a.s.l. (i.e., ∼4.7 km above the NSEC).

#### Paroxysmal Phase Length

The 25-26 October episode produced magma jets rising on average between 300 and 400 m (and maximum height of 500 m), a height range commonly observed for lava fountains at Etna (Carbone et al., 2015). However, the 8 hduration of the paroxysmal phase (characterized by continuous fountaining) was much longer than usual. Lava fountain episodes exceeding 3 h, in fact, were only 12, 13 and 20%, respectively, of the total number of episodes during the 2000 (64 episodes; Alparone et al., 2003), 2011–12 (25 episodes; Behncke et al., 2014), and 2013 (20 episodes; De Beni et al., 2015) sequences. In particular, the paroxysmal phase duration is much longer than during the 2011–12 (from 0.5 to 4.5 h with mean duration of 1.54 h) and the 2013 episodes (from 0.2 to 8 h with mean duration of 2.6 h for 19 episodes).

#### Cone Growth

For the cone constructed during the 25 episodes of the 2011– 12 sequence from the NSEC, Behncke et al. (2014) evaluated by field surveys a growth of ∼190 m and a DRE volume of 9.5 × 106 m3 (proximal tephra; vesicularity of 50%). These values are on average equivalent to a DRE volume of 0.38 × 106 m3 per episode (**Table 1**). We also derived a DRE volume for the 2013 cone (De Beni et al., 2015) of 20.48 × 106 m3 , corresponding to a DRE volume of 0.54 × 106 m3 per episode (averaging 19 of the 21 episodes). By elaborating DEMs from different periods, for the 25-26 October activity we assume, in terms of erupted DRE volume, the average of the 2011– 13 period (44 episodes) of 0.45 × 106 m3 (**Table 1**). Such a value is in good agreement with any average episode of the 2011–12 (slightly higher) and of the 2013 sequences (lower), which was mostly characterized by large-scale fountaining episodes (e.g., Andronico et al., 2015; Poret et al., 2018a,b in revision).

#### Distal Tephra

The distal tephra fallout was characterized by a TEM of 2.25 ± 0.76 × 108 kg, which is equivalent to a mean MER of 7.81 × 103 kg/s. These values are slightly lower than the 4–5 September 2007 (3.9 × 108 kg−1.1 × 104 kg/s) and 12–13 January 2011 (1.5 × 108 kg−2.5 ± 0.7 × 104 kg/s) lava fountaining episodes (Andronico et al., 2008a, 2014b), but a little higher than those calculated for the 24 November 2006 paroxysm (1.9 × 108 kg−5 × 103 kg/s; Andronico et al., 2014a), characterized by alternating powerful Strombolian activity and pulsating fountaining. Physical parameters, therefore, would suggest that, during its 8 h of paroxysmal activity, the style and intensity of the 25-26 October 2013 eruption was that of a smallscale lava fountain (Andronico et al., 2015). Expressed in DRE magma volume they provided TEM and MER equivalent to 0.083 × 106 m3 and 2.89 m3 /s respectively (**Table 2**).

The componentry analysis shows that nearly all of the ashsized products (∼97 %) were constituted by juvenile material (of which sideromelane particles were >87 %). In agreement with the physical parameters, such high values are common during fountaining activity (e.g., Andronico et al., 2008a) as it is well suggested by the steady eruptive dynamics without any evidence of significant variations in the magma flux throughout the long paroxysmal phase, i.e., absence of significant syn-eruptive phenomena (like inner sliding and collapses within the eruptive conduit). This also permits to establish that the tephra forming the distal fallout is almost all represented by new, fresh magma.

#### Lava Flows

The 25-26 October 2013 lava flows show unusual geometric features when compared with other lavas of the 2011–2013 effusive activity (Behncke et al., 2014; De Beni et al., 2015). The morphology of the two lava fields reflects the different slopes of the NSEC cone, i.e., the smooth south flank and steep east flank, over which lava flows propagated from the respective vents up to only 1.3 and 1.5 km in around 8 h, covering a total area of 0.59 km2 (**Table 3**) with average thicknesses of 3 m. In spite of the longer paroxysmal phase, the 25-26 October lava-covered area is comparable to areas measured for the lava flows emitted during the sequences in 2011-12 (25 episodes) and 2013 (11 episodes), which were respectively between 0.29 and 1.32 km2

(mean value of 0.76 km2 ) and between 0.02 and 1.3 km2 (mean value of 0.55 km2 ) (data estimated from Behncke et al., 2014 and De Beni et al., 2015) (**Table 3**). Conversely, the lava flow length is significantly lower than the mean length measured in 2011–12 (2.94 km for 25 episodes) and 2013 (2.06 for 11 episodes). The DRE volume of lava emitted on 25-26 October 2013 is 1.47 × 106 m3 (**Table 2**), higher than the mean DRE values of the 2011–12 and 2013 sequences of 0.90 × 106 m3 (averaged for 25 episodes) and 1.04 × 106 m3 (19 episodes), respectively (**Table 4**).

The 25-26 October lava flows were emitted at a DRE effusion rate of 51.11 m3 /s (**Table 2**), a value much lower than episodes during the 2011–2013 sequences, where the mean effusion rates are estimated 161.62 m3 /s in 2011–12 (25 episodes) and 88.72 m3 /s in 2013 (19 episodes) (**Table 4**).

Comparison between the 25-26 October event and selected events or averaged values of the physical parameters/durations of the 2000, 2011–12, and 2013 sequences shows that the study event has a similar areal extent and DRE volume on the one hand, and minor lava flow length and DRE effusion rate on the other hand (compared to the 2011–13 paroxysmal sequences). Conversely, the fallout deposit erupted as distal tephra is only slightly smaller than other studied deposits, and the volume increase of the cone during this event is comparable to the average volume increase derived from the 2011-2013 sequences. However, the longer duration of the paroxysmal phase along with a significant lower DRE effusion rate during the 25-26 October event seem to indicate a relatively slow emplacement of the lava flows. We infer that their limited capability of advancing could be related to the low content of gas in the lava, possibly fed from a residual, degassed magma batch inherited from the first paroxysmal sequence in 2013, which started on 19 February and ended on 27 April. We conclude that, although eruptive style and dynamics are fully overlapping and comparable with previous events, the 25-26 October 2013 event is a small-scale paroxysmal episode as confirmed by the longer duration of the paroxysmal phase associated with a lower magma output.

#### Magma Budget Estimation From Explosive and Effusive Deposits

In general, the estimation of the magma budget provided by a volcano during its history is a good proxy for understanding its magma dynamics, how the feeding system works (e.g., Wadge et al., 2006), and how magma (degassed or undegassed) is transferred to the surface (Allard, 1997). The peculiarity of a volcano, which erupts as frequently and often displays episodic eruptions as Etna, also bears significant implications for volcanic hazard assessments. Estimating the magma budget available during paroxysmal events, in particular from the SEC or NSEC, and as long as this value is constant or at least similar within the same paroxysmal sequences, may be therefore crucial for assessing their potential impact with time on the territory. In other terms, how much is the magma erupted as tephra (distal plus proximal) and lava? Does the ratio between these two amounts change or remain constant as suggested by recent studies based on the analysis of borehole dilatometers (Bonaccorso et al., 2014)? While the lava flows may be evaluated both using standard methods (mainly direct mapping and subtraction of DEMs extracted from remote sensing data; Coltelli et al., 2007; Behncke et al., 2014) and new techniques such as infrared data from geostationary sensors (Ganci et al., 2012) and near real-time production of DEMs (De Beni et al. under review), calculating the portion of magma fragmented "explosively" is more difficult. The most common models used for the fallout deposits, in fact, do not sufficiently take into account the

We reported: the total lava covered area and total lava length for the two sequences, and the average values per episode.

TABLE 4 | Volcanological parameters of the NSEC cone during the 2011–12 and 2013 lava fountaining sequences, and the single 25-26 October 2013 paroxysm.

We reported: the total lava volume (bulk and dense rock equivalent-DRE) and the correspondent average values per episode.

coarse-grained tephra which does not enter into convection within the eruption column. Conversely, during low- to midintensity explosive eruptions, like lava fountaining and powerful Strombolian events at Etna, which typically produce low, <10 km-high eruption columns a.s.l., the contribution of very coarse clasts that coevally fall back into and closely around the active vent of a rapidly growing scoria cone is significant.

Considering the contribution of proximal tephra and lava flows emitted during the 2011–12 and 2013 sequences, we quantified the average DRE magma volume erupted for episode. For the 25 episodes of the 2011–12 sequence, we obtain a total DRE volume of 31.9 × 106 m3 , i.e., 1.28 × 106 m3 per episode, a value similar to the mean DRE volume of 1.47 × 106 m3 estimated averaging 19 of 25 episodes by Ganci et al. (2012) using a geostationary sensor (SEVIRI). Considering 138,600 s for the combined duration of the paroxysmal activity of all eruptive episodes (Behncke et al., 2014), we obtain a total mean DRE magma eruption rate (MGER) of ∼230.16 m3 /s, i.e., 68.54 × 106 m3 /s and 161.62 × 106 m3 /s the contributions from the cone and the lava, respectively, without including the magma volume erupted as distal tephra. By the same calculation averaging 19 of the 21 episodes in 2013, we derived a DRE volume of lava plus cone of 3.00 × 107 m3 , and a DRE MGER of 168.47 m3 /s.

Considering the 25-26 October activity, in terms of DRE erupted volumes we estimated 15.63 m3 /s for the cone and 51.11 m3 /s for the lava flows, respectively (**Table 2**). Such estimations indicate that the cone increase was on average a fourth than during the 2011-13 episodes (**Table 1**), while the lava output value is almost a third of the mean value obtained for the 2011-12 episodes from Behncke et al. (2014), and ∼half the value obtained for 19 of 21 episodes in 2013 by De Beni et al. (2015) (**Table 4**), suggesting that the lava effusion from the studied episode occurred at a lower rate compared to the averaged values of the 2011-13 paroxysmal activity. In addition, we obtained 0.083 × 106 m3 of magma DRE volume erupted as distal tephra, which accounts for a further DRE MGER of 2.89 m3 /s (**Table 2**); for the 2011-12 lava fountains, this parameter was only approximately estimated by Behncke et al. (2014) and not calculated by Ganci et al. (2012). Summing together proximal tephra, lava flows, and distal tephra, for the 25-26 October paroxysmal activity we have a total mean DRE MGER of 69.63 m3 /s (**Table 5**). This value is significantly lower than 230.16 and 146.28 m3 /s derived by Behncke et al. (2014) and De Beni et al. (2015) for the 2011-12 and 2013 episodic activity (lava and cone contribution without distal tephra).

## Magma Dynamics During Paroxysmal Activity at Etna

In the literature, the eruptive mechanism invoked for describing sustained lava fountains from the summit craters of Etna (e.g., Allard et al., 2005; Andronico and Corsaro, 2011) is the separate ascent of a bubble foam layer related to a gas-melt separation. This is the collapsing foam model, which implies the accumulation of volatiles at the top of the plumbing system, the formation of a foam layer (e.g., Jaupart and Vergniolle, 1988, 1989; Vergniolle and Jaupart, 1990), and its dramatic collapse emptying the magma reservoir. Allard et al. (2005) postulated that the 14 June 2000 lava fountain at Etna was supplied by a shallow magma body on top of which "repeated growth and collapse of a bubble foam may have triggered periodic lava fountaining." They also measured distinct chemical compositions of the gas driving the magma jets and the lava flows, suggesting someway the separation from the parcel already stored in the shallow reservoir resulting into the emission of fragmented magma and degassed magma by distinct vents. Polacci et al. (2006) inferred that the foam collapse process may be followed and accompanied by syneruptive, although minor, volatile degassing. Consequently, the occurrence of lava fountaining suggests in general the availability of a pressurized volume of magma (composed of solid, liquid and abundant volatiles) in the shallow reservoir.

The 25-26 October lava fountain had a longer duration but lower mean MGER than the respective mean values calculated for the 2011–12 and 2013 sequences. Volcanological data, camera recordings and direct observations show that the intensity of the paroxysmal activity was relatively low. Noteworthily, most of the magma emitted was discharged as lava flows and proximal tephra. The distal tephra, in fact, represents only around 4% of the total erupted DRE volume, while the total volume of pyroclasts (proximal plus distal tephra) is lower (around 36%) than the magma volume erupted as lava flow (**Table 1**). This means that approximately three-quarters (i.e., 73%) of the magma was emitted as lava flows, i.e., it was already degassed and/or crystallized for the most part. The remaining >25 % of magma was erupted explosively, and mostly (>80 %) fragmented and

We reported: the MGER DRE total cone volume, the MGER DRE total lava volume, the sum of the previous values and the MGER DRE tephra (see Table 2). DRE, dense rock equivalent; MGER, magma eruption rate.

deposited as coarse tephra clasts above the NSEC cone. **Figure 9** summarizes the DRE magma volume (%) distribution as derived from our work in terms of lava, distal tephra and proximal tephra deposits.

In other terms, on 26 October 2013 the fraction of magma erupted as tephra and transported in convection above the NSEC was rather small compared to the magma emitted as lava flows and cone-building tephra. The importance of the extruded lava mass and of the effusive process in general is also suggested by the detachment of blocks from the NSEC summit and carriage downward. Noteworthy, the presence of a dark eruption column above the volcano which then becomes progressively lightercolored for tens of km away from Etna could be related to high amounts of water vapor (**Figures 3c,d**), but foremost points to a relatively low discharge of tephra, the fine-grained part of which diluted in the atmosphere.

Although the eruption column that typically accompanies the common lava fountains at Etna is the most striking eruptive phenomenon, and possibly the one, which distinguishes this activity from the "fire-fountaining" activity at Hawaiian volcanoes, its contribution in terms of magma volume, at least for the studied event (but probably also for other small-scale lava fountain episodes), was only a small fraction of the total budget of erupted magma.

In conclusion, we infer that before 25 October, as already proposed for past episodic eruptions (Andronico and Corsaro, 2011; Behncke et al., 2014), a cooling magma batch was resident below the NSEC. This batch was extruded for the most part as lava flows, possibly as result of degassing during the 6 months-long residence in the shallow plumbing system since the previous paroxysmal episode on 27 April. However, it opened the plumbing system to the rise of more gas-rich magma, which was erupted during the following five lava fountains episodes (e.g., the exceptionally powerful 23 November 2013 episode; Bonaccorso et al., 2014; Poret et al., 2018a) up to 2 December 2013, plus two powerful Strombolian events in the following month. Similarly, for the 64 paroxysmal episodes that occurred in January-June 2000, Andronico and Corsaro (2011) showed that the ascent of new, more primitive and volatile-rich magma triggered the start of this long episodic eruption, progressively intruding into the SEC reservoir and mixing with the resident and more evolved magma, finally determining the composition of erupted volcanics. Noteworthy, in the case of the 2000 sequence, the more primitive magma was emitted only during the last events.

The latter hypothesis is also in agreement with two recent studies from Viccaro et al. (2014) and Giuffrida et al. (2018) about the triggering mechanisms of paroxysmal episodes at Etna. Viccaro et al. (2014) infer that textures and compositional zoning in selected plagioclase crystals of lavas can explain the relations between the duration of the paroxysms and dynamics of pre-eruptive magmatic processes at depth. They proposed a recharge by more mafic, gas-rich magma in case of prolonged, ramp-shaped tremor increase, and conversely a prominent role of gas injections into the residing system when recording a sudden increase of the seismic amplitude before the climax of the eruption. The 25-26 October 2013 episode,

indeed, shows intermediate features between the previous endmembers, because the tremor increase accompanying the initial Strombolian phase was neither long nor short, while the paroxysmal phase was clearly longer than usual. At the same time, by studying lithium gradients in plagioclase crystals, Giuffrida et al. (2018) suggest that, during the 2011–2013 paroxysmal episodes, the fast recharge of gas-rich magma would have caused the re-activation of the magmatic system after quiet eruptive periods, thus enhancing powerful eruptions. We therefore cannot exclude the presence of two distinct mechanisms which concurred to open a short sequence of 8 paroxysmal events after 6 months of quiescence. We finally suggest that the rise of a new mafic, gas-rich magma batch triggered the 25-26 October episode and provoked the extrusion of degassed resident magma, accompanied by the resumption of the Strombolian activity and sustaining a long-lasting paroxysmal phase.

## CONCLUDING REMARKS

The ascent of magma during paroxysmal eruptions and its partition to produce disproportionate values of lava/tephra ratios as well as a wide range of MERs have been never properly investigated and debated at Etna. In spite of various research on single paroxysmal episodes in recent years, to date only this study quantifies, on the base of field data, the DRE magma volumes erupted as proximal tephra, distal tephra, and lava flows. The 25-26 October 2013 lava fountain deposited a proximal tephra deposit of 0.45 × 106 m3 on top of the NSEC cone (causing ∼30 m growth of the cone), dispersed a distal tephra fallout deposit of ∼0.083 × 106 m3 , and emitted ∼1.47 × 106 m3 of lava flows. Considering 8 h of paroxysmal

activity, the mean rate at which all the erupted magma was supplied is 69.63 m3 /s. The estimation of the total mass and the relative percentages of lava and tephra may address the evaluation of the eruption dynamics during paroxysmal episodes at Etna, and gain insights into the wide spectrum of eruption intensity (MER) and style (powerful Strombolian to large-scale lava fountain events) observed in the last years at the summit craters. Most importantly, our study highlights that, at least for ordinary paroxysmal episodes, the magma budget distribution is strongly balanced in favor of lava flows and proximal tephra (accumulating around the NSEC cone and enhancing its fast growth), while the distal tephra fallout seems to have a rather modest importance in quantitative terms (i.e., only <5% of the total erupted magma). A similar distribution of magma (predominant output of lava flows/coarse tephra around the cone than distal tephra) is coherent with the eruptive behavior observed during the paroxysmal phase of most of the lava fountain episodes at Etna, i.e., the start of lava fountaining and formation of the eruption column, which are accompanied by high-rate lava effusion.

In the future, further and more reliable estimations of the different contributions of magma discharged as lava and tephra during paroxysmal events are desirable. The application of new measurement techniques (like drones and groundbased photogrammetry) is likely to allow swift calculation of detailed and reliable volume values, not only of lava flows but also of the proximal tephra which accumulates on top of the cone; for example, by subtraction of DEMs focused on the NSEC cone just before and after an eruption. Moreover, coeval and punctual geochemical measurements of the gas plumes are recommended to correlate the total SO2 flux and in situ measurements of gas phases (e.g., SO2, HCL, HF, and their relative ratios) with the eruptive activity, thus allowing to validate the total volume of erupted magma and possibly, its partitioning in terms of lava and tephra.

Understanding the relationships between effusive and explosive deposits (i.e., degassed vs. gas-rich magma) erupted during single paroxysms, finally, may be greatly useful to predict the evolution of such episodic eruptions and their likely conclusion, and thus also to improve the assessment of the

#### REFERENCES

potential volcanic hazard from future tephra fallout over the densely populated slopes of Etna.

#### AUTHOR CONTRIBUTIONS

DA conceived the research and drafted the article. DA and BB made field observations during the whole duration of the lava fountain. DA carried out survey to map the tephra fallout from the proximal area to the volcanic slope, and collected representative samples. BB and ED carried out field campaigns to map the lava flow field. ED produced the cartographic maps for the lava and tephra deposits, and elaborated the obtained maps to calculate the volume of the cone and of the lava deposits. AC analyzed the products of the tephra fallout in the laboratory. SS estimated the volume of the tephra fallout. ML surveyed the distal fallout deposit and collected the relative ash samples. MLC collaborated during the laboratory analysis of the tephra samples. All authors discussed the results and reviewed the manuscript, strictly collaborating for the interpretation of the results obtained. BB also gave a substantial contribution to the draft of the manuscript.

#### FUNDING

Osservatorio Etneo—Sezione di Catania Istituto Nazionale di Geofisica e Vulcanologia.

#### ACKNOWLEDGMENTS

We are strongly indebted with the INGV staff, which maintained the video-camera network during the period of the eruption, and in particular our colleagues who were in the field to make observations from different viewpoints at different time during its long paroxysmal phase. Finally, we are really indebted to reviewers, Prof. Marco Viccaro and Dr. Luis E. Lara, whose criticism and observations gave more clarity to the paper. We also thank the Editor, Prof. Clive Oppenheimer, for his continuous support and the Chief Editor, Prof. Valerio Acocella, for the final editing to the revised paper that has given greater completeness and self-sustainability to our work.

Fierstein, J., and Nathenson, M. (1992).Another look at the calculation of fallout tephra volumes. Bull. Volcanol. 54, 156–167. doi: 10.1007/BF00278005

Ganci, G., Harris, A. J. L., Del Negro, C., Guehenneux, Y., Cappello, A., Labazuy, P., et al. (2012). A year of lava fountaining at Etna: volumes from SEVIRI. Geophys. Res. Lett. 39:L06305. doi: 10.1029/2012GL051026

system from micro-analytical studies of crystals. Earth Sci. Rev. 171, 289–322. doi: 10.1016/j.earscirev.2017.06.003

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Andronico, Behncke, De Beni, Cristaldi, Scollo, Lopez and Lo Castro. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Assessing Impact to Infrastructures Due to Tephra Fallout From Öræfajökull Volcano (Iceland) by Using a Scenario-Based Approach and a Numerical Model

Sara Barsotti1 \*, Dario Ingi Di Rienzo1,2, Thorvaldur Thordarson2 , Bogi Brynjar Björnsson1 and Sigrún Karlsdóttir1

1 Icelandic Meteorological Office, Reykjavík, Iceland, 2 Faculty of Earth Sciences, University of Iceland, Reykjavík, Iceland

#### Edited by:

Sonia Calvari, Istituto Nazionale di Geofisica e Vulcanologia (INGV), Italy

#### Reviewed by:

Marco Viccaro, Università degli Studi di Catania, Italy Eisuke Fujita, National Research Institute for Earth Science and Disaster Prevention, Japan

> \*Correspondence: Sara Barsotti<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Volcanology, a section of the journal Frontiers in Earth Science

Received: 02 June 2018 Accepted: 19 October 2018 Published: 13 November 2018

#### Citation:

Barsotti S, Di Rienzo DI, Thordarson T, Björnsson BB and Karlsdóttir S (2018) Assessing Impact to Infrastructures Due to Tephra Fallout From Öræfajökull Volcano (Iceland) by Using a Scenario-Based Approach and a Numerical Model. Front. Earth Sci. 6:196. doi: 10.3389/feart.2018.00196 Mt. Öræfajökull is one of the most dangerous volcanoes in Iceland with potential for a VEI6 eruption and the generation of many severe associated hazards. It is not a frequently erupting volcano with two eruptions in the last 1100 years, in 1362 and 1727–28. During the 1362 eruption 10 km3 of freshly fallen tephra was emitted, the eruption plume reached the stratosphere and was dispersed offshore toward mainland Europe. In this study we investigate the possible impact due to tephra fallout to critical infrastructures in Iceland namely – roads, airports, electrical power-lines – in case of a new eruption at Öræfajökull of similar intensity as in 1362. The analysis is done by running several times the VOL-CALPUFF dispersal model to simulate the dispersal of ash in the atmosphere and its deposition on the ground. The resulting maps show the probability of exceeding critical thickness of the tephra fall. Critical infrastructures have been added to the analysis to get a quantitative assessment of the potential impact. The results indicate that in case of an event similar to the 1362 eruption, the tephra fallout could be expected over most of the country, with higher likelihood on the eastern side. The tephra fallout is likely to have a severe impact in the proximity of the volcano, generating a deposit with a load of up to 1000 kg/m2 . The likelihood of failure for more than 160 km of the electrical power-line and for critical driving conditions on about 900 km of the main ring road is between 50 and 100%. The probability that the tephra fall will affect three of the main domestic airports is higher than 50%. An eruption of this magnitude is likely to affect commuting and communication between the greater Reykjavík area, where the government resides, and the rest of the country. Our analysis also reveals the limitations of current knowledge and understanding of the Öræfajökull volcano and highlights the need for further studies on past activity to better characterize its future behavior.

Keywords: Öræfajökull volcano, Iceland, tephra fallout, 1362 eruption, hazard assessment, numerical model, critical infrastructure

## INTRODUCTION

feart-06-00196 November 13, 2018 Time: 12:17 # 2

There are about 30 active volcanic systems in Iceland and about half of those have featured tephra-producing eruptions (Thordarson and Höskuldsson, 2008). Volcanic eruptions are common in Iceland and have a recurrence interval of 2–5 years (Thordarson and Larsen, 2007). Basaltic eruptions are the most common volcanic events, and among them explosive within-glacier eruptions are most frequent because the most active central volcanoes are capped by glaciers (e.g., Katla, Grímsvötn, Bárðarbunga). Less frequent are explosive eruptions featuring more evolved magmas, such as dacite and rhyolite, that typify central volcanoes such as Öræfajökull and Hekla (Larsen and Eiríksson, 2008a,b). Highly active volcanic systems, as Hekla, Katla, Bárðarbunga and Grímsvötn, have explosive eruptions rates of 82, 97, 90, and 95%, respectively (CIV, 2017).  Tephra dispersal and fallout is by far the most widespread hazard affecting local as well as distal regions. Ash clouds and tephra fallout can cause severe health issues (Baxter, 1990; Horwell and Baxter, 2006), affect important infrastructure like as electrical supply systems (Wilson et al., 2012), the national and international transportation network (Guffanti et al., 2009; Wilson et al., 2012), sensitive buildings (Spence et al., 2005), human health and life stock, vegetation and eco-system (Thorarinsson, 1979; Wilson et al., 2012; Ágústsdóttir, 2015).

It was during the infamous eruption at Eyjafjallajökull in 2010 when a persistent northwesterly winds carried the ash-rich plume toward Europe for more than a month (Baerbel et al., 2012; Gudmundsson et al., 2012). On that occasion a prolonged closure of the European airspace resulted in severe economic impact estimated to be € 1.3 billion in the first week of the eruption (Bolic´ and Sivcev, 2011 ˇ ). The southernmost tip of Iceland experienced heavy tephra fallout, which impacted the local residence in various ways. Situation of low visibility happened often during the eruption as well as in its aftermath because of resuspension of ash (Petersen, 2010; Karlsdóttir et al., 2012). Air quality was often poor and affected the health of population living closest to the volcano (Carlsen et al., 2012). Relocation of life stock became essential due to heavy tephra fallout (Karlsdóttir et al., 2012; Thorvaldsdóttir and Sigbjörnsson, 2015). Specific investigations were done during the Eyjafjallajökull eruption to assess the effect of ash contamination on electrical power plants (Rarik, 2010). Eruptions at Hekla volcano have also been investigated to assess their effect on the environment and ecosystem (e.g., Frogner et al., 2006). Heavy tephra and lapilli fallout during the eruption at Heimaey (Vestmannaeyjar Island) in 1973 destroyed and damaged several houses (Williams, 1983; Spence et al., 2005; Gudmundsson et al., 2008), some of those were restored after extensive cleaning effort to remove the tephra fall deposit.

In the period 2013–2016 the Icelandic Government, together with the International Civil Aviation Organization (ICAO), supported several projects aimed at assessing in a quantitative manner the long-term volcanic hazard in Iceland, to be used for more in-depth risk analysis. A specific component of these projects was the investigation of the impact of tephra fallout in Iceland for both medium-size, more frequent, and large, less frequent, explosive eruptions. Öræfajökull volcano was selected as the low probability but high-impact scenario.

Here, we focus on the AD 1362 Öræfajökull event as the most extreme scenario for a regional tephra fallout hazard and a preliminary risk assessment. The sparse eruption records for Öræfajökull introduces uncertainties in the volcanological scenario considered and we address this issue by using a numerical model to investigate the sensitivity of model results to variations in the volcanological input parameters.

The results of this study are presented through probabilistic hazard maps. These type of maps have become a commonly practiced representation of volcanic hazards and helps with the visualization of the footprint of the volcanic phenomena that may impact the surroundings of a volcano (Haynes et al., 2007; Nave et al., 2012; Calder et al., 2015; Thompson et al., 2015). Different types of hazard maps exist in literature. They can be produced on the basis of geological data, historical records and/or numerical model results (Calder et al., 2015; Loughlin et al., 2017). They can refer to a past eruption, to a specific hypothetical eruptive scenario or to a distribution of scenarios. If based on numerical results they can show the results from a single specific simulation (deterministic map) or from a multitude of scenarios. In the latter case the maps are often representing the impact of a specific hazard as a spatial probability and we refer to them as "probabilistic hazard maps." Hazard maps are often used by volcano monitoring institutions to inform their stakeholders (e.g., decision makers institutions, general public, emergency managers, land-plan managers) about areas prone to be affected by specific hazards in case of an eruption. On a map it is easy to visualize extent of borders plus location of sensitive infrastructures, roads, towns and villages and, therefore, put the hazards into a spatial context that can be perceived more effectively by the users. Volcanic hazard maps have been produced for several volcanoes using numerical models and are applied for long-term hazard and risk assessment at particular volcanoes. For example, hazard maps have been produced for pyroclastic density currents at Mt. Vesuvius and Napolitean area (Esposti Ongaro et al., 2002, 2012; Sandri et al., 2018); for lava flows at Etna (Favalli et al., 2005; Tarquini and Favalli, 2013), Nyiragongo volcano (Favalli et al., 2009) and Lanzarote (Felpeto et al., 2001); for volcanogenic floods at Öræfajökull volcano (Pagneux et al., 2015). Probabilistic hazard maps for tephra fallout have been produced for Mt. Etna (Scollo et al., 2013), Campi Flegrei (Costa et al., 2009), Tarawera volcano (Bonadonna et al., 2005), Indonesian volcanoes (Jenkins et al., 2012), Santorini volcano (Jenkins et al., 2015). Most recently probabilistic maps for hazard due to ejection of ballistic have been produced for Mt. Chihshin in North Taiwan (Nurmawati and Konstantinou, 2018).

Here, we have produced probabilistic hazard maps to investigate the potential impact at a national level of a VEI6

FIGURE 1 | Map showing the locations of Öraefajökull volcano and the main sites mentioned in the text. On top of it the reconstructed deposit for the 1362 AD eruption is displayed.

(Newhall and Self, 1982) eruption at Öræfajökull. The impact of tephra fallout on infrastructures is calculated via a numerical model that simulates the atmospheric dispersal of tephra and its deposition. GIS-referenced layers with information on the powerline network, roads, airports are included to assess potential disruptions to services and commuters. Possibly, these results will support the national authorities in planning and designing mitigation actions that are necessary to reduce the risks posed by a future eruption at Öræfajökull. A more complete hazard assessment should also include the investigation of smaller eruptions (VEI4) as they appear to be most frequent events at Öræfajökull as reported in Section "Geological Background." The choice to look first at a VEI6 eruption was dictated by the need to quantify the potential damage the worst case scenario would cause to the society nowadays. The indications provided by the Icelandic Civil Protection addressed this research and supported the worst case scenario as the reference scenario for Öræfajökull.

#### Geological Background

Öræfajökull is an ice-capped stratovolcano located in South-East Iceland on the southern margin of Vatnajökull glacier (Thordarson and Larsen, 2007; Gudmundsson et al., 2008; Sharma et al., 2008; Larsen et al., 2015). It is about 20 km in diameter with a 3 by 4 km ice-filled caldera which rises to a summit of 2110 m a.s.l. (**Figure 1**). The volcano is part of the intraplate Öræfajökull Volcanic Belt, situated to the east of the current plate margins and possibly represents an embryonic rift (e.g., Thordarson and Larsen, 2007; Thordarson and Höskuldsson, 2008).

The Öræfajökull central volcano has featured two explosive eruptions in historical times (e.g., Thorarinsson, 1958). The most recent was a small icelandite eruption of VEI4 in 1727–1728 CE (e.g., Larsen et al., 1999; Larsen et al., 2015). This was preceded by a much larger rhyolitic Plinian (VEI6) eruption in 1362 CE. Studies on the tephra stratigraphy in soils around the volcano have revealed five prehistoric silicic explosive eruptions at Öræfajökull and all are assumed to be smaller in magnitude and intensity than the 1362 CE event (Gudmundsson, 1998; Larsen et al., 2015).

In Iceland VEI6 events are infrequent and only three such events are known during the Holocene: the before mentioned 1362 CE event at Öræfajökull along with the two largest Holocene silicic explosive eruptions in Iceland, the 3 ka H3 and 4.2 ka H4 events at the Hekla volcanic system (e.g., Larsen and Thorarinsson, 1977; Stevenson et al., 2015).

The 1362 CE Öræfajökull eruption is the largest rhyolitic eruption in Iceland since settlement in the 9th Century, with an estimated volume of freshly fallen tephra of 10 or ∼6 km3 of compacted tephra and about 2 km3 when calculated as dense rock (Thorarinsson, 1958). Sharma et al. (2008) obtained a smaller volume of 2.3 km3 assumed to equal 1.2 km3 calculated as dense rock. Early stage pyroclastic density currents and intercalated jökulhlaups, along with the subsequent tephra fall, inundated the then prosperous farming district "Litla Hérað" causing fatalities (e.g., Thorarinsson, 1958; Jónsson and Valdimar, 2007; Thordarson and Höskuldsson, 2007). The reconstructed tephra dispersal is shown in **Figure 1**, where isopachs (i.e., lines of

The mass estimated to match the reconstructed deposit with the modeled results is also shown.

equal thickness) for the 1362 eruptions are shown as black lines, and the dashed lines indicate the inferred dispersal over the sea. About half of the country received >1 mm of ash as a consequence of the tephra fall from this eruption. Close to the volcano up to 20 cm of ash has been accumulated over an area of 1000 km2 , peaks in the deposit thickness are found in Grófarlækur (40 cm) at about 10 km from the summit volcano (Thorarinsson, 1958) and between Hnappavellir and Fagurhólsmýri where the thickness reached 2 m (Jónsson and Valdimar, 2007; Sharma et al., 2008) (**Figure 1**). Tephra from Öræfajökull has been identified in Western Europe (e.g., Pilcher et al., 2005) and in Greenland ice-cores (e.g., Palais et al., 1991).

The magma erupted in 1362 is a medium-K alkali rhyolite that is extremely homogeneous in composition. It is thought to be extracted from the topmost portion of a compositionally stratified magma storage zone, although there is no geophysical indication of any upper crustal (<15 km depth) magma storage zone beneath Öræfajökull at present. In contrast to the East Volcanic Zone and other rift zone rhyolites, which demonstrably are generated by partial melting of hydrated metabasaltic crust (e.g., Sigmarsson et al., 1992), the Öræfajökull and other silicic magmas generated by intraplate volcanoes are inferred to have evolved from a more mafic parent via fractional crystallization

(e.g., Selbekk and Tronnes, 2007; Martin and Sigmarsson, 2007; Sigmarsson et al., 2008).

## Recent Unrest

From September 2017 Öræfajökull volcano has showed clear indications of unrest. The reinvigoration of the geothermal activity beneath the volcano was corroborated by elevated seismicity, gas release and the formation of a cauldron in the middle of the caldera. Consequently, the Icelandic Meteorological Office raised the color code for Öræfajökull to yellow1 in November 2017. The aviation color code was turned back to green on 4th May, 2018 when the main monitoring parameters indicated a stable situation with no immediate hazard to the aviation. However, given the persistent potential for local hazards, the Icelandic Civil Protection decided to maintain the level of Uncertainty for Öræfajökull (IMO, 2018a).

In the period November 2017 to May 2018, a priority has been placed on improving the real-time monitoring around Öræfajökull. At the same time, hazard, and risk assessments performed for glacial outbursts originating from this volcano (Pagneux et al., 2015) were provided to the Department of Civil Protection and Emergency Management of the National Commissioner of the Icelandic Police to finalize the evacuation plans for the Öræfi district (NCIP, 2017). The preliminary study on tephra fallout and its potential impact on infrastructures at a national level, here presented, was also finalized.

At the time of writing the seismicity is elevated and higher than in November 2017. The deformation data (from both cGPS and InSAR analysis) reveal clear indication of an on going inflation process likely due to injection of new magma at depth (IMO, 2018b).

## MATERIALS AND METHODS

The methodology used for the hazard assessment in this study is shown in **Figure 2**. It consists of four main steps: (1) identification of the scenario of interest (this can be defined on the basis of an Event Tree outcome or from literature); (2) initialization of VOL-CALPUFF model (Barsotti et al., 2008) by using a synthetic scenario (selected eruption source parameters to be used as a model input) as well as a range of meteorological scenarios); (3) execution of several numerical runs by using different starting times; (4) statistical processing of the results from multiple runs. The obtained probabilities, as visualized in the final maps, are those called "conditional probabilities," i.e., conditioned to the occurrence of that specific eruptive scenario.

#### Scenario Definition

Available publications (Thorarinsson, 1958; Jónsson and Valdimar, 2007; Thordarson and Höskuldsson, 2007; Sharma et al., 2008), provide the general framework for the eruption source parameters used in the dispersal simulations. Calculating the total mass emitted during the Plinian phase of the eruption is a first order estimate because large portionof tephra fell onto the sea (see **Figure 1**). Using the values given in Thorarinsson (1958) the total mass is estimated to be 4.8 × 1012 kg (**Table 1**).

A GIS referenced reconstruction of the original map by Thorarinsson (1958) is used to recalculate the erupted tephra mass by (1) assuming a constant thickness between the different isopachs and (2) by interpolating between two successive isopachs assuming a linear trend. These two estimates give an erupted tephra mass of 4.29 × 1012 and 6.96 × 1012 kg, where the latter is about factor of 1.5 larger than that obtained from Thorarinsson (1958) data (see **Table 2**).

TABLE 3 | Tephra fallout conditions investigated to cause insulator flashover and pole and line damages for wet and dry ash (as reported in Wilson et al., 2012).

1http://en.vedur.is/about-imo/news/a-new-ice-cauldron-in-oraefajokull-volcano

In order to run a dispersal model several input parameters need to be defined. VOL-CALPUFF model requires initial values of parameters like the vertical exit velocity of the volcanic mixture, the vent radius, the total grain size distribution (TGSD). As we currently have no information on the values for these variables, they are obtained via iteration by running VOL-CALPUFF model with a range of input values for the above mentioned variables. The output values selected are those that produced the best fit between the modeled and the observed deposit distribution as well as the eruption column heights. In an attempt to reproduce the scenario reported in **Table 1** and to match the isopachs as depicted in Thorarinsson (1958), the VOL-CALPUFF dispersal model was run multiple times using a vent radius from 150 to 300 m; a gas mass fraction from 1 to 5% and three different TGSDs (one peaked at phi = −3, one peaked at phi = 0 and a bi-modal distribution with two peaks at phi = −3 and phi = 2, respectively). By a comparison of the model results and the original isopach map and by constraining the top plume height between 24 and 34 km, the best fit was obtained by using the following input parameters:

A forward run of the plume model with these input parameters produced the following values for plume height and tephra emitted mass (to compare with values reported in **Table 1**):

This synthetic scenarios is in a reasonable agreement with anticipated duration of the Plinian phase and plume height. However, there is a larger discrepancy in terms of mass. The new numerical simulation results suggest that in order to match the original isopachs, the DRE volume is 5.5 km3 . This volume is 2.75 times larger than that provided by Thorarinsson, and 1.9 times larger than the value calculated with the GIS interpolation (**Table 2**). Considering the huge uncertainites affecting the meteorological conditions during the eruption, the real extension of the deposit and the few observational data available, we considered valuable in this study to use those input parameters obtained through the matching procedure to perform the probabilistic assessment.

FIGURE 4 | Model results are overlayed with the deposit isopach for the 1362 eruption. The simulation has been run by using a wind field characterized by a westerly wind to match the prevailing direction of the deposit.

The tephra dispersal from the synthetic scenario is compared with the isopachs reconstructed for the real event in **Figure 4**. The gray filled contours are the model results whereas the solid black lines are the deposit isopachs as reconstructed by Thorarinsson (1958). A westerly wind has been selected to run the dispersal and match the general feature of the deposit pattern. For the same simulation (i.e., by assuming the same input condition) the model reproduced a super-buoyant plume with height of 30.6 km above the vent.

#### Numerical Approach The Numerical Model

Plume ascent is described solving plume theory equations (Bursik, 2001) to compute column height as a function of volcanological source input data and wind field parameters. The latter are relevant for simulating weak plumes that are strongly affected by wind shearing. During plume ascent the heaviest particles fall from the column and a lighter mixture continues its upward motion, entraining air up to a null-vertical velocity altitude where only lateral dispersion takes place. The plume initially decelerates due to higher density compared with surrounding atmosphere, but due to heating of entrained air (mixed by turbulent motions) the mixture can eventually become lighter than air. Buoyancy effect can cause the mixture to accelerate upward until an equilibrium is achieved.

The dispersal code VOL-CALPUFF originates from the CALPUFF model and it has been modified to reproduce some processes specific of a volcanic eruption, e.g., plume rise phase and the dispersal of a distribution of solid particles. It is a hybrid model in which the plume rise phase is described with a Eulerian approach, whereas the ash cloud transport is solved in a Lagrangian framework. Along the plume and at its top the material is released as a series of diffusing packets (puffs) containing an initially assigned amount of particulate matter which varies during the transport due to gravitational fallout. Since its development the VOL-CALPUFF model has been applied mostly at Mt. Etna to reconstruct past explosive events (Barsotti and Neri, 2008; Barsotti et al., 2008), as an ash dispersal forecasting tool (Barsotti et al., 2008) and to estimate potential hazards posed by volcanic ash to human health and ground infrastructures (Barsotti et al., 2010). In the past years VOL-CALPUFF has been also applied to other active volcanoes to

produce forecasting maps of ash dispersal during eruptive crises at Redoubt Volcano (Alaska) in 2009, Eyjafjallajökull (Iceland) and Mount Merapi (Indonesia) in 2010 (Barsotti et al., 2011; Spinetti et al., 2013) and Grímsvötn (Iceland) in 2011.

#### Monte Carlo Simulations

Numerical models can be used to investigate the behavior of a specific process process (Kavanagh et al., 2018) as for example dispersal of ash (Bonadonna et al., 2011; Folch, 2012), lava flow invasion (Del Negro et al., 2005; Favalli et al., 2005, 2009), maximum distances of pyroclastic flow (Esposti Ongaro et al., 2012; Dufek, 2016). Each simulation needs specific input conditions to characterize the volcanological scenarios to be investigated. For volcanic ash dispersal simulation the eruptive source parameters (ESP) as plume height, particle size distribution, mass flow rate need to be quantified (Mastin et al., 2009). As we do not know in advance about the next eruptive conditions, like the weather and the ESP, a way to treat this uncertainty is to reflect this into a probabilistic analysis. Looking into a range of eruptive scenarios it is possible to get a statistics that would investigate and reflect the uncertainty in the assumption made analyzing a single scenario. By running the model with several sets of starting conditions enables us to estimate the probability that a specific area will be affected by a certain type of hazard. A method widely used to achieve such a result is called Monte Carlo approach and is based exactly on the assumption that a model could be executed multiple times for as many initial conditions as required for producing an ensamble modeling (Sparks and Aspinall, 2013). The application of this approach for a probabilistic analysis of a simulated volcanic process can be found in Cioni et al. (2003), Hurst and Smith (2004), Bonadonna et al. (2005), Macedonio et al. (2008), Costa et al. (2009), Barsotti et al. (2010), Scaini et al. (2012, 2014), Scollo et al. (2013), Biass et al. (2014), Bonasia et al. (2014).

A Monte Carlo simulation is performed for this study by running several times the dispersal model VOL-CALPUFF for a fixed volcanological scenario (the 1362 AD reference eruption) and by using several years of meteorological data. Each simulation is performed using the same input data given in Section "Scenario Definition." The simulations are initiated with different starting date and times of the day, over a period of 10 years. This is done to ensure the randomness of the analysis and to avoid bias in the results due to the daily variations of the atmospheric parameters as we only investigate the effect of the wind field statistics on the tephra dispersal. The potential variability associated with the uncertainty in the volcanological

Thorarinsson (1958) are reported as a reference with black lines.

scenario is not considered. A total of 500 simulations were performed in order to obtain a convergence for the modeled tephra dispersal distribution.

Forecast data produced by the Europen Centre of Mediumrange Weather Forecast is used to run the VOL-CALPUFF model. The meteorological data for the probabilistic maps have been downloaded from the ERA-INTERIM archive and cover a period of 10 years, from 1980 to 1991. This dataset is produced by re-analyzing the forecast, i.e., by assimilating the observational data into the original weather forecast (Dee et al., 2011). This means that this dataset provides a quite complete and verified description of the 3D atmospheric fields over the period of interest. On the other side the horizontal resolution of the meteorological data is of 0.7◦ (i.e., about 35 and 77 km in the longitude and latitude respectively), making the spatial resolution of this dataset a bit coarse for the domain considered in this project. The statistics for the investigated area around Öræfajökull produced by using the ECMWF data over this time period is shown in **Figure 5**.

In **Figure 5** each windrose shows the direction of provenance for the wind field, identified by the sector, and the wind velocity, identified by the color. Each windrose is produced for a specific altitude expressed in pressure from 1000 mb (surface level) to 20 mb (∼25 km). At low levels (1000 mb, a in **Figure 5**) the wind is weak, with wind velocity between 5 and 10 m/s. The prevailing direction at this altitude is location dependent, and has a NE component over Vatnajökull. Higher up in the atmosphere the wind field is more uniform and less affected by the surface topography. The velocity tends to increase moving higher up to 300 mb (9–10 km, f in **Figure 5**). Further up the wind velocity decreases to increase again at 20 mb (i in **Figure 5**) where there is a clear W-E directionality, with Easterly winds weaker than Westerly ones.

## Impact of Tephra on Infrastructure

Tephra can have wide type of impacts either nearby and far away the volcano. As reported in literature, at ground level volcanic tephra can cause:

FIGURE 8 | Probabilistic hazard map for tephra loading higher than 1000 kg/m2 given an eruption at Öræfajökull as 1362 AD. The isopachs reconstructed by Thorarinsson (1958) are reported as a reference with black lines.

TABLE 4 | Likelihood of receiving specific amounts of tephra is calculated for five main towns in the country.

Composition of the tephra, its grain-size distribution and presence of precipitation might enhance some of these hazards, as for example roof collapse conditions, damages to electrical infrastructure and contamination of water and vegetation. Wet tephra can reach higher load due to the contribution of rain that remains trapped in the deposit (Macedonio and Costa, 2012). This means that tephra fallout might have a different impact on buildings if it rains during or shortly after the eruption. Similarly, wet conditions might affect the conduction properties of tephra enhancing its effect in flashover events (Wilson et al., 2012). Finally silicic tephra can be strongly toxic for humans, pollute water supplies and poison grazing animals (Thorarinsson, 1979; Cronin et al., 2003).

Here, we investigate the potential impact of tephra fall on three key infrastructures, because if damaged, it can result in prolonged disruptions to the local population. The primary/direct hazards due to tephra fall onto airports, roads and power-lines are here considered.

#### Roads

Tephra on roads is an issue for traffic safety. It can result in the reduction of tire friction, obscure road markings, cause blockage of engine air intake filters and reduce the visibility for drivers. Few studies have properly investigated the thresholds of tephra fall capable to trigger critical driving conditions. Recently, Blake et al. (2017) presented results of laboratory experiments on how the properties of tephra fall (thickness, particle composition, particle size) along with precipitation can compromise skid resistance.

These results showed that for a 1 mm-thick fallout the skid resistance is below the "difficult site" safety levels. In case of thicker deposit (>5 mm) the conditions are more favorable, with skid resistance above the "safe" value. The study concludes that tephra fall thickness of 5 mm is the critical limit for which mitigation actions need to be taken to guarantee safe driving conditions on roads. At the same time Blake et al. (2017) presents a list of road disruptions occurred at volcanoes worldwide and shows how in several cases deposit several-10s-of-mm-thick have also been causing difficulties to the ground transportation. In this study we use a threshold value of 3 mm to represent disruptive conditions on the roads.

#### Airport

As reported in Guffanti et al. (2009) and Prata and Tupper (2009) the primary hazard to airports is ashfall, which can cause loss of visibility, create slippery runways, infiltrate communication and electrical systems, interrupt ground services, and damage buildings and parked airplanes. The skid resistance analysis perfomed by Blake et al. (2017) can be partly applied to airfield and runways, even though no clear thresholds exist for this environment with each airport operating authorities responsible for mantaining the runways functional and secure. Some critical conditions described for roads can be also applied to airfield and runways with few-mm ash deposit- to be considered a condition that can be critical for safe operations. Here, we use a threshold value of 1 mm to investigate the impact on airports.

#### Power Lines

Several papers discuss the vulnerability of electrical instrustructure to volcanic ash (Bebbington et al., 2008; Wardman et al., 2012; Wilson et al., 2012) and are used to identify correlation between specific thresholds of ash thickness and the level of impact (low, medium, and high) it would have on electrical infrastructure (Wilson et al., 2012). Critical infrastructure as power-lines have been investigated for damages due to:

The effect of tephra fall on this type of infrastructure depends on three main factors: the grain-size, the tephra load and presence of precipitation (wet deposition).

In this study we only look into those critical conditions associated with fine tephra because, as shown in **Figure 3**, this is the assumption we use for an eruption at Öræfajökull.

**Table 3** summarizes the results as reported in the mentioned paper for fine tephra. For a deposit larger than 5 mm in wet conditions there is a high likelihood of insulator flashover. Damages to tower, poles and lines are highly likely for thicknesses larger than 100 mm in dry conditions and 5 mm in wet conditions. In the following analysis we have been considering the likelihood to get deposit larger than 10 and 100 mm in dry conditions. The assumption is that these limits can be considered as of danger in wet condition for a fully wet deposit (i.e., all the voids in the deposit are full of water).

## RESULTS

## Probabilistic Hazard Maps

We performed numerical simulations of volcanic tephra dispersal by using VOL-CALPUFF code. The simulations yielded tephra thicknesses and concentrations data points over Iceland. We then generated probabilistic hazard maps for tephra loading at given thresholds by adopting a Monte Carlo approach. Here, we present three probabilistic hazard maps (**Figures 6–8**) for tephra loading thresholds of 1.0 kg/m2 (equivalent to about 0.1 cm under the assumption of a deposit of 1000 kg/m3 ), 100 kg/m2 (equivalent to about 10 cm) and 1000 kg/m2 (equivalent to about 1 m).

Comparison of these three maps reveals that the higher deposit load threshold investigated, the smaller is the area likely interested by this load. This is because in each simulation the thick deposit impacts a proximal area around the volcano, whereas thinner deposit can spread over a larger domain. In addition, in all maps there is a general eastward trend that reflects the meteorological statistics, as shown in **Figure 5**. At high altitudes the stronger westerly wind dominates the distribution of tephra toward the East.

In **Figure 6** the isoline of 0.1 cm thickness produced by Thorarinsson (1958) falls entirely within the 25–50% probability of reaching this deposit thickness. All the country, apart from the northern part of the Westfjords, has at least a 5% probability of exceeding a ground load threshold of 1 kg/m2 .

In **Figure 7** the isoline of 10 cm produced by Thorarinsson (1958) falls entirely within the 75–100% probability of reaching this deposit thickness. About half of the country has at least a probability of 1% of exceeding a ground load of 100 kg/m2 . This area includes eastern and South-Eastern

FIGURE 11 | Map showing potential tephra fallout impact to power-lines in case of a 1362-like eruption at Öræfajökull.

Iceland, most of the highlands and parts of the North– East.

In **Figure 8** the likelihood of exceeding a threshold of 1000 kg/m2 is confined to a small area closest to the volcano edifice. Part of the proximal area enclosed within the 20 cm isoline falls within the 5–25% probability of reaching a deposit thickness of 1 m.

The potential impact of tephra fallout can also be investigated at specific locations by performing an analysis site by site. For example **Table 4** shows the likelihood of exceeding specific tephra ground loads at five key population centers in Iceland (**Figure 1**). All of these population centers have a likelihood higher than 5% of experiencing tephra fallout deposit exceeding a thickness of 1 mm. Höfn (97.8%) and Egilsstaðir (61.4%) are the most vulnerable centers as they are closest to the volcano. Höfn can experience very high loads (10 and 100 kg/m2 ) and the likelihood is about 95 and 73%, respectively.

Reykjavík is the less exposed town due to its distance from the volcano and mainly because it is located upwind the dominant wind direction. Akureyri (in the North), despite its long distance, is showing intermediate values.

#### Impact Maps

An analysis to investigate how these results translate into an impact assessment has been done by quantifying the level of disruption due to tephra fallout on roads, airports, and powerline network. The type of disruption has been defined with regards to: (1) kilometers of roads potentially exposed to critical driving conditions, (2) number of airports potentially affected by tephra fallout on the runways, and (3) kilometers of powerline potentially affected by tephra fallout that can cause either flashover and physical damage to the lines. For each of this infrastructure a map has been created and the results are shown in **Figures 9–11**.

When looking at the possible disruption to road traffic (**Figure 9**) we can see that up to 268 km of the main road system (either paved and unpaved according to the National Land Survey IS50V database v3.4, 2012) could be affected by tephra fallout thicker than 3 mm with a likelihood between 75 and 100% (red road sector). There is a 25% likelihood that more than 2000 km of the road network can be affected by this condition (the yellow, orange, and red road sector). The main road extending from Vík, in the

South, to Húsavík, in the North, passing through the East part of the country will be impassable (see **Figure 1** for site locations).

The analysis performed for the airports network includes all the main airports in the country for which flights are regularly scheduled (as reported from the National Land Survey IS50V database v3.4, 2012). The analysis reveals that the airport in Hornafjörður will be affected by a deposit thicker than 1 mm with a probability between 75 and 100%. The likelihood for two more airports (Egilsstaðir and Vopnafjörður) to receive this amount of tephra is between 50 and 75%. These three airports are located in the Eastern sector of the country. The international airport of Akureyri shows a likelihood between 25 and 50%, whereas the domestic airport in Reykjavík is calculated to have a likelihood between 5 and 25%. The International airport of Keflavík has potential to be affected by tephra fallout from Öræfajökull volcano with a likelihood between 1 and 5%. All the main airports of Iceland (14 in total) have >0% likelihood to be affected by tephra fallout capable to create disruption to the operations of landing and taking off.

**Figure 11** shows a zoomed domain around Öræfajökull, as the main impact on the power line network is assessed to be quite proximal to the volcano. Only the network owned by the principal Icelandic provider Landsnet, which operates >33 kV powerage line, has been considered here. For this investigation a threshold of 100 mm of ash has been adopted, as given in **Tables 3**, 4. We assume this condition, corresponding to a load of about 100 kg/m2 , to correspond to critical conditions for powerlines. The results show that up to 115 km of powerline network will be exposed to such a load with likelihood of 75 to 100%. This is the part of the powerline passing nearby the volcano at a minimum distance of about 9.5 km (red line sector). Additional 45 km will be exposed to such load with a likelihood between 50 and 75%. So that overall more than 160 km of the powerline network can be damaged by an eruption at Öræfajökull and the likelihood for such an occurrence is greater than 50%. Almost the entire powerline network feeding the Öræfi district is vulnerable to these effects.

## DISCUSSION AND CONCLUSION

An eruption at Öræfajökull is likely to have a significant impact in Iceland, its nearby Nordic countries and, Europe. As evidence shows, such an eruption can produce pyroclastic density currents that will inundate the communities at the base of the volcano, glacial outbursts can also cover a substantial area (possibly >300 km2 ) around the volcano. On top of that a heavy tephra fallout is expected to occur over a period of about 1 day (e.g., Thorarinsson, 1958; Jónsson and Valdimar, 2007; Thordarson and Höskuldsson, 2007; Gudmundsson et al., 2008).

The volcanological scenario selected for this study represents the worst case scenario known to have happened at Öræfajökull, the 1362 AD eruption. The model results and the probabilistic hazard maps clearly indicate how, in case of a similar event, the associated tephra fallout might have effects all over the country. However, the location of the volcano, in the southernmost tip of the island, together with the prevailing westerly wind at high altitude favor the tephra fallout and ash dispersal to occur mainly toward the East. This will possibly represent a serious issue for the air traffic over the European air space, as occurred in 2010 when Eyjafjallajökull volcano erupted (Reichardt et al., 2018).

This study highlights two main issues: firstly, there are no places in the country completely safe from receiving ash generated from an eruption like 1362 at Öræfajökull (**Figure 6**); second, the tephra fallout can have a very severe impact in the proximity of the volcano with up to 1000 kg/m2 of ash expected up to a distance of only 25 km the vent (**Figure 9**). Most of the main towns in Iceland have likelihood higher than 1% to receive an amount of ash larger than 10 kg/m2 (**Table 4**).

The impact analysis performed for three different types of infrastructures (roads, airports, and power lines) reveals the vulnerability of the country in case of such an eruption. **Figure 9**, which shows the kilometeres of road network potentially affected by critical driving conditions, suggests that the entire Eastern part of the country will be hardly reached by car either through the Southern and the Northern ring road section. The main town of Egilsstaðir and the very popular localities of Höfn, Skaftafell and Jökulsarlón will be cut off the main viable connections, with important implications for either inhabitants and tourists potentially trapped in the area due to very low visibility conditions and unsafe driving conditions of paved roads.

The results show also that the likelihood for the airports in Hornafjörður (next to Höfn) and Egilsstadir to be disrupted by the tephra fallout is >50% (**Figure 10**). Consequently, connection with the Eastern part of the country via air or road will be very difficult and dangerous during and shortly after the eruption. An eruption like 1362 will affect the communication and commuting between the East part of the country and the Greater Reykjavík area (i.e., the capital and its surroundings).

Failure of the electricity provision can be expected due to damages to powerlines during and shortly after the eruption (**Figure 11**). These data need to be seriously evaluated when planning for mitigation actions and evacuation plans. The impact of a similar scenario would be even more dramatic if an eruption will take place during the winter time, when the daylight time is very short and the need for electricity to illuminate is higher and essential for the daily activities of the society and its economy.

The results shown so far are just the beginning of a long study that is needed to fully address the volcanic hazard at Öræfajökull volcano, but still enough to identify the potential critical scenario that could raise if it will erupt again with a similar intensity as in 1362. The unrest phase declared in November 2017 has been putting emphasis on the need of a proper multi-hazard assessment for Öræfajökull and this study is a contribution to it.

## AUTHOR CONTRIBUTIONS

feart-06-00196 November 13, 2018 Time: 12:17 # 15

SB worked on the definition of the volcanological scenarios, performed the simulations, processed the data and wrote the interpretation, and discussion of the results. DDR worked on the volcanological scenarios and performed part of the simulations, processed the data and contributed to the discussion of the results. DDR, BBB contributed to the interpretation of the results, used GIS to recalculate the masses, and prepared all the graphical maps. TT worked on the definition of the volcanological scenarios and contributed to the discussion of the results. SK contributed to the discussion and interpretation of the results.

## REFERENCES

## FUNDING

This study has been benefiting from the GOSVÁ project (2012–2018) funded by the Icelandic Government together with the International Civil Aviation Organization (ICAO).

#### ACKNOWLEDGMENTS

A great thank goes to the Steering Commitee of the project for fruitful and constructive discussion. Marco Viccaro, the reviewer and the Chief Editor Valerio Acocella are kindly thanked for their help in improving the paper.

a multiphase rhyolitic eruption at Tarawera, New Zealand. J. Geophys. Res. 110:B03203. doi: 10.1029/2003JB002896

Smellie and M. G. Chapman (London: Geological Society, Special Publications 202), 81–90.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Barsotti, Di Rienzo, Thordarson, Björnsson and Karlsdóttir. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Paroxysmal Explosions, Lava Fountains and Ash Plumes at Etna Volcano: Eruptive Processes and Hazard Implications

Sonia Calvari 1 \*, Flavio Cannavò1 , Alessandro Bonaccorso1 , Letizia Spampinato1 and Alessandra G. Pellegrino1,2

1 Istituto Nazionale di Geofisica e Vulcanologia - Osservatorio Etneo (INGV-OE), Sezione di Catania, Catania, Italy, 2 Dipartimento di Scienze Biologiche Geologiche e Ambientali, Università degli Studi di Catania, Catania, Italy

#### Edited by:

Roberto Sulpizio, Università degli Studi di Bari Aldo Moro, Italy

#### Reviewed by:

Geoff Kilgour, GNS Science, New Zealand Susanna Jenkins, Earth Observatory of Singapore, Singapore

> \*Correspondence: Sonia Calvari<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Volcanology, a section of the journal Frontiers in Earth Science

Received: 06 April 2018 Accepted: 13 July 2018 Published: 15 August 2018

#### Citation:

Calvari S, Cannavò F, Bonaccorso A, Spampinato L and Pellegrino AG (2018) Paroxysmal Explosions, Lava Fountains and Ash Plumes at Etna Volcano: Eruptive Processes and Hazard Implications. Front. Earth Sci. 6:107. doi: 10.3389/feart.2018.00107 Lava fountains have a major impact on the local population since they cause ash plumes that spread several kilometers above and hundreds of kilometers away from the crater. Ash fallout is responsible for disrupting airports and traffic on the motorways well beyond the area of the volcano itself, as well as affecting the stability of buildings and causing public health issues. It is thus a primary scientific target to forecast the impact of this activity on local communities on the basis of parameters recorded by the monitoring network. Between 2011 and 2015, 49 paroxysmal explosive episodes occurred at two of Mt Etna's five summit craters: the New South-East Crater (NSEC) and the Voragine (VOR). In this paper, we examine the features of the 40 episodes occurring at the NSEC during 2011–2013, and of the 4 events at VOR in December 2015. We study these paroxysms using geophysical monitoring data, characterize the episodes, and analyse all available data statistically. Our main results are two empirical relationships allowing us to forecast the maximum elevation of the ash plume from the average height of the lava fountain, useful for hazard assessment and risk mitigation. For Etna, and using the examples described in this paper, we can infer that wind speed <10 m s−1 generally results in strong to intermediate plumes rising vertically above the crater, whereas wind speed >10 m s−1 is normally associated with weak plumes, bent-over along the wind direction and reaching lower elevations.

Keywords: Etna volcano, paroxysmal explosive activity, lava fountains, ash plume, thermal cameras

## INTRODUCTION

Lava fountains are a mixture of liquid clots and droplets (the pyroclasts) and magmatic gases (H2O, CO2, SO2, CH4, N2O, HCl, HF, and CO, in order of decreasing abundance; Allard et al., 2005). They are generated when large volumes of volatiles rapidly exsolve from the magma during its ascent and decompression along the conduit (Wilson et al., 1995). The mixture of gas and pyroclasts emerges from the vent as a jet (Wilson et al., 1995), expanding vertically into an eruptive column, where three main regions can be identified. The lowest portion is the gas-thrust region, and comprises the lava fountains (Sparks et al., 1997; Bonaccorso and Calvari, 2017).  Above the gas-thrust is the convective region, where atmospheric mixing occurs (Sparks et al., 1997; Bombrun et al., 2018), causing a decrease in ascent velocity and lateral wind transport (Carey and Sparks, 1986). The uppermost portion is the umbrella region, where the volcanic plume reaches the zone of neutral buoyancy and spreads out laterally (Carey and Sparks, 1986; Bonadonna and Phillips, 2003). Within the buoyant umbrella region, plume rise rates are less than 15 m s−1 (Patrick, 2007; Patrick et al., 2007). Volcanic plume dynamics depend closely on the interaction with the atmospheric wind field. They are classified as strong or weak depending on whether the characteristic plume rise velocity is greater or much smaller than the cross wind speed, respectively (Carey and Sparks, 1986). Thus, strong plumes vertically rise (**Figure 1A**), while weak plumes tend to bend over in the downwind direction (**Figure 1C**) which then start spreading sub-horizontally around the neutral buoyancy level (Carey and Sparks, 1986; Bonadonna and Phillips, 2003). The plume shape has significant effects on the forecast of the fallout area, which is ideally distributed all around the vent in the case of a strong plume and zero wind, and elongated downwind for a weak plume and strong wind. The stronger the wind, the narrower the resulting fallout region from the plume (Andronico et al., 2009). Interestingly enough, the wind speed has important effects not just on the fine-grained portion of the plume, but even on the bombs and ballistics, being able to transport large ejecta for several kilometers down wind, causing damage to cars, solar panels, and windows (Andronico et al., 2015). It is difficult to univocally forecast plume height because it is related to the conditions and properties of the erupting system, such as magma discharge rate (Carey and Sparks, 1986), exsolved volatiles from magma, magma/water interaction (Schmith et al., 2018), and vent geometry (either elongated fissure or near-circular conduit; Wilson and Head, 1981; Wilson et al., 1995). Given that discharge rate and ultimately erupted volume have been measured from the lava fountain height (Calvari et al., 2011; Carbone et al., 2015; Bonaccorso and Calvari, 2017), there must be a direct connection between lava fountain height and maximum plume elevation, although this last parameter is also affected by wind speed (e.g., Bonadonna et al., 2015).

Paroxysmal explosive basaltic activity, especially when occurring at frequently erupting volcanoes, can have a strong impact on the population, air traffic and public health (Andronico et al., 2008, 2009; Martin et al., 2009; Baxter, 2010; Horwell et al., 2017). During 2011–2015, Mt. Etna volcano produced 49 paroxysmal lava fountain episodes, evolving into volcanic plumes that spread for more than 9 km above sea level (a.s.l.) into the atmosphere, and with ash fallout as far as hundreds of km from the vent (e.g., Dellino and Kyriakopoulos, 2003; Azzopardi et al., 2013; Athanassiadou, 2016). It is thus of paramount importance to monitor this activity and analyse the data with the aim of developing routines or tools that might be useful for hazard assessment and impact prediction (Cannavò et al., 2017).

This paper seeks to analyse the data recorded during Etna's recent lava fountaining episodes in order to find possible relationships between measurable parameters. We have used images from the network of fixed thermal and visual cameras installed on the volcano slopes by the Istituto Nazionale di Geofisica e Vulcanologia-Osservatorio Etneo (INGV-OE), and gathered from these several parameters such as lava fountain height, exit velocity, duration of the episode, erupted volume and plume height. Among these, estimates of tephra fall volume are of paramount importance for hazard assessment and risk mitigation during future explosive events, since they are typically based on eruption intensity (eruption rate) and magnitude (total erupted mass; Walker, 1980).

We present an analysis of 44 of the 49 paroxysmal episodes occurring at Etna between 2011 and 2015; the remaining 5 could not be analyzed due to poor visibility from the ground cameras. We have used thermal and visual images recorded by the monitoring network of fixed cameras installed by INGV-OE (Andò and Pecora, 2006; Coltelli et al., 2017) to gather information about the total erupted mass and eruption rate of each episode. We have statistically analyzed the dimensional parameters extracted from the monitoring cameras, such as height of the lava fountain and ash plume vs. time, duration of each episode, erupted volume of pyroclastics (obtained from the camera images) and of lava flows (taken from available published data), time-averaged discharge rate (TADR, here referred only to the pyroclastic portion of the erupted material; Harris et al., 2007) and maximum exit velocities of both the lava fountain and plume. Following Harris et al. (2007), instantaneous effusion rate (IER) is the discharge rate from a vent averaged over a smaller unit of time when compared to TADR, which we consider as averaged over the whole duration of the explosive event. We carried out a rigorous correlation analysis on the measured parameters and tried to find a link among parameters.

## ETNA'S ERUPTIVE ACTIVITY

Before 2010, lava fountain episodes at Mt. Etna (**Figures 1**, **2**) were an exceptional activity, normally preceding major flank eruptions, such as those in 1991-93, 2000-01, and 2002-03 (Calvari et al., 1994, 2001; Andronico et al., 2005; Andronico and Corsaro, 2011). However, between January 2011 and December 2013, a sequence of 44 such events occurred from a pit on the east flank of the South-East Crater (SEC), one of the four summit craters of Etna (**Figure 2b**). This activity built up a new summit cone (the New South-East Crater, NSEC, **Figure 2b**; Behncke et al., 2014; De Beni et al., 2015), but no major flank eruptions have followed this activity yet. Bonaccorso and Calvari (2013) explained this apparent paradox by the balance between intruded and extruded magma, which at Etna volcano has been found steady during the last 3–4 decades (Wadge and Guest, 1981; Harris et al., 2011, 2012). This steady 2013 at 17:36.

output can result either in a few, high-output rate effusive phases, or in several, low-output rate lava fountain events (Bonaccorso and Calvari, 2013). Whatever the cause, it is clear that paroxysmal explosive activity at Etna is becoming increasingly common, thus requiring a greater effort to forecast its possible evolution.

Lava fountains are always associated with ash plumes, often disrupting the air traffic and driveability on the main highway along the east coast of Sicily, where frequent winds blow the ash and cause problems not only to traffic but also public health, agriculture, aviation, and building stability (Andronico et al., 2008; Martin et al., 2009; Baxter, 2010; Horwell et al., 2017). Eruptive columns and ash plumes, associated with the lava fountains, cause a more widespread fallout of ash transported well beyond the boundaries of Sicily and southern Italy and into Greece (Dellino and Kyriakopoulos, 2003) and Malta (Azzopardi et al., 2013).

## Lava Fountaining From NSEC, 2011–2013

The 44 lava fountain (LF) episodes at the NSEC between 2011 and 2013 were preceded by deep magma uprising (Del Negro et al., 2013). We distinguished three main phases, separated by two long eruptive pauses (**Table 1**). The first phase lasted from 11 January 2011 to 24 April 2012, and comprised 25 LF events displaying similar features. These are: duration of the Strombolian phase preceding the paroxysmal episode (hours to days); duration of the LF event (hours); maximum height of the LF (800–1,500 m); extension of the lava flow field associated with the explosive event (1–4 km; Behncke et al., 2014); and erupted volume of pyroclastics and lava flows (∼2.5 × 106 m3 ; Calvari et al., 2011; Gouhier et al., 2012; Bonaccorso and Calvari, 2013). The first eruptive pause lasted from April 2012 to February 2013 (**Table 1**). During this pause, only mild and intermittent eruptive activity within the Bocca Nuova summit crater (BN; **Figure 2b**) occurred between July and October 2012, producing a small cinder cone and lava flow field within its crater (Cannata et al., 2015; Slatcher et al., 2015; Spampinato et al., 2015). During the period preceding the second LF phase, ground deformation data from the permanent GPS network revealed an inflation of the volcano due to the pressurization of its plumbing system (e.g., Bruno et al., 2016), whereas gas fluxes showed that no volatile-rich magma was supplied to the shallow (1–2 km below sea level, b.s.l.) volcano feeding system (Spampinato et al., 2015). This suggested that the BN eruptive activity was produced by magma already residing within the shallow plumbing system (Spampinato et al., 2015). A depressurization of the plumbing system was observed when LF activity resumed again in 19 February 2013 (**Table 1**), giving rise to the second eruptive phase that lasted until 27 April 2013, producing 13 additional events (Cannata et al., 2015; Spampinato et al., 2015; Bruno et al., 2016). There was another eruptive pause between April and the end of October 2013, followed by six additional paroxysmal explosive events until the end of that year. This last and third phase was very different from the previous, and comprised the most explosive end-member of the series, occurring on 23 November 2013 (Bonaccorso et al., 2014; Andronico et al., 2015; Carbone et al., 2015; Corradini et al., 2016; Montopoli, 2016). This very short but violent paroxysm produced an exceptionally high LF of over 2.5 km above the crater, and was the only one at the NSEC that did not produce lava flows (Bonaccorso et al., 2014). The third phase also displayed a unique episode (26 October 2013, **Table 1**) during which two craters were erupting together, with the North-East Crater (NEC) producing an ash column that flanked that of the NSEC (Andronico et al., 2018, this volume). As a result of the proximal accumulation of ballistics from the LF around the vent, a new cinder cone (the NSEC) grew on the east flank of the SEC, with a volume of 19 × 106 m3 (Behncke et al., 2014) by April 2012, and growing up to 50 × 106 m3 by October 2014 (De Beni et al., 2015). The shallow storage feeding these paroxysmal events, imaged from volcanic tremor and long-period seismic events, is very shallow and located within the volcanic pile between sea level and the ground surface (Cannata et al., 2015) and should represent just the final portion of the uprising path. We know very little about the episode occurring in 2014 (**Table 1**) because of poor weather conditions.

## Lava Fountain From VOR, 3–5 December 2015

In January 2015, ash plume emissions started from the Voragine crater (VOR), after 16 years quiescence at this crater (Corsaro et al., 2017; Cannata et al., 2018). Explosive activity from this crater increased in intensity and became almost regular since 27 October, climaxing with a sequence of four paroxysmal events between 3 and 5 December (Vulpiani et al., 2016; Bonaccorso and Calvari, 2017; Neri et al., 2017). The four explosive events interrupted the usual trend of inflation observed at Mt. Etna (Aloisi et al., 2017). They were characterized by decreasing LF maximum heights (from 4.1 to 1.1 km), by lack of associated lava flows, by durations varying between 65 and 114 min, and by a total erupted volume of at least 10 × 106 m3 dense rock equivalent (Bonaccorso and Calvari, 2017). The eruption columns and ash plume rose vertically with velocities between 150 and 200 m s−1 , and the transition from the convective to the umbrella region was observed at ∼8 km height (Vulpiani et al., 2016). In about 10 min the ash plume expanded up to 10–15 km above sea level (a.s.l.; Vulpiani et al., 2016; Corsaro et al., 2017). Tephra fallout accumulated mostly within the summit crater, resulting in a fountain-fed rheomorphic lava of ∼7.2 × 106 m3 (Corsaro et al., 2017; Neri et al., 2017). The strain data constrained the depth of the magmatic source

(Continued)

**152**

at ∼ 1.5 km b.s.l. (Bonaccorso and Calvari, 2017). This is much deeper than the magmatic source at sea level inferred for the lava fountain episodes at the NSEC (Bonaccorso et al., 2013a,b; Cannata et al., 2015). The fast transfer from the deepest levels of the plumbing system of basic, undegassed magmas has been considered as the crucial triggering factor leading to the development of exceptionally violent volcanic phenomena (Cannata et al., 2018). Also the time scale for completely refilling the VOR system and renewing magma is very fast, being as short as a few tens of hours (Pompilio et al., 2017).

## CAMERA RECORDED DATA

We have used the thermal and visual cameras from the INGV-OE monitoring network (**Figure 2c**) to acquire information about the duration of the eruptive activity, such as the start and end of the Strombolian activity preceding the paroxysm, of the lava fountain, of the lava flows, exit velocity, as well as the height of lava fountains. **Table 2** reports the characteristics of the cameras used in this study. Additional information can be found in Andò and Pecora (2006) and Coltelli et al. (2017). Details on the time and kind of eruptive activity for each episode, as well as its duration, are listed in **Table 1**, which also contains the erupted volumes of pyroclastics for each event. These have been obtained, following the approaches of Calvari et al. (2011) and Bonaccorso and Calvari (2017), by integrating the height of the LF obtained from the saturated portion of the thermal images over time at 1 min time lapse, and using the formulas:

$$\mathbf{U} = \left(2\mathbf{g}\mathbf{H}\right)^{0.5} \tag{1}$$

$$\mathbf{V} = \mathbf{U} \times \mathbf{A} \times \mathbf{t} \tag{2}$$

where U is the mean fluid exit velocity at the vent, H is the LF height, g is the gravity acceleration, V is the fluid volume (gas + pyroclastics) erupted by a LF, A is the vent surface area, and t is the duration of the LF. Vent surface area is calculated considering a circular vent with a diameter of 30 m (Calvari et al., 2011) and supposed constant. Considering a circular shape for the vent area, and multiplying the mean fluid exit velocity at the vent by the vent surface area integrated over time, we transformed the 2D view given by the camera obtaining a 3D volume estimation. The images from the thermal cameras (EMOT, ENT, EMCT, and MBT; **Figure 2c**, **Table 2**) have been used to measure the maximum height of the LF, considered as the saturated portion of each image. Saturation varies with respect to distance between 50◦C for the most distant camera (ENT, 15 km; **Figure 2c**) to 100◦C for the closest (EMOT, 3 km; **Figure 2c**). The error in the height measurement is ±50 m. We carefully measured the height of the saturated portion, avoiding errors due to clouds interference and to the variable position of the eruptive vent, given that the NSEC, responsible for this eruptive activity, grew by more than 190 m during the lapse of time considered here (Behncke et al., 2014). The heights of the LF at 1 min intervals have then been used to calculate the total fluid erupted volume from Equations (1) and (2), which includes both gas and pyroclastics (Calvari et al., 2011; Bonaccorso and Calvari, 2017). To obtain the volume of pyroclastics from the total volume

§Data averaged from Ganci et al. (2012) and Behncke et al. (2014).

Bonaccorso

and 16 March 2013 (De Beni et al., 2015) divided by the number (8) of eruptive episodes. \*\*Data from Scollo et al. (2014). §§Data from Scollo et al. (2015). ◦◦Data from

from

Vulpiani et al. (2016).

=Data from this work. Gray lines indicate poor weather conditions and lack of visibility from cameras.

 and Calvari (2017).

+Data from De Beni et al. (2015).

⋄Data from

Behncke et al. (2014). ++Derived from the total volume erupted between 19 February

Corradini et al. (2016). #Data

> Guerrieri et al., 2015." " Data from

TABLE 2 | Characteristics of visual and thermal cameras of the INGV-OE monitoring network used in this paper.

See Figure 2c for their location.

of erupted fluid (gas + pyroclastics), we have considered 0.18% as the ratio between the volumes of magma and fluids within the eruptive column as typical for Etna's fountains (e.g., Donnadieu et al., 2016). This 0.18% between the pyroclastic portion of the LF and the total erupted fluids (gas + pyroclasts), derives from the average volumetric growth of the NSEC cone during the first 12 paroxysmal episodes. In fact, most of the LF fallout accumulates all around the vent, causing the fast growth of the NSEC cinder cone. The NSEC cone volume was measured by Behncke et al. (2014) on 31 August 2011, after 12 LF episodes from the NSEC, at ∼9.5 × 106 m3 dense rock equivalent (DRE). This volume is ∼0.18% of the total volume of fluids (gas + pyroclastics, ∼5.4 × 109 m3 ) erupted during the first 12 LF events and obtained from Equations (1) and (2). Considering this a constant, we applied the 0.18% to all pyroclastic volume calculations of the LF episodes from the NSEC reported in **Table 1**. For the four paroxysmal episodes of VOR, this value has been calculated individually for each event, showing a significant variability and a general gradual decline (from 0.22 to 0.09%; Bonaccorso and Calvari, 2017). The maximum error on the pyroclastic volume calculation is reported in **Table 1**, and is calculated from the LF height error (±50 m) and its propagation when using Equations (1) and (2).

The lava flow volumes in **Table 1** are taken from Ganci et al. (2012), Behncke et al. (2014) and De Beni et al. (2015), or are averages of more values when available. TADR is related only to the explosive component of each episode and averaged over the whole duration of the paroxysmal explosive event.

Note that the starting and ending times in **Table 1** (in UT format) have been obtained from careful analysis of the images of the closest cameras (3.0 km from the vent, see **Table 2**), and thus refer to surface expression of the eruptive activity. As such, they could be significantly different from those obtained by geophysical data that refer to changes detected in the source region (e.g., Carbone et al., 2015; Spampinato et al., 2015; Bonaccorso and Calvari, 2017). In addition, given the change from Strombolian to LF activity is generally gradual and characterized by a transitional phase (e.g., Spampinato et al., 2008), we have taken as the starting time of the LF activity the moment when tephra fallout began spreading well outside the crater, in a way to minimize errors in volume calculations. In fact, we have observed that during the transitional phase, the jet of magma is highly collimated, and tephra fallout mainly falls back within the crater, with very little contribution to the volume accumulated on the outer flanks of the cone and beyond.

For four of the 44 eruptive episodes at NSEC between January 2011 and December 2013, we lack data because poor weather conditions hampered camera visibility throughout the entire episode. These are: 8 October 2011, 21 February 2013, 11 November 2013, and 2 December 2013 (**Table 1**). However, in the next section we also considered these episodes when analyzing the relationships between pause interval and total erupted volume. When instead cloud coverage was partial, we have interpolated the missing values. In addition, the 19 September 2011 event (**Table 1**) allowed visibility for only 7 min out of 65, thus in this case most of the erupted volume has been extrapolated on the basis of available data and must be considered a minimum, given that the climax of the LF has been lost. We have analyzed 6333 frames to extract dimensional parameters of the LF. Of these, 352 have been interpolated because of cloud cover, making a 5.6% data interpolated. The amount of interpolated frames for each LF is indicated in **Table 1**. However, it is worth noting that the two events of 18 April 2013 and 20 April 2013 occurred in poor weather conditions, allowing only 1/3 and 2/3 of their duration to be captured by the cameras. These two volumes must be considered as minimums. **Figure 3** displays the first LF episode from NSEC (A) occurring on 12 January 2011 and (B) the first episode from VOR occurred on 3 December 2015, with the outline of the volcano and the two vent positions. Note the very different size of both the LF (in white-red) and eruptive plumes (bluish to pink).

The ECV visual camera (**Figure 2c**, **Table 2**) has been used to extract the growing height of the eruptive plume, following the method proposed by Scollo et al. (2014). They calibrated the ECV images by projecting the plume on a vertical plane passing through the craters and parallel to the wind direction, taken from the daily wind direction forecast, thus obtaining a grid for the plume height estimation with an estimated error on the height of ±0.5 km. The ECV camera has a field of view allowing a maximum height of 9 km a.s.l., thus higher plumes are not visible. Fifteen out of the 44 episodes occurred during daytime and in good weather conditions, also enabling us to analyse the plume height growth together with LF height growth. These data are displayed in **Figure 5**, with the limit of 9 km a.s.l. due to the camera field of view. Note that **Table 1** reports maximum plume heights, thus most of the values greater than 9 km are taken from existing literature.

In summary, we have analyzed the images of 40 eruptive episodes from the NSEC and 4 from VOR in order to obtain:

22:05 UT, and (B) the first lava fountain episode from VOR on 3 December 2015 and taken at 02:53 UT. The red circles show the position of vent, and the red dotted line is the volcano profile. The white vertical bar in A is for scale. The lava fountain portion is in white-red, the ash plume in bluish-pink. Note the different size of the two episodes, with the VOR lava fountain and plume being much taller and wider.

(1) the duration of the Strombolian activity preceding each episode; (2) the duration of the LF; (3) the total volume of erupted pyroclastics; (4) the time-averaged discharge rate (TADR), referred just to the LF portion of each episode and averaged over the entire duration of the paroxysmal explosive phase; (5) the average and (6) maximum LF heights; (7) the maximum exit velocity of LF; (8) the duration of the lava flow output; (9) the plume type; (10) the maximum elevation of the plume (when lower than 9 km) and the plume growth for 15 events displayed in **Figure 5**; (11) the maximum exit velocity of the plume. Plume type is differentiated into a strong plume when vertical, intermediate when slightly inclined (∼70◦ ) downwind, and weak when markedly bent in a downwind direction (**Figure 1**). These results are listed in **Table 1**. Lava flow volumes and maximum plume heights are taken from available published data as specified in **Table 1**.

## STATISTICAL APPROACH

To make the most of the substantial set of paroxysms, we statistically characterized the measured volcanological features with the final aim of analyzing possible correlations among them. Indeed, conducting a correlation analysis that assumes the data follow a normal distribution when, in fact, the data are nonnormal, can lead to inaccurate results (e.g., linear correlation analysis is not robust in the presence of non-normality; Kowalski, 1972). To avoid this potential error, we first determined the right distribution of each considered feature. For the analysis, we took into account the following nine features, listed in **Table 1**:

Each feature was tested for a normal distribution. We used the Shapiro-Wilk (SW) test that is a widely used test of normality in frequency statistics (Shapiro and Wilk, 1965). The null-hypothesis of this test is that the population is normally distributed. The Shapiro-Wilk test is based on the correlation between the data and the corresponding normal scores and performs better than the Kolmogorov-Smirnov (KS) test in recognizing samples from non-normal distributions (Steinskog, 2007). We adopted a confidence level for the test of 5%.

For all the considered features, the test rejected the nullhypothesis of normality, thus we estimated the best distribution among the following ones: t location-scale (Rinne, 2011), inverse Gaussian (Chhikara and Folks, 1989), exponential (Balakrishnan and Basu, 1996), logistic and log-logistic (Balakrishnan, 1992). To this set of alternative distributions, we also added the normal distribution to confirm the SW test outcomes. Each distribution was fitted using maximum likelihood estimation. The Akaike information criterion (AIC) was then applied to choose the distribution that better describes the data. The AIC estimates the relative quality of statistical models for a given dataset (Burnham and Anderson, 2003). Given a collection of models for the data, AIC estimates the quality of each model, relative to each of the other models. Thus, AIC provides a means for model selection. The AIC corrected (AICc) for finite sample sizes is defined as:

$$AIC\mathfrak{c} = -2\log\left(L\left(\hat{\theta}\right)\right) + 2k + (2k+1)/(n-k-1) \tag{3}$$

L being the model likelihood of the estimated model parameters (θˆ), k the number of parameters, and n the sample size. For all the estimated best distributions, we calculated the first and second moments when existing and finite. Results are reported in **Table 3**, where it is then possible to read the estimated characteristic scales of the features.

We performed a correlation analysis to infer possible hidden relationships among the considered features. Correlation is a bivariate analysis that measures the strengths of association between two variables and the direction of the relationship. In terms of the strength of relationship, the value of the correlation coefficient varies between +1 and −1. When the value of the correlation coefficient lies around ±1, then it is said to be a perfect degree of association between the two variables. As the correlation coefficient value goes toward 0, the relationship between the two variables will be weaker (Hazewinkel, 2001). Commonly, Cohen's standard is used to evaluate the correlation coefficient to determine the strength of the relationship, or the effect size, where correlation coefficients between 0.10 and 0.29 represent a small association, coefficients between 0.30 and 0.49 represent a medium association, and coefficients of 0.50 and above represent a large association or relationship (Cohen, 1988).

There are different types of correlation methodologies. Pearson correlation (parametric test) is the most widely used correlation statistic to measure the degree of the relationship between linearly related variables (Pearson, 1895). For the Pearson correlation, both variables should be normally distributed. Other assumptions include linearity and homoscedasticity. Linearity assumes a straight-line relationship between each of the variables in the analysis and homoscedasticity assumes that data are normally distributed about the regression line (and both variables have a finite variance). These assumptions make the Pearson coefficient unreliable in case of non-normally distributed data. In order to mitigate the distribution effects, we used the Box-Cox transformation (Box and Cox, 1964) to convert the feature distributions into normal ones. In addition, we used the non-parametric tests of Spearman rank (Spearman, 1904) and Kendall rank (Kendall and Gibbons, 1990) correlations to measure the degree of association between two variables without any assumptions about the distribution of the data.

Based on the statistical analysis of the parameters obtained from analyzing the camera images, considering all the three cited coefficients, the correlation analysis confirms that the only significant/large associations are those for the self-explanatory feature pairs:

Nevertheless, it is important to underline that the above statistical analysis is valid for the considered set of volcanological data and should not be generalized outside the considered time interval without further accurate studies.

## LAVA FOUNTAIN AVERAGE HEIGHT VS. PLUME MAXIMUM HEIGHT

We investigated the relationship between the average height of the LF (HF) and the maximum height that reaches the plume (HP). The average heights of the LF are derived from the thermal images at 1 min intervals. Regarding the measurement of the maximum heights of the plume, in the current literature the work presenting the most complete and detailed catalog is the one published by Scollo et al. (2014), which reports the height measurements of the plume of 19 LF from 10 April 2011 to 27 April 2013. For each event considered, in their **Table 1** the Authors reported three heights respectively estimated from: (1) ECV camera, which, however, cuts the image to 9 km a.s.l.; (2) the MODIS satellite; (3) the SEVIRI satellite. The heights measured with greater precision are those seen by ECV, which records the peak in elevation when below 9 km. Then, with less precision, is the value obtained by MODIS, which has less frequency of satellite passages (and therefore fewer images; 1 every 15 min with 1 km of spatial resolution in the thermal infrared; Corradini et al., 2016). Finally, with even less precision, is the value obtained by SEVIRI, which has a spatial resolution of 3 km and records images every 15 min (Ganci et al., 2012). To obtain a more robust dataset, we therefore used the values obtained from ECV measurements, when HP was less than 9 km in height, and for other cases (HP > 9 km) we considered the value from MODIS and alternatively, if this was not available, the value of SEVIRI (**Table 4**). For a few cases (**Table 4**) where ECV camera images saturated at 9 km a.s.l., but with the satellite passage providing a lower height because their passage was temporarily distant from the LF climax, we used the ECV camera height. We did not consider the 19 September 2011 and 08 October 2011 episodes due to cloudy weather. Moreover, for the LF of 12 January 2011, we used the plume height reported in Calvari et al. (2011); for the LF of 23 October 2011, we considered the plume height measured by SEVIRI and reported in Guerrieri et al. (2015); for the LF of 23 November 2013, we considered the height inferred by Corradini et al. (2016) by using a multi-sensor approach considering different satellite instruments such as Meteosat, MODIS and SEVIRI, and finally for the four events of the VOR of 3–5 December 2015, we used the plume heights reported in Vulpiani et al. (2016), obtained by three radar systems (at C and X band) operated by the Italian Department of Civil Protection (**Table 4**).

The relationship between HF and HP has some limitations. In fact, it considers a heterogeneous set of data because heights from ECV camera saturate at 9 km a.s.l., and in these cases it is necessary to consider this underestimated value or to resort to the satellites that have a greater resolution error, when often their passage in time is not perfectly synchronous with the maximum height reached by the plume. Despite these limitations, the relationship presents an interesting correlation (**Figure 4**) with a good R2 of 0.78 for the overall average (Equation 4), and 0.91 for the VOR events only (Equation 5). This provides a useful tool to make a rough forecast of the maximum height reached by the plume once the fountain, well after its initial phase, stabilizes its explosive power and thus its average height, which can be accurately measured in real-time by the cameras monitoring system.

The equations expressing these relationships are:

$$\mathbf{y} = 5.26\mathbf{x} + 6, 830 \text{ R}^2 = 0.78 \tag{4}$$

$$\mathbf{y} = 2.00\mathbf{x} + 11,286 \text{ } \mathbb{R}^2 = 0.91 \tag{5}$$

This equals to:

$$\text{Hp} = 5.26 \text{Hf} + 6,830 \text{ (NSEC)} \tag{6}$$

$$\text{Hp} = 2.00 \text{Hf} + 11,286 \text{ (VOR)} \tag{7}$$

#### TIME-DEPENDENT RELATION BETWEEN LAVA FOUNTAIN AND PLUME HEIGHTS

For 15 of the 49 eruptive events, occurring during daytime and with good weather conditions, we also managed to obtain the growing height of the plume over time from the analysis of the images recorded by the ECV camera (**Table 2**). This allows a view of the eruptive plume up to about 9 km a.s.l., estimated using the scheme proposed by Scollo et al. (2014).

The comparison between height of the LF and height of the ash plume is shown in the graphs of **Figure 5**. Although with the limit of the 9 km caused by the camera field of view, from the analysis of **Figure 5** it is possible to observe different behaviors. A first class of events comprises 7 LF (10 April 2011, 20 August 2011, 29 August 2011, 8 September 2011, 18 March 2012, 18 April 2013, and 4 December 2015; **Figure 5**), where the ash plume stays above the 9 km threshold for most of the duration of the LF. These are all strong plumes (**Figure 1A**) formed in wind conditions variable between 2.69 and 8.69 m s−1 and with TADR, calculated only for the explosive portion, between 96 and 157 m3 s −1 , and with the outsider event of 4 December 2015 which displayed the exceptional value of 561 m3 s −1 (**Table 1**). A second TABLE 4 | Episodes (first column) and dates (second column) of lava fountains of NSEC (events from 1 to 42) and VOR (events from 46 to 49) with the mean height of the fountains (HF ) as obtained from the cam analysis in this study (third column) and maximum height of the plume (HP) (fourth column) as reported in literature with the sensor utilized (fifth column).

class of events comprises the two events of 9 July 2011 and 12 August 2011 (**Figure 5)**, where the plume remained above 9 km for a short period of time (30 and 52 min, respectively) during

Calvari et al. (2011) for the event of 12 January 2011, and from Guerrieri et al. (2015) for the event of 23 October 2011. The four HP values of the VOR December

2015 events (red squares) are taken from Vulpiani et al. (2016). The five yellow diamonds show data taken from ECV, with the 9 km a.s.l. limit.

the climax reached by the lava fountains, and decreased soon after and while the LF was still on going. These two events are both strong plumes formed in wind conditions of 3.99 and 7.99 m s −1 , and having TADR of 110 and 99 m3 s −1 , respectively, thus comparable to the range of the first class. A different pattern is shown by the 6 events of 4 March 2012, 12 April 2012, 3 April 2013, 12 April 2013, 27 April 2013, and 23 November 2013, where the height of the ash plume increases and decreases together with the height of the LF (**Figure 5**). It is worth noting that the 27 April 2013 plume decline after the end of the LF is not reported in **Figure 5** because it has been observed from thermal images and not from the ECV camera, given that it occurred during the evening and hence with no visibility from this camera. These are intermediate to weak plumes (**Figures 1B,C**) formed in wind conditions comprised between 3.22 and 21.78 m s−1 and displaying TADR between 23 and 171 m3 s −1 . Thus, from an overall view of these data, it appears that conditions of strong wind (above ∼10 m s−1 ) or low TADR (below 100 m3 s −1 ), or a combination of these two factors, might favor weak to intermediate plumes that decrease in height soon after or together with the end of the LF. Conversely, conditions of high TADR (above 100 m3 s −1 ) and slight wind (below ∼10 m s−1 ), or a combination of these two factors, involves the formation of strong plumes, with plume heights above 9 km maintained even during the phase of waning LF height (**Figure 5**).

The effect of wind speed on the plume style is clearly visible when comparing the LF events on 12 May 2011 and on 19 July 2011 (**Figure 6**). It is worth noting here that wind speed is sampled every 3 h, thus we lack details of wind variability during an event, and we consider it constant throughout a paroxysmal event. The first event occurred in low wind conditions (7.13 m s −1 ), with LF maximum height of 880 m, and average of ∼300 m, and with a TADR of 135 m3 s −1 (**Table 1**), allowing the gradual development of a strong plume rising vertically from the vent. As the eruptive phase proceeded, we observed at first, between 00:00 and 01:02, intermittent minor ash emissions (**Figure 6A**) accompanied by a Strombolian to transitional phase, followed between 01:02 and 01:17 by the formation of a weak plume (**Figure 6B**) inclined westward (left) by about 45◦ ; between 01:18 and 01:26 (**Figure 6C**) it passed to an intermediate plume with inclination of ∼70◦ westward; and from 01:26 to 3:21 giving rise to a strong plume rising vertically from the vent (**Figure 6D**) and spreading laterally (**Figure 6E**); between 3.21 and 3:44 the plume again became transitional (**Figure 6F**) until the end of the eruptive event. An opposite example is given by the event of 19 July 2011, which had a similar TADR (121 m3 s −1 , **Table 1**) and a slightly greater vertical extension of the LF height (1,040 m maximum, 487 m average; **Table 1**), but whose plume formed in conditions of stronger wind (15.05 m s−1 ). The activity started as pulsating and increasing ash emissions between 00:00 and 00:29 (**Figure 6G**), passing to weak plume (bent about 45◦ eastward, **Figure 6H**) between 00:29 to 01:38 (with LF about 800 m high), and again at pulsating and decreasing ash emissions between 01:38 and 02:18 (**Figure 6I**), until the end of the event. Thus, **Figure 6** displays the importance of wind speed and IER in determining the plume style, passing from no plume to weak, intermediate and then to strong plume with increasing IER (from 43 to 67, 85 and 130 m3 s −1 , respectively) and in conditions of low wind (7.13 m s−1 , **Table 1** and **Figures 6A–E**). Vice-versa, the opposite change occurs when IER vanishes (**Figures 6E,F**, from 119 to 0 m3 s −1 ). On the contrary, if wind speed is strong (15.05 m s−1 in the case of the 19 July 2011 event, **Table 1** and **Figures 6G–I**), then also a powerful LF ∼800 m high (with IER up to 158 m3 s −1 , **Figure 6H**) is unable to form a strong plume

(A,G) is for scale. See text for further details.

(**Figures 6G–I**). It is worth noting that the two events reported in **Figure 6** have very similar TADR and durations (averaged over the whole duration of the explosive event, 135 and 121 m3 s −1 , and durations of 8,220 and 7,800 s, respectively; **Table 1**) but wind speed is two-fold during the second case (7.13 vs. 15.05 m s −1 ; **Table 1**).

## DISCUSSION

In this paper, we have analyzed 44 paroxysmal explosive events occurring at Mt. Etna between 2011 and 2015, with the aim of analyzing the eruptive process and finding simple relationships among easily measurable parameters that might be useful for hazard assessment and risk mitigation. Using the images recorded by the network of thermal and visual cameras (**Figure 2c**, **Table 2**), we have characterized each event, extracting the parameters listed in **Table 1**. They are: the vent that produced the event (NSEC, VOR or NEC); duration of the Strombolian phase preceding the LF event; starting and ending time and duration of each LF event; volume of erupted pyroclastics and its maximum error; TADR; maximum and average height of the LF; maximum speed of the LF; the lava flow duration when associated to the explosive event; lava flow volume (taken from available published data, as detailed in **Table 1** caption); plume type (either strong, intermediate or weak, see **Figure 1**); maximum plume heights (from available published data, detailed in **Table 1** caption, or from our measurements using the ECV camera image); maximum velocity of plume growth; and wind speed.