Airspace Contamination by Volcanic Ash from Sequences of Etna Paroxysms: Coupling the WRF‐Chem Dispersion Model with Near‐Source L‐Band Radar Observations

1. Introduction

1.1. Volcanic ash hazards to aviation

During explosive eruptions, large quantities of volcanic material (gas, ash, aerosols) are injected into the atmosphere as an ash plume, reaching tropospheric and sometimes even stratospheric heights. The subsequent ash cloud drifts according to the dominant winds in altitude, generating hazards for aviation, human health, and the environment (pollution) both in the short and long terms. Whereas lapilli and coarse ash (> 1 mm in diameter) fall out within an hour, fine ash (< 1 mm, atmospheric residence > 30 min) and especially very fine ash (< 30 μm, atmospheric residence from three hours to several days) are agents for volcanic cloud hazards to aircraft and aviation services [1]. Impacts can be engine power loss, blockage of sensors, impaired vision owing to windscreen damage, contamination of aircraft ventilation and pressurisation systems, health risks to crew and passengers, abrasive damage of aircraft exterior, malfunction of electric/electronic elements with potential short-circuits, manoeuvre for cloud avoidance, flight disruptions, deposits of volcanic ash on runways, limited ground operations, airspace and airport closure, and all subsequent implications for travelling passengers [2]. Because very fine ash is often dispersed thousands of kilometres from the source volcano, it may remain for weeks in the atmosphere in amounts hardly detectable by aircraft, Air Traffic Control and weather radars, or even visually. Ash products from satellite imagery used for real-time monitoring (e.g., Jiménez-Escalona et al. [3]) sometimes also fail to detect ash due to (i) overcast weather, (ii) mask effects from very fine ash acting as cloud condensation and ice nuclei resulting in a mixed plume/cloud containing liquid water, SO₂ and ice, as reported by Marchese et al. [4] for the 3 December 2015 paroxysm of Etna, and (iii) dilute particle concentrations below detection threshold in distal areas. Even when ash is detected, communication delays to aircraft from the hazard occurrence time remain another limitation. In our study, we model the dispersion of very fine ash with particles less than 10 μm in size (PM10), as these are also well-known to affect air quality and human health (sub-10 μm inhalable volcanic ash, including the subset of respirable volcanic ash below 4 μm) (e.g., Eychenne et al. [5]). Very fine volcanic ash, SO2 and aerosols further impact the environment and climate through radiative forcing (e.g., Sellito et al. [6]).

The modelling of volcanic plume transport has attracted increasing attention to the natural hazard research community because of the impact of volcanic emissions dispersal on aviation and the economy. The global economic impact of the eruption of Eyjafjallajökull in 2010 has been estimated by Oxford Economics to be ~US$ 5 billion, with about US$ 2 billions of direct financial loss to the aviation industry. According to EUROCONTROL [7], it could be more than halved owing to improvements achieved in the air sector regarding crisis management and optimised operational performances in about a decade. For example, cost-based trajectory optimization based on dispersion models and flight planning software decreases cost and disruption of air traffic deriving from volcanic events [8].

The number of reported aircraft encounters with volcanic ash clouds causing various degrees of airframe or engine damage amounted to 129 between 1953 and 2009 [9] and to 113 between 2010 and 2014 [10]; yet, not all encounters are reported and therefore their number remains underestimated. While the most damaging encounters have occurred within 24 h of eruption onset and/or within 1000 km of the source, less safety-significant but still economically damaging encounters have occurred at greater distances and extended times [9]. In parallel, the total number of passengers transported each year also increases, including by a factor of 2.8 in the two decades preceding 2019 when it reached over 4.5 billion [2].

As demonstrated in this paper, Mount Etna emissions present significant hazards to regional air traffic, in addition to significant economic losses. The closest international aeronautic infrastructure, the Catania airport, must frequently limit its services because of Etna ash emissions, leading to numerous rerouted or cancelled flights and airport closures. About 11 million annual commercial passengers have transited through Catania airport in 2022 (6.4 and 7.1 million in 2013 and 2015, respectively), for about 76000 flights in 2022 involving 86 routes to domestic and European destinations (https://cataniaairport.com/statistics/, accessed on 3 May 2023).

1.2. Volcanic Ash Advisories for Etna

With the purpose of keeping aviators informed about volcanic hazards, the International Civil Aviation Organization (ICAO) established nine Volcanic Ash Advisory Centres (VAAC) covering the world as part of the International Airways Volcano Watch (IAVW). Each of the nine centres under the IAVW responds to reports of volcanic ash within their region and provides the aviation community with information on the ash cloud's current location and forecasts of future extent in case of significant volcanic eruptions. Ash warnings issued by the VAACs are in the form of standardised Volcanic Ash Advisory messages (VAA, sometimes including Graphics: VAG) based on Volcano Observatory Notices for Aviation (VONA), aircraft in flight (AIREP) and governmental agencies utilising meteorological ground stations and satellite imagery (SIGMET, NOTAM, ASHTAM) (https://www.ssd.noaa.gov/VAAC/vaac.html).

The Toulouse VAAC of Météo-France (https://vaac.meteo.fr/), with London VAAC as a backup, is in charge of the nowcasting of volcanic emissions dispersal in a very large zone of the globe including the Mediterranean basin hosting very active volcanoes like Etna and Stromboli. Between 2006 and 2022, the Toulouse VAAC has released 792 Volcanic Ash Advisories (VAA, most including Graphics - VAG) for Etna, including 22 and 29 during/following the eruptions from the two sequences analysed in this paper (Tables 1, 2, A1). The reaction time of the VAACs depends mainly on the on-site volcano observatory that must timely provide a VONA with critical information about the on-going activity and input parameters. At Etna, the all-year 24/7 monitoring is mainly carried out by INGV-Osservatorio Etneo, transmitting VONA, but also by information to the Italian Department of Civil Protection (https://www.protezionecivile.gov.it/en/) through specific communication procedures during crises.

As for volcanic hazards to aviation, present challenges include the detection of volcanic ash clouds, the accurate forecasting of their dispersion, and the timely and targeted communication of this information. Dispersion model simulations may contribute to evaluate the space-time distribution of volcanic constituents (solid particles, gases, aerosols), including for ash dispersal nowcasting as implemented by the VAACs during eruptions. Understanding engine and airframe tolerances to ash ingestion and gas effects would further inform the operational risk management of airlines.

1.3. Improving ash dispersal modelling: coupling WRF-Chem and time-varying Eruption Source Parameters

Despite the knowledge gained with remote sensing data and the numerical modelling of volcanic eruptions, the identification of hazard areas around active volcanoes still remains a considerable scientific challenge [11]. The parameters that most influence the forecast of ash column loading distributions for mapping flight hazard areas are the free-troposphere turbulence levels, precipitation threshold for wet deposition, and eruption source parameters such as plume height and mass eruption rate (MER), and the onset and end times of the paroxysm [12,13,14]. Egan et al. [15] have implemented a simplified volcanic ash aggregation scheme into the Weather Research Forecasting model with Chemistry (WRF-Chem; Grell et al. [16]) for computing the aggregation rate coupled to the atmospheric environment at each model step. Harvey et al. [13] have listed the observation of plume height through mobile radars among the top research priorities for reducing uncertainties in forecasts of long-range transport of volcanic ash. WRF-Chem model has been recently utilised [17] to model the volcanic plumes associated with Etna’s paroxysmal episodes on 3–7 December 2015 that led to the closure of the nearby Catania International Airport. It emerged that the WRF–Chem model reasonably reproduced the distribution of SO₂ and volcanic ash, as compared with multi-platform data from satellite sensors and Doppler radars. More recently, Rizza et al. [18] have demonstrated the need for time-dependent Eruption Source Parameters (ESP) together with the correct specification of injection height to properly describe the transport of Etna volcanic ash in the Mediterranean basin. The online coupling between chemistry/aerosols with meteorology represents an innovative development of current volcanic ash transport and dispersion models (VATDM) as it avoids any space-time interpolation by considering the two-way feedback for all the processes (meteorology, chemistry, and physics) at the time step level, in a context of time-varying ESP [18].

The main aim of this paper is to document the large-scale contamination of airspace by a succession of closely interspersed explosive eruptions, providing insight into hazards to air traffic and atmospheric dispersion of volcanic pollutants. It also demonstrates the usefulness and high potential for nowcasting of ground-based near-vent radar observations to provide the most realistic source term estimates of volcanic ash that will undergo long-range transport. To these purposes, we use MER retrievals from a L-band Doppler radar (carrier frequency 1.274 GHz) monitoring Etna (VOLDORAD-2B, here noted V2B; [19, 20]) in input to WRF-Chem atmospheric model simulations to evaluate the regional-scale dispersion and accumulation of very fine volcanic ash released by successive paroxysms. This full coupling methodology, using measured time-evolutive source parameters, accurately reproduces the atmospheric dispersal from volcanic transient emissions, as attested by satellite imagery from various gas and ash sensors for the 23 November 2013 paroxysm [18]. The higher time resolution in the definition of the meteorological conditions and the specification of eruption conditions allows the reduction of related uncertainties present in the modelling of transport and dispersal of tephra in the atmosphere.

In the present study, we analyse two test-case eruptive sequences of Etna, the first one including 6 successive paroxysms from the New SouthEast crater (NSE) that occurred in 38 days between October 26 and December 2, 2013, the second sequence including 4 paroxysms from the Voragine crater (VOR) that occurred in less than 3 days between December 3 and 5, 2015. The methods and data used in this study are described in Section 2, which also describes the details of the two eruptive sequences, the volcanic emission estimates determined using the V2B radar data and the model setup. In section 3, we discuss the modelling results, in terms of meteorological and ash transport characteristics in the numerical model domain. Conclusions and aims for future works are presented in Section 4.

2. Materials and Methods

2.1. Description of the two eruptive sequences

2.1.1. The 2013 sequence: October 26 to December 2

The first simulation reproduces the sequence of paroxysms that occurred at the New South East (NSE) crater from 26 October to 2 December 2013. After a 6-month quiescence, a new paroxysm began on 26 October at NSE crater. For this event, the SO2 mass detected by the OMI (Ozone Monitoring Instrument) sensor onboard Aura satellite peaked at 0.77 kilotons (https://so2.gsfc.nasa.gov/index.html, accessed on 05 May 2023). The October-December sequence includes six paroxysms featuring usual lava fountain feeding many kilometres high tephra plumes, ranging in duration from 3.2 to 10 hours and from 2.9x10⁴ to 3.9x10⁵ kg s^-1 in average MER (peaks: 2.1x10⁵ to 8.9x10⁶ kg s^-1) as retrieved from VOLDORAD-2B radar (Table 1). Only the last three paroxysms punctually exceeded 10⁶ kg s^-1, the shortest and most violent of the series occurred under strong westerly wind on November 23 [18, 21, 22]. This paroxysm averaged nearly 4x10⁵ kg s^-1 during 193 min including 3.5x10⁶ kg s^-1 sustained during its 16 min climax phase which released 80% of the total erupted mass of pyroclasts at maximum velocities over 200 m s^-1 [20]. The six successive paroxysms released in total about 18.9 x 10⁹ kg of tephra with respective mass contributions of 18, 40, 24, 5, 8, and 5 %.

Plume top altitudes for this sequence (8.1–11.1 km a.s.l.) can be found in Corradini et al. [23]. The Toulouse VAAC Météo-France released 23 VAA for Etna during the period of the six paroxysms for flight levels up to FL380 (11.58 km a.s.l.) with wind speed estimates up to 35 m s^-1.

2.1.2. The 2015 sequence: December 3 to 5

An exceptionally rapid sequence of eruptive episodes occurred in December 2015 generating powerful eruptions columns during paroxysms, tall lava fountains, Strombolian activity, lava flow emission and persistent ash emission, involving all four summit craters at different times [24, 25]. This activity culminated with four paroxysmal events on December 3–5, 2015 from the Voragine (VOR) crater constituting the second sequence of our analysis (Table 2, Figure 2). As most paroxysms at Etna, they were preceded by Strombolian activity and displayed lava fountains and subsequent high tephra plumes. These paroxysms were short (less than 1.5 h) and not accompanied by lava flows and, as compared with the 2013 sequence, they were characterised by exceptional eruptive power and rapid occurrence rate. They rank among the most violent in the last two decades, in part due to the depth of their magma source region 1.5 km b.s.l. [26,27].

The first episode on December 3 was the most powerful, with lava fountain velocities inferred from V2B averaging 200 m s^-1 during the climax and peaking at 250–290 m s^-1. An ash plume ascent velocity of 5.6 m s^-1 (10 km in height above crater level in about 30 min) is inferred from the timing of the V2B near-source power and velocity increase and the plume top altitudes from MSG-SEVIRI HOTVOLC data and vertical temperature profiles [28], the maximum plume top altitude of 13.5 ± 1.5 km a.s.l. being reached around 03:00 UTC. The following paroxysms gradually decreased in MERs, total mass of erupted pyroclasts, and lava fountain average height [26].

These paroxysms caused the atmospheric dispersion and widespread ash fallout to more than 100 km away, with disruptive consequences to aviation and aeronautical infrastructures and services. In particular, the Catania international airport reported repeated closures of airspace sectors between 3 and 9 December 2015, and at least 57 cancelled flights, 56 diverted flights and tens of modified flights on 4–5 December. The Toulouse VAAC Météo-France released 29 VAA in 7 days (December 3–9) for flight levels up to FL350 (10660 m a.s.l.) and wind speeds up to 15 m s^-1, while INGV released 16 VONA in four days (December 2–5) reporting on these paroxysms.

The regional environmental impact of volcanic emissions from this short series of paroxysms can be apprehended from estimates of released products amounts. According to Bonaccorso and Calvari [26], the total production of pyroclasts was 10x10⁶ m³ DRE (Dense Rock Equivalent, i.e. about 25x10⁹ kg) with respective mass contributions of 41, 27, 19, and 13 % for the successive events. The V2B estimates point to 11x10⁹ kg with successive mass contributions of 84.1, 6.5, 4.8, and 4.6 %. Differences among estimates can stem from both the numerous hypotheses used by the former to derive magma volumes from strainmeter and camera measurements and from the possible underestimation by V2B of the last three weaker paroxysms due to the marginal position of the Voragine crater with respect to the radar beam. The SO₂ mass detected by the OMI sensor peaked at 20 kilotons in early December (https://so2.gsfc.nasa.gov/index.html, accessed on 16 May 2023).

The definition of eruption source parameters (ESP), namely the volcano location in the grid domain, start and end of each single paroxysm (UTC time), the total grain size distribution (TGSD), MER, and the injection height was pre-processed as described in the following paragraphs.

2.2. The VOLDORAD-2B (V2B) Doppler Radar System

2.2.1. Description

Between 2009 and 2022, INGV-OE monitored Etna tephra emissions using a volcano Doppler radar of OPGC (VOLDORAD-2B), jointly operated from the Montagnola station about 3 km south of the New SE Crater. This L-band (23.5 cm-wavelength) pulse radar probed 13 atmospheric volumes, 150 m-deep, right above the summit craters using a fixed-pointing antenna beam. The echo power and Doppler velocity parameters (level 1 data) in each range bin are computed at intervals of 0.22 s from 3 incoherent integrations of the Doppler spectra (raw data). These near-source measurements of emitted tephra are unique as they allow maximum ejection velocities at the base of the eruptive column and MER (level 2 products) to be determined, potentially in real-time [19,20] and to a satisfying level of uncertainty as revealed by comparisons with other sensors [29,30,31]. During the 2011–2015 period, the V2B radar successfully recorded about 50 paroxysms of vigorous lava fountains feeding tephra plumes up to over 13 km high. Level 1 data are available from an open-access database of the Observatoire de Physique du Globe de Clermont-Ferrand (http://voldorad.opgc.fr/bddtr.php, accessed on 16 June 2023; [32]), as well as more elaborated (level 2) products upon request, such as MER, onset and end of fountaining phase, and eruptive crater.

2.2.2. Data elaboration

We used the methodology of Freret-Lorgeril et al. [20] to compute the MER from the V2B data. They showed that the product of the measured echo power and radial velocities of rising tephra in the beam volumes located immediately above the eruptive crater is proportional to the MER (i.e., a MER proxy) and correlates to the observed plume heights for many paroxysms at Etna. Their calibration constant was then found by comparing the radar proxies with MER derived from a plume ascent model matching observed plume heights and including wind vertical profiles [33].

In this study, the MER was computed from range bins displaying the maximum power, at 3135 and 3285 m in slant distance from radar for paroxysms originating from the New SE Crater during the first sequence (26 October–2 December 2013, Table 1) and from 3 range gates for the second sequence of paroxysms originated from the Voragine Crater (Table 2). The V2B MER data were elaborated to be compatible with the WRF-Chem input system; in particular, they were averaged on a 30-second interval to be integrated in the numerical model at each time step. The corresponding time series of the MER are depicted in Figure 1 and Figure 2 for the 2013 and 2015 eruptive sequences respectively.

Figure 1. Mass Eruption Rate from V2B radar system ingested in WRF-Chem volcanic package for the 2013 sequence: (a) full sequence, (b) NSE2 paroxysm (c) NSE3 paroxysm (d) NSE4 paroxysm (e)NSE5 paroxysm, and (f) NSE6 paroxysm; Julian day 290 corresponds to October 17, 2013.

Figure 1. Mass Eruption Rate from V2B radar system ingested in WRF-Chem volcanic package for the 2013 sequence: (a) full sequence, (b) NSE2 paroxysm (c) NSE3 paroxysm (d) NSE4 paroxysm (e)NSE5 paroxysm, and (f) NSE6 paroxysm; Julian day 290 corresponds to October 17, 2013.

Figure 2. Mass Eruption Rate from V2B radar system ingested in WRF-Chem volcanic package for the 2015 sequence of paroxysms from Voragine crater: (a) full sequence; (b) VOR1; (c) VOR2 and VOR3; (d) VOR4. Julian day 335 corresponds to December 2, 2015.

Figure 2. Mass Eruption Rate from V2B radar system ingested in WRF-Chem volcanic package for the 2015 sequence of paroxysms from Voragine crater: (a) full sequence; (b) VOR1; (c) VOR2 and VOR3; (d) VOR4. Julian day 335 corresponds to December 2, 2015.

2.3. Determination of the eruption height

In case of missing experimental data for the model injection height, the zero-order semi-empirical relationship of Mastin et al. [34] linking column height with MER observations can be utilised. It is expressed by the following eq.1 [35], and reported in Figure 3:

$$h=5\times {10}^{-4}{\left(\frac{MER}{\rho }\right)}^{0.241}$$

(1)

where h is expressed in metres (above the vent), ρ denotes the magma density (2600 kg m^-3) and MER is obtained by the above-described elaborations of the V2B radar data and reported in Figure 1 and Figure 2.

For the range of mass eruption rates considered in the two series of eruptions [0–6] x 10⁶ kg s^-1, the corresponding eruption height has an upper limit of 12–13 km (above the vent), in agreement with Mastin et al. [34] (their Figure 1).

The emission profile is specified according to the vertical profile derived from the umbrella cloud model in which 75% of the erupted mass is constrained in the umbrella cloud and the remaining 25% is linearly distributed from the umbrella base down to the vent [35].

2.4. Model description and setup

For this study, the WRF-Chem model version 4.3.1 has been utilised in a numerical domain covering the central Mediterranean region, with 310 × 320 grid points and a horizontal grid spacing of 6 km (1854 × 1914 km², Table 3). Boundary and initial conditions are at 1-degree resolution and provided from NCAR/NCEP Final Analysis from Global Forecast System (FNL from GFS). The first sequence simulation started at 00:00 UTC on October 25 and ended at 00:00 UTC on December 04, 2013. The second sequence run started at 00:00 UTC on December 2 and ended at 00:00 UTC on December 08, 2015.

The full set of parameterisations of the WRF-Chem model are reported in Table 3.

Based on the WRF setup recommended by Shi et al. [36] and Rizza et al. [37]the following physics schemes are utilised: the Mellor–Yamada–Janjic (MYJ; [38]) and the Eta similarity theory [39] schemes are used to describe the planetary boundary-layer (bl_pbl_physics = 2) and the surface layer (sf_sfclay_physics = 2) processes, respectively. The land surface exchange coefficients and the associated fluxes of momentum, heat, and humidity are represented by the Noah-MP Land Surface Model (sf_surface_physics = 4; [40]).

The radiative schemes are parameterized using the Goddard radiation model [41] for shortwave (ra_sw_physics = 5) and longwave (ra_lw_physics = 5) components. The microphysics parameterization considers the new one-moment Goddard four-class ice (4ICE) scheme developed by Lang et al. [42] and implemented in the early versions (3.4.1) of the WRF package (mp_physics = 7). This scheme has prognostic variables for cloud ice, snow, graupel, and hail.

In WRF-Chem, the total grain size distribution (TGSD) consists of ten bins of volcanic ash (vash_#) particles with a grain diameter range starting from 2 mm down to less than 3.9 μm [35] and expressed in mixing ratio units ($\mu$g kg^-1). The corresponding mass fraction percentage is described by the E1 distribution and reported in Table 4. It has been recently introduced inside the WRF-Chem volcanic package by Rizza et al. [18] and derived from the Poret et al. [22] analysis of the experimental data from the November 23, 2013 Etna eruption.

3. Results and discussion: hazard maps

3.1. Atmospheric circulation during the events: main mid-tropospheric flows

In this section we describe the general synoptic situation in terms of prevailing mid-troposphere flows, i.e., at 500 hPa, which also represents an important flight level (see next section). In Figure 4 and Figure 5 are reported the geopotential height and the wind vectors derived by WRF-Chem simulations. Data are shown in correspondence of the eruption timeline for both the 2013 and the 2015 sequences; for each event, the chosen times are close to the paroxysms and are reported in Table 5.

Although studied as a unique sequence, the 2013 eruptions lasted 38 days and the atmospheric configurations continually changed; it nevertheless generated some recurrent flows. For some of the six 2013 eruptive events several similarities are in fact noticeable in terms of synoptic configuration and prevailing upper-level winds. For NSE2, NSE4 e NSE5 (Figure 4b, 4d, 4e, respectively) some upper-level troughs, although of different amplitude and depth, dominated the central Tyrrhenian favouring south-westerly winds on southern Italy and on the Sicily region in particular.

For NSE6 (Figure 4f) a trough is also visible, but located south of the Strait of Sicily, between the southern coasts of the region and north Africa. This configuration favoured the advection of south-easterly upper-level winds over eastern Sicily.

An upper-level ridge is visible on southern Mediterranean for the NSE3 case (Figure 4c); this structure is connected to a trough on the Balearics, and the main mid-tropospheric flows over north-eastern Sicily followed the curvature of the ridge, i.e., with currents toward east-southeast.

Finally, a different (from the previous ones) configuration is noticeable for NSE1 (Figure 4a), characterised by a wide high-pressure area over the whole Italian Peninsula. In this situation, the upper-level currents over Sicily were weaker and mainly directed towards the southwest, also driven by a lower pressure area over North Africa.

For the 2015 sequence (Figure 5), that lasted about 3 days, the main synoptic configuration remained quite stationary and has been characterised by a very wide high-pressure system persistently dominating the whole Mediterranean Basin, only slowly moving westward during the sequence. Some local (and weak) high- and low-pressure areas, visible on Figure 5, developed over the Mediterranean during this sequence, maintaining relatively weak mid-tropospheric currents; these winds flowed over north-eastern Sicily firstly toward east (Figure 5a, 5b) and later, after the formation of a local low over north Africa, toward north (Figure 5a, 5b).

To highlight some general differences between the two sequences it can be seen how the 2013 eruptive events were characterised by a succession of different baric configurations, most of the time characterised by the presence of low-pressure areas over the central Mediterranean. On the contrary the 2015 events, that occurred in a few days, were characterised by a prevailing and quite stationary wide high-pressure system, with smaller pressure gradients in the Mediterranean. The resulting circulations at 500 hPa were also quite different, with mid-tropospheric currents on average in the range 20-40 m s^-1 for the 2013 (Figure A2) and less than 10 m s^-1 for the 2015 (Figure A3) sequence, respectively. These results, based on WRF-Chem simulations, represent a robust meteorological basis for the evaluation of the regional-scale ash transport after volcanic eruptions (see next section).

3.2. Analysis of the volcanic ash dispersion patterns

The objective of this work is to evaluate the contamination of airspace by a succession of Mount Etna explosive paroxysms and the implications for hazards, in particular with respect to air traffic in the Mediterranean area. For this purpose, we define the vertically integrated ($\overline{p_{10}}$) mass, that is calculated considering the ash bins of the E1 distribution with diameter less than 10 μm (Table 4), as follows:

$$\overline{p_{10}}={\int }_0^{top}<PM10>{\rho }_{air}dz$$

(2)

$$PM10=0.5\times \mathrm{v}\mathrm{a}\mathrm{s}\mathrm{h}\_10+\mathrm{v}\mathrm{a}\mathrm{s}\mathrm{h}\_9+0.5\times \mathrm{v}\mathrm{a}\mathrm{s}\mathrm{h}\_8$$

(3)

where PM10 is the sum of the finest ash bins (units μg kg^-1) and ${\rho }_{air}$ is the density of the air (kg m^-3). Brackets denote the time average over the simulation period and the vertical integration (represented by the overbar) is performed considering the whole vertical domain (top = 18 km). The resulting $\overline{p_{10}}$ units are expressed in μg m^-2 which can be considered as a columnar ash concentration.

The model simulations produce hazard maps of $\overline{p_{10}}$ to inform on the atmospheric ash dispersal of PM10 but it should be reminded that time- and spatially-averaged representations presented below are intended to guide aviation in terms of hazard appraisal and risk mitigation of Etna paroxysms. These representations should not be confused with instantaneous ash concentrations that could be met by aircrafts or measured by instruments at a given moment. To circumvent the lack of vertical information in columnar ash concentration representations, we derive ash concentrations at specific flight levels in Section 3.2.2 by performing a WRF-Chem data interpolation (3D) for a given pressure level (specialised function of NCAR Command Language named wrf_user_intrp3d). To further image the ash dispersal with a better time resolution, in Section 3.2.3 we also compute the daily average ash concentration at specific flight levels for four selected paroxysms.

3.2.1. Time-averaged columnar ash concentration maps

Figure 6 and Figure 7 show the results for the 2013 and the 2015 sequences, respectively. Both figures were built considering the time average on the full simulated sequences. For the 2013 sequence, a successive spatial filter is applied to leave out days without ash in the domain.

These two figures give a quantitative estimation of the global contamination of airspace by a succession of eruptions. The average spatial pattern outlined by the 2013 sequence (Figure 6) shows basically three affected regions in terms of vertically-integrated columnar ash concentrations: a red-shaded region around 10 g m^-2 within a radius of approximately 100 km, a yellow-shaded region of the order of 0.5 g m^-2 within a ~ 500 km radius, and a green/blue-shaded area with columnar ash content lower than 0.1 g m^-2 within a ~ 800 km radius. This average spatial pattern is overprinted by the dispersal from each event oriented by its specific wind conditions, conferring the starfish-like overall pattern.

The distribution of columnar ash concentrations from the 2015 short sequence (Figure 7) presents a similar starfish-like spatial pattern, with remarkable differences consisting in: (i) a more extended red-shaded region centred on Sicily and south Italy and also more intense in term of magnitude, within a radius of approximately 150 km from the eruptive crater and (ii) a strong eastward upper tropospheric transport caused by the VOR1 eruption.

As highlighted in the previous section, these differences are mainly related to the different mid-tropospheric transport conditions during the two sequences and also by the very rapid succession of eruptive episodes in December 2015, that generated powerful eruptions columns up to the stratosphere.

This analysis may be considered a preliminary indication of potential hazard zones around Mount Etna, which can be re-evaluated regularly as meteorological conditions change.

3.2.2. Time-averaged ash concentration maps at specified flight levels

To evaluate the contamination of airspace and the related hazards to aviation, it is necessary to point out that the forecasts by VAACs (https://www.ssd.noaa.gov/VAAC/vaac.html) use the following concentration thresholds to define the ensuing levels of ash contamination: low (200–2000 μg m⁻³), medium (2000–4000 μg m⁻³), and high (> 4000 μg m⁻³).

The flight level (FL), represents the altitude of an aircraft at the standard air pressure, expressed in hundreds of feet (1 FL = 100 ft). It may be considered as a proxy for elevation, and it may be associated with a given pressure level as reported in Table 6. These levels correspond to the upper airspace characterised by ICAO, as the zone where most air transport operations take place.

Figure 8 displays the time-averaged PM10 concentration for the 2013 (upper row) and 2015 (bottom row) sequences, at flight levels FL180 (panels a,d), FL240 (panels b,e) and FL300 (panels c,f) corresponding to the pressure levels of 500, 400 and 300 hPa (Table 6). The label-scale of Figure 8 is saturated at 4000 μg m^-3 (corresponding to 10^3.6) which represents the threshold defined by ICAO for the high level of contamination and the consequent prohibited airspace.

The transport of volcanic ash at these flight levels has different patterns for each respective sequence. The effects of the low-pressure weather conditions that occurred during the 2013 sequence is particularly evident in Figures 8a-c from the circular to elliptical stretched vortices (800 to over 1600 km in diameter) that transported volcanic ash counter-clockwise at large distances. This is particularly evident for the 2013 sequence for FL180 (Figure 8a) and FL240 (Figure 8b) flight levels at which ash remains concentrated mainly along these vortices in narrow bands generally less than 100 km wide. Ash concentrations in these bands range from about 1 to over 4000 μg m^-3. Being well localised under these high wind conditions, ash can remain concentrated over long distances, up to 25 μg m^-3 (10^1.4) still at 1000 km eastward for the most violent paroxysm of the sequence that occurred on 23 November 2013 (NSE4). At FL300 (300 hPa) in the upper troposphere, vortices are still present, but ash is more widely dispersed in the north direction owing to tropospheric wind at 300 hPa (Figure 8c).

Contrastingly, the 2015 sequence, occurring in prevalently high-pressure conditions and weak winds (see Figure A3), caused high ash concentration at all flight levels, as outlined by the red shaded area in Figures 8d-f. At FL300 (Figure 8f), mostly the intense eastward transport of very fine ash from VOR1 by tropospheric wind at 300 hPa is visible, due to the combination of the eruption intensity and high-pressure low wind conditions. Considering this is a time-averaged map, this potentially indicates high ash hazards at all flight levels during - and for days following - this eruptive sequence. The concentration of very fine ash appears more localised spatially going upward into the upper troposphere, although still reaching no fly zone level.

What can be inferred from this analysis on the PM10 maps at these flight levels is that the long-range transport is strictly correlated with the synoptic circulation. In particular, synoptic patterns with relatively weak tropospheric currents are conditions leading to the accumulation of volcanic ash in the atmosphere after an eruption. This seems a rather obvious conclusion, but it is important to highlight the fundamental role played by the meteorological component in VATDM models and the need for their online coupling with the numerical aerosol packages.

This result is in part confirmed by analysing the VAA issued by Toulouse-VAAC that are reported in Table A1: a total of 19 VAA were issued in the period December 3 to 6, of which 12 were classified as “Red”, 6 “Orange”, and 1 “Yellow”.

3.2.3. Daily average ash concentration maps of selected paroxysms

In the following Figure 9 and Figure 10 we have analysed the daily average (from 00:00 UTC until 23:59) ash concentrations of two powerful paroxysms of the first sequence having the most contrasted durations on November 23 (NSE4, 3.2h) and 28 (NSE5, 8.3h), and the most contrasted in MER for the second sequence, on December 3 (VOR1) and 5 (VOR3-4). It is calculated using the following quantity:

$$p_{10D}=\left[\mathrm{P}\mathrm{M}10\right]$$

(4)

where square brackets denote daily time average on the selected day. The left column of both figures (panels a-d) refers to FL180, the central to FL240 (panels b-e) and the right column to FL300 (panels c-f). The contour shading matches the above-mentioned four hazard levels used by aviation for volcanic ash concentration reported in Table 7.

For the NSE4 paroxysm, the geopotential height at 500 hPa (Figure 4d) shows a wide low-pressure area on the central Mediterranean basin favouring south-westerly strong upper-tropospheric winds (> 30 ms^-1) on southern Italy. The very fine ash was rapidly transported toward Albania (in about 6h) and the Balkan regions at relatively low concentrations (Figure 9a).

The very fine ash after the NSE5 paroxysm (Figure 9d,e,f,) was transported eastward in a straight line by the intense zonal currents at 500 hPa (Figure 4e).

These maps indicate low potential risks for aviation at FL180 to FL300. The shading contours (blue-green) indicate in fact the low hazard level. On the other side all the 21 VAA by Toulouse VAAC reported in Table A1 indicate “NIL” alert level (NIL = no ash cloud is produced) for the whole period.

The corresponding analysis of the 2015 sequence is depicted in Figure 10. Panels a-b-c in the upper row report the daily p_10D concentration at FL180, FL240, and FL300 for December 3, 2015 (VOR1), while the lower panels (d-e-f) describe the equivalent for December 5, 2015 (VOR3, and VOR4).

For both cases, the p_10D daily patterns, in agreement with the high pressure meteorological conditions (Figure 5a,d), indicate a red-shaded region in proximity of Mount Etna for FL180 and FL240 (panels a-b-c-d). FL300 (Figure 10c) relative to VOR1 again reveals the presence of p_10D ash fraction at upper tropospheric level (about 9200 m) in the form of wave-like pattern which is the artificial effect of the averaging (daily) of pulsed ash plumes. To better visualise these maps, a high-resolution Figure A1 is provided in the appendix with a zoom on the Mount Etna area. This is consistent with the 19 VAA issued by the Toulouse VAAC during the period 3 to 6 December 2015, 12 of which were classified as “Red” and 6 as “Orange”.

3.2.4. Vertical distribution of very fine ash

Animated videos (gif files) showing the fast volcanic PM10 space-time evolution of the successive emissions (often 'pulse-like') due to the transient nature of their source can be found at https://doi.org/10.25519/GG1K-2H34 for the 2013 and 2015 sequences. Figure 11 are extracted snapshots for different paroxysms informing on the vertical distribution of PM10 along cross-sections indicated in Figure 6 and Figure 7. The ash cloud thickness is a particularly valuable information for aviation hazard assessment and countermeasures in case of ash encounter, and for satellite mass retrievals as generally a non-measured average thickness is assumed. For instance, the ash cloud inherited from the 26 Oct. 2013 tephra plume (NSE1) is seen slowly drifting to the SW after about 1.5 day between 1.9 and 3.7 km a.s.l., i.e., with an inferred vertical thickness of about 1.8 km (Figure 11a). After the top part of the VOR1 plume had raised above 12 km a.s.l. in the minutes following the eruption onset, some very fine ash from the middle region of the eruptive column is observed after 7 h rising to the NE above 11.8 km a.s.l. with an inferred vertical thickness of about 2.5–3 km (Figure 11c), while the lower part of the eruptive column feeds another PM10 cloud, 700 m thick, drifting at about 2 km a.s.l. laterally to the SW 1.5 day after the onset (Figure 11d). The modelled ash cloud from the 4 Dec 2015 evening paroxysm has a vertical thickness of 1.4 km at an altitude of 5 km a.s.l. when the following paroxysm occurs (Figure 11f). Note that an ongoing vertical plume emission coexisting with ash from a previous plume is visible in Figure 11e–f.

Another interesting feature seen in Figure 11b,e is ash stagnating for several days at ground level. PM10 concentrations are commonly above 500 μg m^-3, as during four days in late October to early November 2013 following NSE1. PM10 concentrations at ground level can even punctually exceed the aviation threshold as observed following VOR2 around the occurrence time of VOR3. The repeated ash emissions of Etna recurrently exposing surrounding inhabitants, animals, and plants to volcanic PM10 and gas may therefore potentially represent health hazards over the long term.

4. Conclusions

We evaluate for the first time the contamination of airspace by very fine ash from a succession of ash plumes produced by Mount Etna and the related hazards to aviation in the Mediterranean Basin. To this aim, the atmospheric dispersal of sub-10 μm (PM10) ash was simulated using the WRF-Chem model coupled online with meteorology and aerosols and offline with Mass Eruption Rates (MER) derived from near-vent Doppler radar measurements and inferred plume altitudes. The pioneering direct input of systematic observations, best resolved in time and thus improving the model accuracy, is all the more so valuable as paroxysms are transients involving time-varying eruption source parameters. PM10 ash concentration hazard maps are provided for two sequences of paroxysms with widely varied volcanological conditions and contrasted occurrence interval and meteorological synoptic patterns: 6 ash plumes produced in 38 days during low-pressure, high-wind conditions in October - December 2013, and 4 ash plumes released in 3 days during high-pressure, low-wind conditions on 3-5 December 2015. Average MERs covered two orders of magnitude (2.9 x 10⁴ – 1.7 x 10⁶ kg s^-1), durations spanned from less than 1 hour to nearly half a day, and plume top altitudes ranged from 8.1 to 13.5 km a.s.l. The synoptic conditions at 500 hPa were quite different, with mid-tropospheric current speeds on average in the range 20–40 ms^-1 and less than 10 ms^-1 for the 2013 and 2015 sequences, respectively.

We analyse the ash PM10 dispersal in terms of time-averaged columnar ash concentration ($\overline{p_{10}}$), of concentration at specified flight levels (FL180, FL240, and FL300) averaged over the entire sequence interval (PM10), and of daily average concentration (p_10D) during selected paroxysm days at these flight levels. Consistently with the Volcanic Ash Advisories issued by the Toulouse VAAC, they may be considered as first-approach model hazard maps to aircrafts and, although certainly not covering all possible volcanological and meteorological configurations for hazard appraisal at Etna, the following conclusions can be drawn:

• The very fine ash from sequences of Mount Etna paroxysms is shown to easily contaminate the airspace around the volcano within a radius of about 1000 km in a matter of days. The airspace from many countries around the Mediterranean Basin is impacted, including most of southern Europe from the Balearic Islands westward, south of France, the whole Italy, Greece and western coast of Turkey and the Balkans eastward, to beyond the Alps northward and to Malta and the African northern coast (from Algeria to Libya) southward.
• Low-pressure weather systems favour the trapping and circulation of very fine ash in the whole troposphere within this area, yet at low concentration, generally below 1 μg m^-3. In this meteorological context, high winds tend to stretch ash clouds into ~100 km-wide clouds forming large-scale vortices 800–1600 km in diameter where PM10 ash concentrations can still exceed the aviation hazard threshold up to 1000 km downwind from the volcano, a distance reached in about 10h (e.g., NSE5).
• High-pressure, low-wind conditions tend to favour the accumulation of PM10 ash in a wide atmospheric region surrounding Etna. In this context, closely interspersed paroxysms tend to accumulate very fine ash, more diffusively in the lower troposphere and in stretched ash clouds higher up in the troposphere.
• In all volcanological and meteorological configurations simulated, the lower troposphere appears particularly prone to accumulate diffuse PM10 ash during sequences of eruption; this is likely to particularly affect take-off and landing aircrafts in regional airports.
• High MER paroxysms propel ash up to the upper troposphere, where most of the air traffic occurs, and sometimes also into the lower stratosphere, whichever the weather conditions. High-troposphere ash clouds from Etna appear as a pulsed feature resulting mostly from the short-lived climax phase.
• Daily average PM10 tropospheric ash concentrations commonly exceed the aviation hazard threshold, up to 1000 km downwind from the volcano and up to the upper troposphere for intense paroxysms.
• The thickness of modelled PM10 ash clouds generated from different parts of the eruptive columns ranges from 0.7 to 3 km.
• Potential health hazards may stem from the stagnation for PM10 ash at ground level for several days, commonly above 500 μg m^-3, and sometimes punctually exceeding the aviation threshold.

Improving the accuracy of ash dispersal models is necessary both for nowcasting and forecasting purposes and to compensate for ash detection gaps from monitoring systems. Future WRF-Chem model developments will include refining injection heights and the initial columnar ash distribution. The semi-empirical law from Mastin et al. [34]’s dataset used in our study averages all observed ash plume heights from mainly large silicic explosive eruptions that have occurred in a variety of wind conditions, tending to overestimate height from a given MER (Freret-Lorgeril et al. [20]; their Figure 7). Other models relating MER and plume height advantageously consider the vertical wind profile (e.g., Degruyter and Bonadonna [33]). Wind profiles may come from the meteorological data online-coupled to WRF-Chem and cross-checked from atmospheric balloons measurements achieved several times a day in the region around Etna. Additional datasets for wind profiles at high resolution (30 km) could be considered, such as the ERA5 data available from the Climate Data Store in netCDF format [44]. Plume heights controlled by wind conditions specific to the eruption could thus be systematically computed from the radar-derived MER. This would lead to a more accurate definition of the injection height in the simulations.

As for future research on volcanic ash forecasts, multi-model multi-source term ensemble approaches seem to represent a promising development, as explored for instance by Plu et al. [45]. They suggest that quantiles of ash concentrations can be relevant products for air traffic management, that can be used for route optimization in the areas where ash does not pose a direct and urgent threat to aviation. Cost and disruption of air traffic could be eliminated to a great extent by including the results of dispersion models into flight planning software to apply cost-based trajectory optimizations.