Real-time probabilistic assessment of volcanic hazard for tephra dispersal and fallout at Mt. Etna: the 2021 lava fountain episodes

Starting from February 2021, Mt. Etna (Italy) experienced a period of intense explosive activity with 17 lava fountain episodes between 16 February and 1 April 2021. During the eruptive cycle, the Istituto Nazionale di Geofisica e Vulcanologia-Osservatorio Etneo (INGV-OE) issued 62 alert notifications known as VONAs (Volcano Observatory Notice for Aviation) to inform the aeronautical authorities about the volcanic activity. We present an automated VONA-based workflow aimed at real-time assessment of the volcanic hazard due to tephra fallout at Mt. Etna. When a VONA reporting tephra emission is issued by INGV-OE, numerical simulations accounting for atmospheric and eruptive uncertainties are automatically initialized to produce probabilistic hazard maps of tephra fallout and atmospheric dispersal. We applied the workflow to three lava fountains that occurred during the 2021 eruptive cycle. To test the modelling results, we compared the simulated ground load with field data, and the extent and position of the simulated volcanic cloud with the observed or estimated volcanic cloud from the Toulouse Volcanic Ash Advisory Center. Overall, we found a good match between simulated and observed quantities (tephra loads and volcanic cloud position), especially when accurate information on eruptive conditions (column height and duration) are supplied by the VONAs. Finally, through a statistical analysis, we found that column height and wind field are fundamental in determining tephra ground accumulation. For this reason, these parameters should be constrained by observational data as accurately as possible when performing numerical simulations, especially in the line of developing operational workflows.


Introduction
Among the variety of volcanic processes, explosive eruptions are of particular interest for the scientific community in the light of mitigating their impact on human lives and activities. The eruption columns formed during these episodes are hot mixtures of solid particles and volcanic gases which, entraining and heating ambient air, can easily reach stratospheric heights. The solid particles, known as tephra, cover a wide range of sizes, shapes, and densities. According to their size, such particles are classified as bombs or blocks (diameter ≥ 64 mm), lapilli (2 mm ≤ diameter < 64 mm), and ash (diameter < 2 mm). Once erupted, the largest particles decouple from the volcanic column and deposit on the ground, while the finest particles are carried upwards in the column and enter the atmosphere forming volcanic clouds. On the ground, livelihoods and infrastructures can be extensively damaged by the fallout of pyroclasts of various sizes and shapes (Barsotti et al. 2010;Wilson et al. 2014). In the air, volcanic clouds pose a serious threat to aviation safety with immediate and long-term effects on aircraft functionalities (Song et al. 2016;Giehl et al. 2017).
Volcanic hazard assessment is therefore a fundamental step in the evaluation of the risk associated with the presence of ash particles into the atmosphere and the fallout of tephra on the ground. Nowadays, numerical modelling strongly supports the assessment of various types of volcanic hazards, including the role of tephra particles emitted during explosive eruptions. Indeed, a variety of numerical codes are employed to simulate the release, transport, and deposition of volcanic particles. Some of them are as follows: CANERM (D'amours 1998), NAME (Jones et al. 2007), VOL-CALPUFF (Barsotti et al. 2008), HYSPLIT (Stein et al. 2015), MOCAGE (Guth et al. 2016), ATLAS (Reckziegel et al. 2019), and FALL3D (Folch et al. 2020;Prata et al. 2021). Some of these models were specifically developed for volcanological applications (as FALL3D), while others originally born as generic atmospheric transport models and were later adapted to simulate the transport and deposition of volcanic products (as VOL-CALPUFF and HYSPLIT). Following Folch (2012), we use the term Tephra Transport Dispersal Models (TTDMs) to point out the codes employed to simulate the release and transport of tephra particles in general, while the term Volcanic Ash Transport and Dispersal Models (VATDMs) indicates the codes considering ash particles only (diameter < 2 mm). Volcanic Ash Advisory Centers (VAACs), volcano observatories, and the many volcanological research centres around the world make use of these numerical models to produce forecasts for the dispersal and deposition patterns of volcanic particles. Extensive research is therefore done to make such forecasts reliable and suitable for volcanic hazard assessment.
In recent years, the use of TTDMs has seen a big change, from the simulation of deterministic scenarios to the production of probabilistic forecasts. This means that ensembles of simulations are performed in order to take into account the many sources of uncertainty present in the modelling of the transport and dispersal of tephra (especially for real-time applications). Indeed, uncertainties in the source term (i.e., eruptive conditions), model processes, and meteorological conditions must be considered with the aim of producing reliable forecasts suitable for operational purposes (Dacre and Harvey 2018;Selva et al. 2018;Zidikheri et al. 2018;Prata et al. 2019). Moreover, the production of probabilistic forecasts enables the assimilation of satellite data of atmospheric ash mass into the TTDMs. As shown by Pardini et al. (2020) and Mingari et al. (2022), the correction of the numerical forecasts by assimilating information on the volcanic cloud coming from satellite-based sensors can greatly improve the reliability of the forecasts.
Thus, in case of explosive eruptions, volcanic hazard can be effectively assessed and mitigated if high-quality forecasts of atmospheric ash dispersal and tephra fallout are supplied. High quality means that these forecasts should be probabilistic, done in real time, and produced with the most recent information on the volcanic activity (when available).
In this work, we present an automated workflow aimed at forecasting, in real time, both the atmospheric dispersal and the ground deposition of tephra particles released at Mt. Etna volcano (Sicily, Italy) during intense explosive eruptions. Specifically, as soon as the Istituto Nazionale di Geofisica e Vulcanologia, Osservatorio Etneo-Sezione di Catania (INGV-OE) notifies tephra emission at Mt. Etna, numerical simulations are automatically initialized to produce forecasts for the transport and dispersal patterns of the volcanic cloud and for the fallout of tephra particles in a domain covering the whole of Sicily. The forecasts are also updated as soon as a new notification with new information on the ongoing eruption, or with the indication of the end of the eruptive activity, is released. We developed the workflow while the February-April 2021 paroxysmal cycle was ongoing. Thus, we performed numerical forecasts for the majority of the eruptive episodes which occurred in that period. However, at that time, the workflow was still under development, and the results achieved were not definitive. Once the development phase was over, we selected three episodes representative of the whole cycle. These are the lava fountains which occurred on 16 February, 28 February, and 4 March 2021. For these episodes, we applied the final workflow by reproducing the timing of the issuing of the alert notifications released during the three episodes. Finally, we tested the numerical results with observations of the tephra deposit produced by the three lava fountain episodes and with information about the extent and position of the volcanic cloud supplied by the Toulouse VAAC during the course of the three explosive eruptions.

Mt. Etna activity
With a total elevation of > 3300 m, Mt. Etna is one of the most active volcanoes in the world. Its eruptive activity has been documented since the Holocene, leading to the creation of a rich database which is constantly updated by the INGV-OE. Among the different types of eruptive activity, in the last decades, Mt. Etna produced hundreds of "paroxysmal" episodes from the summit craters characterized by strong Strombolian and/or lava fountaining activity (Alparone et al. 2003;Corsaro et al. 2017). These episodes are frequently clustered in sequences over weeks to months. Each episode may last hours to a couple of days, starting with mild-to-intense Strombolian activity before progressing into a climatic phase. The latter is characterized by the formation of an eruption column above the volcano followed by the dispersal of a volcanic cloud into the atmosphere (Alparone et al. 2003).
A new sequence of paroxysmal episodes started in 2021 at the South-East Crater (SEC), one of the summit craters of the volcano, producing a total of 17 episodes in the period 16 February-1 April 2021. From the analysis of data coming from a multiplatform satellite system, Marchese et al. (2021) found that some of the 2021 episodes were among the most intense eruptive episodes occurring at Mt. Etna in the last 10 years in terms of thermal anomalies detected from space. Andronico et al. (2021a) classified the 2021 episodes as Lava Fountains (LF) and Large Scale Lava Fountains (LSLF). With these two terms, the authors indicate that they fall within the most common and known paroxysmal activity of Mt. Etna, consisting of sustained lava fountains producing eruption columns rising a few kilometres above the vent. Compared to the LF, the LSLF are more intense episodes leading to the formation of very high (column height > 10 km above the sea level) and tephra-rich eruption columns. At Mt. Etna, the hazard associated with this kind of summit eruption is mainly represented by the tephra fallout and the atmospheric ash dispersal. Catania, the largest city near the volcano and located about 30 km south of it, has a population of over 300,000 and it hosts infrastructures such as an international airport and railways. Thus, when explosive eruptions occur, real-time forecasts should be supplied to the local authorities and population in order to mitigate the risk due to tephra dispersal and fallout.
In the present work, we consider three episodes of the 2021 sequence and we use them as test cases for our automated workflow. The episodes are all classified as LSLF by Andronico et al. (2021a) and took place on 16 February, 28 February, and 4 March 2021. According to the weekly bulletins issued by the INGV-OE, the 16 February 2021 lava fountain was preceded by a gradually increasing Strombolian activity at the SEC. The lava fountaining started at around 17:05 (all times are UTC) (panel (a) of Fig. 1). The magma jets were estimated up to 600-700 m above the vent and produced a volcanic column reaching about 9 km above sea level (a.s.l.) that was drifted southward by the prevailing winds. Fine-grained lapilli fallout affected the city of Catania (at about 27 km from the SEC), while ash also fell in Siracusa further south (∼ 80 km), and lava flows inundated Valle del Bove and Valle del Leone in the east and north-east upper slopes of Mt. Etna, respectively. The INGV-OE surveillance cameras showed that lava fountains were over before 18:00 (INGV-OE 2021a). The 28 February 2021 episode started few minutes after 07:00, when a weak Strombolian activity resumed from the SEC and evolved very quickly into lava fountaining around 08:00. The activity reached its peak in intensity between 08:20 and 08:30, with the fountains reaching a height of 700 m and an eruption column rising up to 11 km a.s.l. (panel (b) of Fig. 1). The fountains were sustained for about 30 min and ended almost abruptly around 08:33. Tephra fallout was reported to the south-east sector, significantly covering the town of Zafferana with a few-cm thick lapilli layer (INGV-OE 2021b). The 4 March 2021 paroxysm consisted of Strombolian explosions started at 00:20, which gradually increased in intensity causing moderate, discontinuous tephra emission for a few hours. At about 08:00, this activity shifted to lava fountaining associated with a sustained eruption column of 11 km a.s.l. that dispersed tephra particles towards NE (panel (c) of Fig. 1). Abundant lapilli fell on the town of Linguaglossa, while Reggio Calabria (∼ 70 km) was covered by a thin layer of ash (INGV-OE 2021c).

Real-time workflow for short-term volcanic hazard assessment
In this section, we present a general description of the automated forecast workflow. The aim of this procedure is to forecast in real time the dispersal and deposition patterns of tephra particles emitted during the course of explosive eruptions. We use the term real time to indicate that the acquisition of the input data driving the workflow is done in real time. Once the data are available, the forecast of tephra transport and deposition for the next hours (up to 12) is released in a relatively short time (from 45 to 90 min).
The workflow is automatically initialized when notifications of tephra emission are released by the volcano observatory of reference. These notifications are in the form of bulletins called Volcano Observatory Notice for Aviation (VONA) and are described below. In this section, we also describe how the workflow operates, the numerical models used, the way in which the numerical simulations are performed (through ensembles of simulations), and the forecast products released. This workflow can be easily adapted for other volcanoes requiring real-time tephra fallout forecasts.

Volcano Observatory Notice for Aviation (VONA)
To mitigate the risk to aviation due to the presence of volcanic ash into the atmosphere, the International Civil Aviation Organization (ICAO) established that, in case of significant pre-eruptive or ongoing volcanic activity, the State volcano observatories should issue volcano-alert notifications in the VONA format (ICAO 2014). VONAs contain basic volcano information (name and location) and a summary of the current volcanic activity with focus on the presence of ash emission. Alert color codes help pilots and aviation authorities to quickly identify the status of a volcano so they can make quick decisions. The level of alert can be "green" (non eruptive state), "yellow" (volcanic unrest above the normal background level), "orange" (significant volcanic unrest with likelihood of eruption), or "red" (eruption is imminent or ongoing with significant ash emission into the atmosphere). Each VONA communicates not only the current color code, but also references to the previous one (if emitted) so that it is immediately evident if the volcanic activity increased (from green to red), decreased (from red to green), or remained on the same level but some conditions changed (for example, two consecutive red VONAs reporting a change in volcanic column height). VONAs are then disseminated to Area Control Centers (ACC), Meteorological Watch Offices (MWO), and Volcanic Ash Advisory Centers (VAACs) through email, fax, telephone, or a website.
It is important to highlight that if possible the VONA indicates the height of the eruption column, which is usually given as maximum column height above sea level. At the INGV-OE, column height can be retrieved in real time through data collected by ground-and satellite-based instruments (Scollo et al. 2019). Ground-based instruments consist of a network of video-surveillance cameras from which visible images of the eruption column are taken. Such images are calibrated through a correction procedure that accounts for the projection of the plume on the camera line of sight due to the wind, and are transmitted to the 24/7 control room of INGV-OE where the volcanologists on duty can use them to estimate the height of the column. However, depending on the meteorological and eruptive conditions, the calibrated images are not always available and, even when there is a measurement, an uncertainty of ± 500 m on the retrieved value has to be considered (Scollo et al. 2014;Corradini et al. 2018). For the eruptions considered in this work, column height was estimated using the images collected by two dedicated visible cameras. The first camera is located in Catania (about 27 km south of the SEC) and the second one near Bronte (about 14 km north-west).
In support of the camera network, the volcanologists of INGV-OE can monitor the eruptions also from the satellite data supplied by the Spinning Enhanced Visible and InfraRed Imager (SEVIRI) on board the Meteosat Second Generation (MSG) satellites . SEVIRI data are automatically downloaded and visualized in the INGV-OE control room at all times. During the explosive activity, the volcanologist on duty can estimate the height of the column in real time through a semi-automatic software processing the SEVIRI images. Overall, in case the volcanologist on duty obtains an estimation of column height, this value is included in the VONA; otherwise, the latter only reports that ash emission is ongoing, without additional information. Finally, it is worth noting that VONAs usually indicate as ash the solid particles emitted during an explosive eruption, which is here used in a general way to indicate all the particles erupted (not only those having diameter < 2 mm). However, in the text, we use the term tephra and not ash to indicate the entire size distribution of the particles feeding the volcanic column.

Automated VONA-based workflow
Numerical forecasts of tephra atmospheric dispersal and ground deposition can be improved if TTDMs are initialized with information coming from real-time observations of the current state of the volcanic activity. Nowadays, eruption timeline (eruption start and end time) and column height can often be supplied in real time during the course of an eruptive episode. Other important properties of the eruptive mixture necessary to initialize the TTDMs (e.g., size, density, and shape of the solid particles erupted, water content, exit velocity of the mixture) are usually retrieved after extensive analyses of the erupted products (Bonadonna and Houghton 2005). Recent works showed that data from ground-and satellite-based sensors (such as thermal cameras and weather radars) can provide in real time the size distribution, velocity, and concentration of the pyroclasts leaving the vent (Pioli and Harris 2019;Freret-Lorgeril et al. 2021;Syarifuddin et al. 2021). However, such methodologies are still under development and not yet operational.
To produce real-time forecasts, we based our automated workflow on the eruption timeline and column heights (if given) supplied by VONAs, while the other eruptive source parameters necessary to run the TTDM are from predefined eruptive scenarios (more details are provided later in the text). We set and tested the workflow considering Mt. Etna as a case study and, for this reason, we refer to the VONAs released by the INGV-OE.
The workflow automatically starts when a red VONA having orange as previous color is issued by the INGV-OE. At this point, a series of forecast cycles are performed in order to produce probabilistic maps of the volcanic cloud in the atmosphere and the tephra deposit on the ground. Such forecasts cover a time slot of 12 h since the emission of each VONA (i.e., for each VONA, we generate forecasts up to 12 h in the future). Thus, the total time covered by the forecasts depends on the release of the VONAs. The workflow is made of a sequence of scripts written in the Python programming language and run on the computer facilities of INGV-Sezione di Pisa. Such scripts allow the detection of new VONAs, the automated execution of the numerical models and the post-processing of the results ending with the production of the probabilistic hazard maps. The state of the simulations are saved every 10 min and can be used as a restart state for a new simulation with different conditions, for example, when a change in plume height is indicated by a new red VONA.
The script checks continuously for the release of VONAs, and in particular for a red VONA issued after a previous orange one, indicating the beginning of an explosive activity producing significant emissions of tephra. This is done by background executing every 2 min a script that controls an email address receiving the VONAs from the INGV-OE. In case a red VONA is found in the inbox, the workflow is activated. The meteorological data necessary to run the simulations is computed in advance twice a day at hours 00 and 12 by downloading and post-processing high-resolution weather forecast data (further details are provided later in the text).
In the paper, we label a sequence of N + 1 VONA emission times during an eruptive episode with a superscript denoting the color and a subscript going from 0 (for example, t Red 0 is the time of the first red VONA) to N (t Or N is the time of the last orange VONA). In the simplest case, two VONAs are issued during an eruption. These are the initial red VONA and a final one, indicating the end of tephra emission in the atmosphere, for example, an orange VONA. In a more complex case, after the issue of the first red VONA, a sequence of red VONAs can be released if the eruptive conditions change (e.g., a variation in column height) or the visibility of the activity varies (e.g., going from day to night). Figure 2 panel (a) shows an example of the forecast cycles performed when only two VONAs are issued during an eruptive episode. Here, t Red 0 represents the time of the red VONA indicating the beginning of tephra emission, while t Or 1 is the time of the first VONA with an update, in this case a change in the level of alert from red to orange, corresponding to the end of tephra released in the atmosphere. In this example, the workflow is activated at time t Red 0 and a first forecast cycle (Initial cycle) is performed using t Red 0 as eruption start time and assuming continuous tephra emission for 12 h. The height of the eruption column is set according to the value given by the VONA (if given) or according to a predefined eruptive scenario. During this Initial cycle, a total of 6 forecasts are produced at hours (with respect to t Red 0 ): +1h (Forecast F 0 1 ), +2h (F 0 2 ), +3h (F 0 3 ), +6h (F 0 4 ), +9h (F 0 5 ), and +12h (F 0 6 ). All the forecast cycles performed when a red VONA is issued produce this sequence of 6 forecasts. When the orange VONA is issued, the workflow searches for the to F 1 4 highlights the final forecasts performed considering the two VONAs. On the contrary, the disregarded forecasts are marked with a dashed cycle (F 0 5 and F 0 6 ). Panel (b) replicates the same scheme of panel (a), but for a more complex case where two red VONAs reporting a change in column height from H 0 to H 1 are issued during an eruptive episode. In this case, three forecast cycles are executed, named Initial cycle (the first cycle with eruption column height H 0 ), Update 1 (which considers the change in column height from H 0 to H 1 ), and Update 2 (where tephra emission ceases at time t Or 2 ) last state (saved with a frequency of 10 min) of Initial cycle produced before t Or 1 (here the state at t Red 0 +7h) and uses this state as initial condition for a new forecast cycle (Update 1). In panel (a), forecasts F 0 1 , F 0 2 , F 0 3 , and F 0 4 of Initial cycle are still valid after the issue of the orange VONA, while forecasts F 0 5 and F 0 6 are disregarded. Four new forecasts every 3 h (from F 1 1 to F 1 4 ), covering a time span of 12 h after the last valid forecast of the Initial cycle, are then produced. Because of the orange level of the VONA, for the simulations of this new cycle of forecasts, no tephra is emitted after time t Or 1 . In case the eruptive conditions change during the episode (e.g., with an increase or decrease in column height), a sequence of red VONAs can be issued and, in this case, the workflow updates the forecasts in real time with the latest information available from the VONA. This example is shown in panel (b) of Fig. 2. Here, the first forecast cycle (Initial cycle) starts when the first red VONA reporting a column height H 0 is released at time t Red 0 . Forecasts F 0 1 , F 0 2 , F 0 3 , F 0 4 , F 0 5 , and F 0 6 are obtained assuming a maximum column height above the sea level equal to H 0 . At time t Red 1 a new red VONA is issued reporting a new column height equal to H 1 . At this point, the workflow searches for the latest state of Initial cycle saved before t Red 1 (here the state saved at t Red 0 + 4h 30min) and starts a new cycle (Update 1) that produces 6 new forecasts (from F 1 1 to F 1 6 ) where the height of the eruption column is H 0 from t Red 0 to t Red 1 and H 1 starting from t Red 1 . This means that the workflow keeps track of all the particles released from t Red 0 to t Red 1 assuming a column height H 0 and, starting from t Red 1 , releases new particles assuming a column height H 1 . In this way, the forecast F 1 1 is done according to the eruption dynamics as described by the two red VONAs issued. The forecasts F 0 1 , F 0 2 , and F 0 3 of Initial cycle are still valid, while F 0 4 , F 0 5 , and F 0 6 are disregarded since they do not account for the change in column height. When the final orange VONA is issued at time t Or 2 , the workflow operates as described for panel (a) and a new cycle Update 2 starts producing new forecasts (from F 2 1 to F 2 4 ) where the tephra emission ends at t Or 2 . In the previous examples, after the completion of each forecast cycle (I nitial cycle to Update N), the workflow stays in waiting mode until a new VONA is found by the script checking the inbox. At this point, the workflow is reactivated and proceeds until a new waiting mode (if the latest VONA found is a new red one) or the end (if the latest VONA found is the orange one). In the eventuality that a new VONA is found while a forecast cycle is still running, the workflow terminates the ongoing simulations and post-processing phase, produces the hazard maps, and skips to the next forecast cycle without concluding all the 6 expected forecasts. For a specific eruptive episode, the workflow terminates when all the 4 forecasts referring to the final VONA indicating the end of tephra emission in the atmosphere are performed.

Numerical models: PLUME-MoM-TSM and HYSPLIT
The numerical model used to simulate the release and transport of volcanic particles in the atmosphere is the HYbrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model (Stein et al. 2015). HYSPLIT, developed by the NOAA Air Resources Laboratory's (ARL), is a widely used software for the computation of both simple air parcel trajectories and more complex scenarios involving the transport, dispersal, and deposition of solid particles and gases. To date, HYSPLIT is the TTDM used by the VAACs of Darwin (Australia), Wellington (New Zealand), Washington, and Anchorage (USA).
To better reproduce the physics of volcanic eruptions, we coupled HYSPLIT with the latest release of the eruption column model PLUME-MoM named PLUME-MoM-TSM (de' Michieli Vitturi and Pardini 2021), in a way that the output of PLUME-MoM-TSM in terms of mass fluxes of tephra particles lost from the column during the rising and mass injected at the neutral buoyancy level are used to initialize HYSPLIT. The coupling between the two codes is automated and based on Python routines which automatically produce the input files of HYSPLIT from the output of PLUME-MoM-TSM. A detailed description and applications of this procedure can be found in Pardini et al. (2020), Tadini et al. (2020), and Andronico et al. (2021b).
The use of PLUME-MoM-TSM allows us to also simulate the initial density-driven spreading of the umbrella cloud at the neutral buoyancy level. Usually, 1D column models do not simulate the gravity current effects for the spreading of the umbrella region or they use simplified approaches (Costa et al. 2016). PLUME-MoM-TSM simulates the umbrella spreading by solving a set of shallow water equations that allow us to describe the regime of gravity current that characterizes the horizontal movement of the volcanic mixture at the neutral buoyancy level. A proper description of this feature expands the use of the model as, for example, umbrella growing can help to constrain the eruptive source parameters especially for large explosive eruptions (Pouget et al. 2013;Constantinescu et al. 2021). Moreover, Webster et al. (2020) showed that to accurately forecast the transport of tephra from large eruptions, the lateral spread of umbrella clouds needs to be represented within TTDMs.
Below the neutral buoyancy level, the tephra particles decoupled from the column enter the atmosphere at a number of heights located from the vent to the neutral buoyancy level itself. At each height, we inject the particles into the atmosphere from a number of point sources positioned at the edge of the horizontal cross-sections of the column. The finest particles reaching the neutral buoyancy height enter the atmosphere from a circular area source representing the position and dimension of the umbrella cloud. The release of tephra particles from this area allows a better reproduction of the volcanic cloud and of the upwind deposit, which is difficult to reproduce when tephra enters the atmosphere downwind only. The detailed description of the definition and calculation of the umbrella cloud spreading is provided in the Appendix.

Ensemble-based forecasting of volcanic particles
In the automated workflow, we simulate the dispersal and deposition of volcanic particles by accounting for uncertainties in both meteorological and volcanological input data. To do this, we perform ensemble-based simulations in which each member of the ensemble is defined by a perturbed version of the original meteorological data and a specific set of Eruptive Source Parameters. Thus, we characterize the tephra deposit and the atmospheric volcanic cloud in a statistical way, by computing the cumulative distribution functions of the quantities of interest and the probabilities of exceeding prescribed thresholds of ground load or atmospheric concentration.
As presented in Pardini et al. (2020), the uncertainty of the meteorological data is modelled by generating N perturbed versions of a reference meteorological data. Each perturbed meteorological data derives from the original one by modifying intensity and direction of all the horizontal wind vectors forming the meteorological grid. We perturb only the horizontal wind field, while the rest of the quantities (i.e., vertical velocity, temperature, humidity profiles) remain unchanged. We define two probability distributions representing the uncertainties in intensity and direction of the reference wind field and we apply Latin Hypercube Sampling (LHS) to generate N near-random samples from the two distributions. Each sample is a couple of numbers representing the perturbations in wind direction and intensity used to modify the original meteorological data. This operation ends with the production of N meteorological files, each used as a member of the ensemble. Compared to traditional Monte Carlo sampling, LHS requires a lower number of samples to well represent the variability of the parameter space, meaning a lower number of members forming the ensemble and thus lower computational times, which is of fundamental importance for a real-time application. A description of LHS is given in the Supplementary Material.
We also consider uncertainties in the volcanological parameters defining the eruptive conditions at the vent. In particular, the height of the eruption column and the Total Grain Size Distribution (TGSD) of the tephra particles erupted are key parameters for the initialization of any TTDM (Pardini et al. 2016;Dioguardi et al. 2020). Indeed, a value of plume height far from the real one produces unreliable forecasts since different column heights (and thus mass flow rates) result in different tephra dispersal and depositional patterns. On the other hand, TGSD does not significantly influence the mass flow rate, but it determines the terminal settling velocity of the tephra particles erupted, which is a key factor in controlling the residence time of tephra in the atmosphere and thus the distance reached. Differently from column height, TGSD is not currently included in the VONA. In the workflow, we define the TGSD as a log-normal distribution through the mean value μ and the standard deviation σ . With the aim of quantifying the effects that such uncertainties have on volcanic hazard, we define column height and TGSD (μ and σ ) through probability distributions. We then sample such space with LHS and we produce a set of N triplets column height, μ and σ , each defining a different eruptive condition at the vent.
To generate the members of the ensemble, we associate to each of the N eruptive triplets a meteorological data randomly chosen from the N previously defined, in a way that all the N meteorological data are used as ensemble members. In this way, we extend the LHS approach adopted for the eruptive parameters to an additional dimension, given by the meteorological data, and we end with N members forming the ensemble. Each member presents a specific set of eruptive and meteorological conditions.

Forecast products
The ensemble of simulations allows us to produce statistical outcomes of the quantities of interest. We look at both the volcanic cloud and the tephra deposit and we produce volcanic hazard maps together with information on tephra load at specific locations.
The hazard due to the presence of tephra into the atmosphere is modelled through maps reporting the probability to exceed a certain level of tephra contamination between two altitudes. According to the ICAO EUR/NAT Volcanic Ash Contingency Plan (ICAO and VACP 2016), three contamination levels are defined: low contamination (volcanic particle mass concentration ≥ 0.2 mg/m 3 , ≤ 2 mg/m 3 ), medium contamination (> 2 mg/m 3 , < 4 mg/m 3 ), and high contamination (≥ 4 mg/m 3 ). These contamination levels aim at guiding the airlines on the decisions to take in case of ash-contaminated fly zones. Indeed, except for the high contamination level, flight operations can be authorized under certain conditions of accurate pre-and post-flight aircraft inspections in the low and medium contamination zones. Thus, we produce three hazard maps reporting the probability of exceeding the thresholds of 0.2 mg/m 3 , 2 mg/m 3 , and 4 mg/m 3 of particle mass concentration at different atmospheric levels. In particular, we defined three atmospheric intervals: from 0 km to 5 km a.s.l., from 5 km to 10 km a.s.l., and from 10 km to 20 km a.s.l. and, for each atmospheric interval, we compute tephra concentration defined as the total mass of tephra inside each 3D grid cell divided by the volume of the cell. For each atmospheric interval, the volume of a cell is defined through its horizontal extents in the x and y directions and the height of the atmospheric interval.
For the deposit, we present two hazard maps showing the probability of exceeding a tephra load of 0.5 kg/m 2 and 5 kg/m 2 . Nowadays, tephra load is the preferred measure of the severity of tephra fallout and it is a common output of the TTDMs. Moreover, to better quantify the hazard due to tephra fallout, we selected 10 key sites which are the most vulnerable towns located within 15 km from the vent (Nicolosi, Zafferana, Linguaglossa, Randazzo, Bronte, and Biancavilla), and the main cities of Catania (about 27 km south of the SEC), Palermo (about 150 km), Siracusa (about 80 km), and Reggio Calabria (about 70 km) (see Fig. 3). We highlight that airports are located in Catania Fig. 3 The map indicates the key sites where the forecast workflow computes the ground load at the different forecast times. We considered the main cities of Catania, Siracusa, Reggio Calabria, and Palermo and the main towns located at the slopes of Mt. Etna (Nicolosi, Zafferana, Linguaglossa, Randazzo, Bronte, and Biancavilla) (Fontanarossa International Airport), Palermo, and Reggio Calabria. For each location, we compute the 5th, 50th, and 95th percentiles of the ground load in kg/m 2 . Finally, we compute the same percentiles of mass eruption rate, total amount of tephra into the air and deposited on the ground.
For each eruptive episode, hazard maps and information on erupted products are uploaded in real time on a web page contained in the git-hub repository: EtnaFountains2021. Each eruptive episode is labelled with the time of issue of the first red VONA referring to it and the web page is updated in real time as soon as new forecast products are produced by the workflow.

Application to the 2021 Mt. Etna lava fountain episodes
This section describes how we set and calibrated the automated workflow for the three case studies at Mt. Etna. In particular, we describe the setting of the numerical models (PLUME-MoM-TSM and HYSPLIT) and the results of a sensitivity analysis performed with the aim of understanding which of the uncertain inputs (i.e., wind parameters, column height, and TGSD) has the greatest control on tephra load at specific locations.

Setting PLUME-MoM-TSM and HYSPLIT
The meteorological data used to run the numerical simulations are supplied by the HydroMeteorological Service of the Emilia-Romagna Regional Agency for Environmental Protection. The files are in GRIB format, contain 72 h of weather forecast with a time step of 3 h and the horizontal grid resolution is 0.045 • (see the Supplementary Material for a detailed description). The probability distributions defining the perturbations on the wind field are normal distributions with standard deviations of 50% for wind intensity and 20 • for wind direction .
The inputs of the column model (PLUME-MoM-TSM) and the TTDM (HYSPLIT) are defined both as deterministic values and probability distribution functions. The distribution of the pyroclastic particle size is described using the Krumbein scale defined as: where D mm is the particle diameter expressed in mm and D 0 is a reference diameter equal to 1 mm. The TGSD is modelled as a normal distribution in the Krumbein scale whose parameters μ (mean value) and σ (standard deviation) are given as uniform distributions ranging from [−1.  (Aubry et al. 2021). In particular, the archive contains the TGSD of some Mt. Etna eruptive episodes spanning from 1998 to 2016. We computed the mean value and the standard deviation of each TGSD supplied by the archive and we defined the ranges of variation of μ and σ by taking the minimum and maximum values computed for the two variables.
The other input variable defined as probability distribution is the height of the volcanic column that we take as the maximum height above the sea level. In case the VONA includes an estimate of the column height, we use this value to initialize the column model by assuming a Gaussian distribution having as mean value the height supplied by the VONA and 500 m as standard deviation. On the contrary, if column height is not supplied by the VONA, we define it as a uniform distribution covering the range [5; 15] km a.s.l.. Such variability has been derived by looking at the column heights of the Mt. Etna explosive eruptions reported in the IVESPA archive with the addition of some data referring to the 2015 explosive events reporting a column height of about 15 km a.s.l. .
We remark that we initialize PLUME-MoM-TSM with the height of the eruption column and not the mass flow rate, which is computed by the model. To do this, we give in input the velocity of the eruptive mixture at the vent (100 m/s) and an interval for the initial radius ([10; 200] m). By sampling the radius interval, the model searches for the value giving the desired plume height. The corresponding mass flow rate is then computed from the initial velocity, radius, and mixture density. Since PLUME-MoM-TSM is a wind-affected model, the resulting mass flow rate depends on the atmospheric profile above the vent (i.e., atmospheric temperature, density, humidity, and wind velocity). This implies that, even if column height does not change during an eruption, variations in the meteorological conditions affect the mass flow rate computed by the model. As shown by Dioguardi et al. (2020), the calculation of mass flow rate from column height is strongly affected by the model used (no wind or wind-affected) and it influences the forecasts produced by the TTDMs (see the Supplementary Material for more details on the setting of the numerical models).
In the workflow, each ensemble of simulations is formed by a total of 21 members. Calibration tests were performed in advance to investigate the sensitivity of the results to the number of members forming the ensemble. In the calibration phase, we searched for the minimum number of members to produce stable results (i.e., results do not change significantly by increasing the number of members forming the ensemble) and to maintain the workflow suitable for real-time applications (i.e., low computational times and the possibility to run the simulations without queuing on the computer facilities of INGV-Sezione di Pisa). The final value of 21 members represents the best choice in terms of numerical results and computational requirements.

Sensitivity analysis
As a preliminary investigation, we performed a sensitivity analysis to study the role that each uncertain model input plays on the variability of model outputs. The investigated model inputs are both Wind Parameters (WP) and Eruptive Source Parameters (ESP). To do this, we employed the toolkit Dakota, which is an open source software used to carry out optimization and uncertainty quantification analyses (Adams et al. 2014). With Dakota, we performed ensembles of simulations considering the uncertain model inputs reported in Table 1, which are the same uncertain inputs used in the workflow for the production of the hazard maps.
For this analysis, we set PLUME-MoM-TSM and HYSPLIT to reproduce the 4 March 2021 lava fountain episode, and we assumed that the eruption started at 8:00 and ended at 10:00 on the same day (i.e., 4 March 2021). This time interval captures the most intense phase of the eruption with an observed column height of 11 km a.s.l.. Tephra dispersal and deposition is modelled until 16:00 on 4 March 2021. The investigated model output is the ground load at two locations chosen in the direction of the main plume dispersal axis. The first point is located at a distance of 7 km from the vent, while the second one is positioned at more than 40 km from it. With Dakota, we performed an ensemble of simulations consisting of a total of 4900 runs. Indeed, the sensitivity analysis requires hundreds to thousands of simulations according to the formula N(dim + 2), where dim is the number of variables (5 in our case) and N is the number of samples for each variables (700) (Adams et al. 2014).  Figure 4 shows the results of the simulations. The left plots refer to GL1, the ground load at the first point, while the right plots refer to GL2, the ground load at the second point. In each column, the first five plots from top to bottom show the variation of GL1 and GL2 as a function of each of the five inputs δ d , δ i , CH , μ, and σ (see Table 1). First of all, we observe a larger variability in the ground load at the farthest point, in the range 10 −5 − 10 3 kg/m 2 , while for the point closer to the vent most of the values exceed 10 −1 kg/m 2 . It appears that, in this analysis, column height is the variable exerting the greatest control on both GL1 and GL2, with lower heights resulting in smaller loading at the two selected locations. This is because column height dictates the mass flux computed by PLUME-MoM-TSM. Furthermore, for the farthest point, a larger spreading in the ground loading values, and thus a larger sensitivity to the other input parameters, is observed for smaller column heights. A weaker correlation emerges between the outputs and μ, while no correlation is present between the ground load and the variables δ d , δ i , and σ . The last two plots of panels (a) and (b) of Fig. 4 show the Probability Histogram and the Cumulative Distribution Function of GL1 and GL2. Probability is 48.5% that GL2 is in the range [0; 2] kg/m 2 , while this value is equal to 21.8% for GL1. The 5th, 50th, and 95th percentiles of GL1 are 0.3kg/m 2 , 8.60kg/m 2 , and 117.70 kg/m 2 , while the same quantities are 0kg/m 2 , 2.30kg/m 2 , and 79.20kg/m 2 for GL2.
To better quantify the fraction of output variability to ascribe to the variability of specific inputs or to the interactions between them, we also performed a variancebased decomposition analysis, in order to compute the Main Effect Index and Total Effect Index (also known as Sobol indices). Considering the output function Y and the ith input x i , the Main Effect Index of the first order (S i ) indicates the fraction of the output variability due to the variability of x i only: where Var x i [E(Y |x i )] is the variance of the conditional expectation and Var(Y ) is the variance of the output function. Higher order interaction indices can be computed and they represent the fraction of output variability due to interaction between two (second order) or more inputs. The Total Effect Index (T i ) indicates the fraction of Y variability due to x i and to its interaction with the other input variables: In both cases, higher indices indicate that the variability in the input parameter has a greater control on the variability of the output.
To compute the Sobol indices through Eqs. 2 and 3, it is necessary to calculate the variances of the output functions. Here, we performed a Generalized Polynomial Chaos Expansion (gPCE) in order to create a surrogate model of the ground load, from which a very large number of output evaluations can be performed with low computational cost (Crestaux et al. 2009;Pardini et al. 2016). Figure 5 shows the main and total Sobol indices computed for GL1 and GL2. From the main indices (panel (a)), it appears that CH is the parameter with the greatest control on the variability of both GL1 and GL2 (in accordance with the findings shown of Fig. 4). After CH , the wind parameters δ d and δ i exert the main role in controlling the values of GL1 and GL2 (especially for the proximal point). The eruptive source parameters μ and σ seem to have a lower importance. From panel (a), it also appears that interactions between eruptive and wind parameters are important in determining the variability of GL1 and GL2. Indeed, the sum of the higher-order order indices contributes to the 60% and 55% of the total variability of GL1 and GL2 respectively. The Total Sobol Indices highlight again the role of CH and its interactions with the other inputs in controlling GL1 and GL2, panel (b) of Fig. 5, and the fact that none of the uncertainties here considered should be disregarded. For this reason, in the following analysis and in the workflow for the assessment of the tephra hazard, the uncertainties in all these parameters are always considered.

Numerical results and hazard maps
As it can be seen in Table 2, a sequence of red VONAs were issued for each investigated eruption. Some of the VONAs report a change in column height, while others confirm the previous height based on new observations (INGV-OE 2021a, b, c).
As a preliminary result, we show how the evaluation of the volcanic hazard is heavily dependent on the information on column height possibly supplied by the VONA, by  (panel (a)) and Total Effect Indices (panel (b)) for the ground loads GL1 and GL2 evaluated at 7 km and 40 km from the vent, respectively. WP indicates the uncertain wind parameters δ i and δ d , while ESP includes the uncertain eruptive source parameters CH , μ, and σ  Fig. 6 panel (a) shows the forecast at 16:30 obtained assuming as column height the value supplied by the first red VONA (6 km a.s.l.). This means that a constant column height of 6 km a.s.l. was assumed from the beginning of the eruption (i.e., t Red 0 =08:11), up to the considered forecast time. With this assumption, the simulated deposit exceeds 5 kg/m 2 in the proximal area only. For this ensemble of simulations, the 50th percentile of the ground load computed in Catania is 2.29×10 −3 kg/m 2 . A different hazard scenario is presented in panel (b), where the hazard map is representative of the same forecast time, but including information also from the second VONA issued at t Red 1 , i.e., from two consecutive red VONAs. In this case, we assumed the height of the eruption column changing from 6 km to 9 km a.s.l. at t Red 1 =8:35 (t Red 0 + 24 min) and then being constant until 16:30. The ground loads for this eruptive scenario are higher than those from the previous one, with the 50th of expected ground load in Catania equal to 5.56×10 −1 kg/m 2 . Panel (c) shows the hazard maps at 16:30 obtained according to the three red VONAs issued at t Red 0 , t Red 1 , and t Red 2 . The last VONA reports a decrease in column height from 9 to 4.5 km at t Red 2 =10:47 (t Red 0 + 2 h 36 min) and this change in eruptive conditions leads to a ground load in Catania lower than the one obtained by considering the first two VONAs only and with the 50th percentile equal to 6.74 × 10 −2 kg/m 2 . Finally, the hazard scenario presented in panel (d) was obtained assuming that no information on the height of the eruption column was supplied by the VONAs. In this case, the workflow assigns to the column height an uniform probability distribution ranging from 5 to 15 km a.s.l.. For this scenario, the 50th percentile of the expected ground load in Catania is 2.78 kg/m 2 .
The different choices for the definition of column height during the eruption (same as Fig. 6) have an effect also on the concentration of tephra in the atmosphere. Figure 7 presents the different hazard maps for the atmospheric ash concentration evaluated at 16:30 on 28 February 2021 corresponding to the four scenarios presented above. It emerges that the scenario where the ash contamination is more extended (both in the horizontal and in the vertical Fig. 6 Hazard maps for the tephra deposit obtained for different assumptions on the height of the eruption column. We present the probability to exceed a ground load of 5 kg/m 2 . The 28 February 2021 eruption is considered as test case. Panel (a) shows the hazard map obtained according to the column height given by the red VONA issued at t Red 0 (6 km). Panel (b) includes information from VONAs issued at times t Red 0 and t Red 1 (9 km), while panel (c) considers information coming from the three red VONA issued at times t Red 0 , t Red 1 , and t Red 2 (4.5 km). Finally, panel (d) presents the hazard map on the ground obtained assuming that no information on the height of the eruption column is supplied by the VONA. In this case, we set a column height ranging from 5 km to 15 km a.s.l. (values chosen from historical records) Fig. 6, but for the hazard maps produced for the presence of volcanic particles into the atmosphere. Panel (a) shows the forecast at 16:30 done considering that the eruption column is 6 km a.s.l.; in panel (b), the same forecast is done assuming that the eruption column increases in height from 6 to 9 km at t Red 1 , while panel (c) considers the complete evolution of the column which decreases to 4.5 km at t Red 2 . Finally, panel (d) presents the forecast produced assuming as column height a uniform distribution from 5 to 15 km a.s.l. directions) is the one obtained when no information on column height is given by the VONAs, for which the column height, and thus the mass flow rate, are overestimated.

The 16 February 2021 eruption
In terms of VONAs, the eruption started at t Red 0 =16:31 with the issue of the first red VONA reporting a column height of 6 km a.s.l.. As a consequence, at that time, the workflow was activated and the first forecast cycle (I nitial cycle) was performed with the production of forecasts from F 0 1 to F 0 5 (see Fig. 8). At t Red 1 =17:05, a new red VONA reporting a change in column height (from 6 km to 10 km a.s.l.) was issued. Because of that, a new forecast cycle (Update 1) started at t Red 1 assuming as column height 10 km a.s.l. and taking as initial condition the last state of I nitial cycle saved before t Red 1 . In this way, the new forecast cycle was performed according to the eruption dynamics as described =17:38, reporting the continuous emission of tephra, but without information on the column height. Thus, numerical simulations were restarted from the last state of the cycle Update 1 saved before t Red 2 and assuming for the column height an uniform distribution between 5 and 15 km a.s.l. (Update 2). With this last information, the workflow produced the new forecasts from F 2 1 to F 2 6 . When the orange VONA stating the end of the explosive episode was issued at time t Or 3 =05:31 on 17 February 2021, the last forecast cycle (Update 3) was performed assuming no tephra emission, but keeping track of the particles emitted from t Red 0 to t Or 3 . In Fig. 8, the forecasts of a cycle at times subsequent to the release of a new VONA, or those still not produced because of its issue, are plotted with dashed lines, while those kept, which takes into account all the information available, at plotted with a solid line.
The complete set of hazard maps obtained for this eruption is available at EtnaFountains2021/16Feb2021. Here, in Fig. 9, we show the probability maps of tephra ground load to exceed 0.5 kg/m 2 (left plots) and 5 kg/m 2 (right plots) for two different times. Furthermore, three contour lines indicate the 5%, 50%, and 95% probabilities to exceed the considered threshold on ground load. In a similar way, Fig. 10 shows hazard maps reporting the probability to exceed the ash concentration threshold of 2 mg/m 3 at three atmospheric intervals. Finally, for this eruption, Table 3 column 16 Feb 2021 collects the cumulative ground load at the end of the deposition process (i.e., 17:30 on 17 February 2021) at the ten key locations. Both the figures and the table show that most of the tephra is expected to be transported to the South, with the sites of Nicolosi, Catania, and Siracusa (all located on the main dispersal axis) having the 50th percentile of cumulative ground load equal to 31.6 kg/m 2 , 24.7 kg/m 2 , and 25.3 kg/m 2 . For this eruption, the ground load at the other locations is negligible. The ash cloud is transported and dispersed by the wind towards South and ash contamination is present also in the atmospheric interval [10; 20] km.

The 28 February 2021 eruption
The results for the 28 February 2021 eruption are presented in a similar way to what was done for the previous case. Figure 11 summarizes all the forecast cycles performed after the issue of the first red VONA at t Red 0 =08:11 on 28 February 2021. In this case, a total of 4 forecast cycles were performed according to the 4 VONAs issued during the eruption (Table 2). In particular, a second red VONA reporting an increase in column height from 6 km to 9 km a.s.l. was issued at t Red 1 =08:35 (24 min after the first VONA). The intensity of the eruption maintained on the same level until t Red 2 =10:47, when the last red VONA was issued showing a column height of 4.5 km a.s.l.. Finally, the orange VONA was released at t Or 3 =14:26. Figures 12 and 13 present some of the hazard maps obtained for tephra fallout and atmospheric dispersal, while the complete set of forecasts for this episode can be found at the page: EtnaFountains2021/28Feb2021. Also for this case, column 28 Feb 2021 of Table 3 collects the cumulative ground load computed at the end of the deposition process at 10 key locations. It appears that the fallout from plume of the 28 February lava fountain Fig. 9 Hazard maps for the 16 February eruptive episode showing the probability to exceed a tephra ground load of 0.5 kg/m 2 (left plot) and 5 kg/m 2 (right plot). The three contours indicate the 5% (black line), 50% (red line), and 95% (blue line) probabilities to exceed the considered threshold on ground load. Panel (a) shows the deposit evaluated at t Red 2 + 1h, while panel (b) presents the deposit after the issue of the orange VONA (t Or 3 + 3h). Each plot also shows the position and the label of the sites where tephra load is evaluated. Load values are collected in Table 3 affected mainly the town of Zafferana (located on the main dispersal axis) with the 50th percentile of ground load equal to 5.12 kg/m 2 . The towns of Nicolosi and Catania (both located on the right of the dispersal axis) show a ground load of 0.471 kg/m 2 and 0.125 kg/m 2 , respectively. The other locations present negligible loads all equal or less than 0.001 kg/m 2 . Following the atmospheric wind, ash is initially dispersed eastward, before moving southward due to a rotation in wind direction from East to South.

The 4 March 2021 eruption
The complete set of forecast cycles performed for the eruption occurred on 4 March 2021 is presented in Fig. 14. Also in this case, the workflow operated according to the VONAs released during the episode and listed in Table 2. As shown by Fig. 14, a total of 5 cycles were performed according to the 5 VONAs issued during the episode (5 red VONAs and the final orange one). The first VONA was issued during the night at t Red 0 =02:28 reporting tephra emission from a column of 5 km a.s.l.. Column height maintained stable until t Red 1 =07:51, when the second red VONA was issued reporting an increase in column height from 5 km to 6 km a.s.l.. At this point, the episode rapidly reached its peak in intensity with the release of the third and the fourth VONAs at t Red 2 =07:55 and t Red 3 =08:17, respectively. Both report an increasing column height that reached 11 km a.s.l. at t Red 3 . This paroxysmal phase lasted Fig. 10 Hazard maps reporting the probability to exceed the ash concentration threshold of 2 mg/m 3 at three atmospheric intervals (0 − 5 km, 5 − 10 km, 10 − 20 km). The three contours indicate the 5% (black line), 50% (red line), and 95% (blue line) probabilities to exceed the considered threshold on ash concentration. The maps refer to the 16 February 2021 lava fountain episode. As done for the deposit, we present forecast times t Red 2 + 1h (panel (a)) and t Or 3 + 3h (panel (b)) about 2 h, and at t Or 4 =10:16, the orange VONA was issued. It is worth noting that the VONA issued at 06:16 was not considered by the workflow since it reported the same column height of the previous red one. The gith-hub page EtnaFountains2021/4Mar2021 contains all the forecasts generated for this lava fountain. Figures 15 and 16 present some of the hazard maps produced for this episode, while Table 3, column 4 Mar 2021, reports the final cumulative ground load at the 10 considered locations. Tephra particles produced by this episode dispersed and deposited towards North East, with the town of Linguaglossa having the 50th percentile of the ground load equal to 9.33 kg/m 2 . Compared to the episodes of the 16 and 28 February, the 4 March lava fountain had the highest column height For each investigated lava fountain episode, we show the 5th, 50th, and 95th percentiles of the cumulative ground load at the end of the deposition process (kg/m 2 ). Labels indicate: 1 Catania, 2 Siracusa, 3 Reggio Calabria, 4 Palermo, 5 Nicolosi, 6 Zafferana, 7 Linguaglossa, 8 Randazzo, 9 Bronte, 10 Biancavilla From our simulations, it emerges that none of the three lava fountain episodes produced tephra fallout affecting Palermo (the biggest town in Sicily). On the contrary, the town of Catania (where the Fontanarossa International Airport is located) can be impacted by tephra fallout depending on atmospheric and eruptive conditions. Indeed, the 16 February episode appears to be the most severe for the town of Catania, while no tephra fallout is foreseen from the 4 March eruption.

Testing the numerical results
In this section, we show how we tested the numerical forecasts by comparing some of the results of the workflow with independent observations of the deposits at selected locations and with the extent of the volcanic clouds as reported by the Toulouse VAAC.

Comparison with the observed deposits
Samples were collected on a measured area selecting sites where the tephra deposit was clearly primary and not weathered by wind, rain, or anthropic agents. The sampling points were chosen in order to have a spatial distribution as regular as possible on the ground deposit and at the same time to highlight the difference in mass, thickness, or covered surface from site to site. Samples were then dried in the INGV-OE Laboratory of Sedimentology at 60 • C and the final weight extrapolated to unit area (1 m 2 ). Figure 17 panel (a) shows the sampling points considered for the three eruptions. Overall, the observed ground load is in the range 5 − 10 kg/m 2 for the sampling points located near the vent and on the main dispersal axis. For the most distal points (up to 40 km from the vent), values less than 1 kg/m 2 are retrieved (the minimum load is 6 · 10 −3 kg/m 2 observed for the 16 February 2021 episode). To compare the observed ground load with the simulated one, we computed at the field sampling points of each episode the 5th, 50th, and 95th percentiles of the ground load and plotted the results in Fig. 17 panel (b). In the panel, the colored diamonds refer to the 50th percentile, while the error bars to the 5th (lower limit) and 95th (upper limit) percentiles. The continuous black line indicates a perfect match between the measured and the simulated loads, while the dashed and point lines correspond to factors 5 (1/5) and 10 (1/10), respectively. For the 16 February 2021 episode (blue diamonds), it can be seen that the simulated ground load computed at the end of the deposition process (at 17:30 on the 17 February 2021) is orders of magnitude higher than the observed one (see later in the Discussions). For the 28 February episode, the simulated load well reflects the measured one (red diamonds of Fig. 17 panel (b)). Indeed, the large majority of the simulated loads lies in the interval 1/5-5. An analogous result emerges for the 4 March episode (black diamonds of  Fig. 17 panel (b)). Also in this case, all the simulated loads (except one) intersect or completely lie in the 1/5-5 interval.

Comparison with VAAC bulletins
When VONAs are issued by the INGV-OE, the Toulouse VAAC is in charge of providing Volcanic Ash Advisories (VAAs) and Graphics (VAGs) to inform about the ongoing volcanic activity. VAAs are text messages produced within 20 min from the notification of an ongoing eruption and provide all the information necessary to define the extent and movement of the volcanic cloud . The data sources of the VAAs supplied for Mt. Etna are images from the ground-based monitoring network of the INGV-OE, satellite imagery, and VONAs. VAGs translate in graphic format the information supplied by VAAs, and contain contour lines showing the horizontal extent of the volcanic cloud together with its total vertical extent in flight levels. Flight levels are specified since such bulletins are primarily used for aviation risk mitigation. Both VAAs and VAGs also supply qualitative forecasts of the movement of the cloud.
Similarly to what was done by Zidikheri et al. (2018) and Gouhier et al. (2020), to test our atmospheric hazard maps, we compared them with the VAGs provided by the Toulouse VAAC during the course of the three lava fountains. In this way, we assess our results against independent information about the position of the volcanic cloud. Moreover, in an operational setting, VAGs are a good metric to test our forecasts since they are supplied in near real time during the eruptive events. It is worth noticing that, if available, satellite retrievals of ash columnar content could be directly used for the assessment of the numerical results. For illustrative purposes, for each episode, we selected the VAG closest in time to one of our forecasts. Results are shown in Fig. 18, where the contour lines of the VAG (red lines) are superimposed on the contour lines provided by our hazard maps showing the 5% probability to exceed 0.2 mg/m 3 of tephra concentration and indicating the extent of the volcanic cloud at different atmospheric intervals (black lines). For each episode, we also indicate the vertical extent of the volcanic cloud as supplied by the VAG. This vertical extent is provided in terms of Flight Levels (FL), where FL 330 corresponds to 10050 m, FL 160 and FL 050 to 4900 m and 1500 m respectively, and FL 360 to 10950 m. For the 16 February lava fountain, panel (a) of Fig. 18 shows the comparison between the extent of the volcanic cloud at 18:00 as supplied by the VAG, and the extent of the cloud provided by our forecasts at 18:10. The same is presented in panel (b), where the cloud produced by the 28 February episode is shown at 12:00 from the VAG and at 11:50 from the numerical forecast. Finally, panel (c) shows the cloud of the 4 March eruption as coming from the VAG at 10:55 and from the workflow at 10:20. Overall, a good match between the numerical forecasts and the clouds provided by the VAACs emerges for the three cases. This is true for both the horizontal and the vertical extensions of the clouds. In particular, we can observe that the episode that produced the highest volcanic cloud is the 4 March 2021 eruption. This is confirmed by both our numerical forecast (the 10 − 20 km level of panel (c) is contaminated by ash) and the VAAC observations that report a maximum height of the cloud of 10950 m. For the 28 February 2021 case (panel (b)), the contour supplied by the VAAC is slightly shifted with respect to the position of the cloud supplied by our workflow. It is worth noticing that VAGs are produced from information coming from satellite-and ground-based instruments and such information are anyway subject to uncertainties.

Discussions
The results presented so far highlight that the automated workflow described in this paper allows us to quantify, in probabilistic terms and accounting for some main sources of uncertainty, the hazard due to tephra fallout and atmospheric dispersal in a reliable way, considerably improving the forecasts with respect to procedures which do not consider the real-time update of the eruptive conditions.
The run times that we achieved are suitable for realtime applications and comparable to the reaction times of the VAACs, which have 20 min since the eruption notification to supply the first VAA/VAG (Gouhier et al. Fig. 14 Sketch of the forecast cycles performed in real time during the 4 March 2021 lava fountain episode. In this case, five VONAs were issued during the eruption, and for each VONA a forecast cycle was performed (from Initial cycle to Update 4) 2020). Indeed, the production of the results of a forecast cycle (i.e., forecasts at +1h, +2h, +3h, +6h, +9h, and +12h from the issue of a red VONA) requires a computational time ranging from 45 to 90 min. In particular, the production of the forecasts at +1h, +2h, and +3h takes from 5 to 10 min each. This means that we have the forecast at +1h 5/10 min after the release of the red VONA, while the forecasts at +2h and +3h are produced 10/20 and 15/30 min after the release of the VONA, respectively. Each of the following forecasts (i.e., +6h, +9h, and +12h) is produced in 20/30 min.
Overall, the numerical results agree well with the observations of the deposits and the information on the volcanic clouds supplied by the VAAC. However, some critical points have to be highlighted and discussed.
The first aspect is whether or not the time interval between the first red VONA and the final orange VONA well represents the real duration of tephra emission from the vent. During volcanic unrest, the volcanologist on duty at the observatory decides if and when to release the VONAs depending on the atmospheric and volcanological conditions. Since VONAs are used for aviation hazard assessment, the volcanologists usually tend to be cautious with the release of a final orange VONA in case the poor visibility of the summit area poses some doubts on the eventuality that tephra emission is ongoing. This can likely lead to overestimating the duration and intensity of tephra fallout from paroxysmal eruptions. An example of such a possibility is represented by the 16 February 2021 lava fountain. This episode took place in the late afternoon and thus with a very low visibility of the summit area of Mt. Etna although there was no bad weather conditions around the volcano. For this reason, after the first red VONA was issued at t Red 0 =16:31, the volcanologist on duty at the INGV-OE waited until the morning of 17 February (when daylight conditions and the visibility of the summit area were good) to release the orange VONA stating the end of tephra emission. Thus, according to the VONAs, about 13 h of tephra emission characterized the 16 February episode. Such dynamics disagrees with the eruption timeline resulting from other sources like the tremor-derived data, suggesting a total duration of the episode of 1.5 h (Andronico et al. 2021a). In Fig. 17 panel (b), we show that, assuming the duration of the episode inferred from the tremor, the simulated and observed ground loads are in much better agreement than assuming 13 h of emission as suggested by the VONAs. In the VONA case, the forecast workflow strongly overestimates the ground loads and then the hazard. A different situation emerges for the other case studies, where the VONA-derived duration and dynamics are in agreement with what is suggested by the tremor. In particular, the 28 February 2021 episode started in the morning (t Red 0 =08:11) and the excellent visibility of the summit craters allowed the volcanologist on duty to well constrain the dynamics of the episode, including the development of the eruption column. Indeed, 6 h of tephra emission was recorded for this episode, with the most intense phase of the eruption (the observed column height was 9 km a.s.l.) characterizing the first 2.5 h (from 08:11 to 10:47), while a weaker activity (column height 4.5 km) preceded the release of the orange VONA at t Or 3 =14:26. These dynamics agree with the tremor-derived data that show a total duration of the episode of 1.7 h, with the most intense phase at 8:30 (a time consistent with the emission of the second red VONA reporting a column height of 9 km a.s.l.) (Andronico et al. 2021a). Similarly, the progression of the volcanic column of the 4 March 2021 episode can be considered well constrained in terms of VONAs. Indeed, despite the episode started during the night (t Red 0 =02:28), the most intense phase of the eruption took place in the morning, with an eruption column of 11 km a.s.l. starting from t Red 3 =08:17. About 2 h later (at 10:16), an orange VONA reporting the end of tephra emission was issued. The tremor-derived duration of this episode is 3.8 h, with the peak in intensity at 9:20. These dynamics fit well with what emerges from the VONAs, where the most intense phase of the paroxysmal activity is from 8:17 to 10:16. A second point to discuss is the importance of evaluating the height of the eruption column in real time during an episode and of assimilating such information in the numerical workflow. Indeed, as it emerges from the Sensitivity Analysis, the height of the eruption column is the input parameter with the greatest control on the ground load, and thus on the reliability of the hazard assessment. According to the findings of Figs. 6 and 7, completely different hazard scenarios (both on the ground and in the air) derive from different assumptions on the height of the eruption column. On the ground, the 95th percentile of the tephra load in Catania ranges from 0.002 to 27.8 km/m 2 depending on the height of the eruption column assumed to produce the forecasts. In the air, ash contamination affects different atmospheric levels depending on the assumed column height.
Overall, in case the poor visibility of the eruptive episode anticipates (or delays) the issue of the first (or last) VONA and/or the height of the eruption column can not be estimated during an episode, the resulting hazard assessment likely tends to overestimate the hazard due to tephra fallout and volcanic cloud transport and dispersal. However, for hazard assessment analyses of this type, overestimating the likely expected or ongoing scenario due to the existence of main uncertainty sources is certainly preferable compared to underestimating it.
With the aim of producing reliable forecasts, the assimilation of observational data could involve not only the height of the eruption column (as it is already done in the present workflow), but also the position and concentration of the volcanic ash cloud as detected from satellite-based sensors. As shown by Pardini et al. (2020) and Mingari et al. (2022), numerical simulations of volcanic clouds can greatly benefit from the information coming from space. Together with data assimilation, machine learning is showing an increasing use in volcanology. Some recent papers show the application of machine learning algorithms to forecast the onset of volcanic eruptions and the changes in eruptive behavior (Watson 2020;Oláh and Tanaka 2022). In our case, data coming from previous paroxysms (seismic, deformation, infrasonic, thermal, geochemical data) could be used to create predictive models of column regime. Such models could then be used to initialize TTDMs allowing the production of forecasts in advance compared to what is done today.
We highlight the importance of atmospheric wind (both direction and intensity) on the evaluation of volcanic hazard. This is evident especially for the more proximal area, where wind direction has the same role of column height in determining the value of the ground load (Fig. 5). This is because the atmospheric diffusion has no time to spread the coarser particles that deposit quickly near the vent as soon as they leave the column. Thus, wind direction strongly controls where these particles deposit and small errors in For the three eruptions, we considered the ground load computed at the end of the deposition process, which means 12 h after the issuing of the last orange VONA. For each episode, colored diamonds indicate the 50th percentile of the ground load, while error bars indicate the 5th (lower limit) and the 95th (upper limit) percentiles. For the 16 February 2021 episode, blue open diamonds indicate the load computed assuming a duration of the episode consistent with the tremor-derived duration (at 18:10 on the 16 February 2021). Dashed and dot lines correspond to 5 (1/5) and 10 (1/10) factors, respectively its estimation greatly affect the hazard assessment near the vent. For this reason, considering the uncertainties in wind direction and intensity for the production of the ensembles of simulations allows us to mitigate this effect on hazard evaluation. We also highlight that our simulations consider only the dry deposition of the tephra particles, without including wet deposition due to rainfall. This latter aspect should be considered for further investigations. Indeed, it has been shown that, for volcanoes located in damp tropical regions, wet deposition, enhancing the removal of tephra from the atmosphere, can significantly affect the forecasts, especially during periods of high rainfalls (Dare et al. 2016).
As a final consideration, we point out the importance of standardizing the way in which the VONAs are written, especially in the light of automated workflows as the one that we propose. For a robust workflow, it is fundamental that the automated reading of the VONA is performed without errors and with the confidence that all the data are properly retrieved, especially those about the height of the eruption column. Recent work has been done in this direction under the European volcanological community networking program EUROVOLC.

Conclusion
We have presented an automated workflow for the realtime hazard assessment for tephra dispersal and fallout at Mt. Etna. The workflow is driven by the Volcano Observatory Notice for Aviation (VONAs) issued by the Istituto Nazionale di Geofisica e Vulcanologia-Osservatorio Etneo during the course of paroxysmal eruptions. Numerical results are processed to produce probabilistic hazard maps and tables showing the quantities of interest (i.e., ground load, atmospheric concentration) able to consider some of the main uncertainties affecting the volcanic phenomenon such as the column height, wind pattern, and TGSD of the pyroclastic mixture. Such results are produced and can be disseminated in real time during the course of the eruptive activity.
We successfully tested the workflow against three lava fountains of the February-April 2021 paroxysmal sequence at Mt. Etna. These are the episodes that occurred on 16 February, 28 February, and 4 March. We found that the simulated ground load is in good agreement with that derived from the analysis of the erupted products. Fig. 18 Comparison between the numerical forecasts and the extent of the volcanic clouds as supplied by the VAG produced by the Toulouse VAAC during the course of the three eruptive episodes. Panel (a) refers to the 16 February lava fountain, while panels (b) and (c) refer to the fountains that occurred on 28 February and 4 March 2021, respectively. In each panel, the black contour lines are from our numerically predicted ash cloud and indicate the 5% probability of exceeding 0.2 mg/m 3 in the three atmospheric intervals: 0 − 5 km, 5 − 10 km, 10 − 20 km. Red polygons indicate the horizontal extension of the volcanic cloud as supplied by the VAG. In each plot, we also indicate the vertical extension of the cloud from the VAG. SFC indicates the surface level, while the top height of the cloud is given as Flight Level. Note that the maximum height predicted by the VAAC for the 16 and 28 February clouds is ≤ 10 km. For this reason, VAAC polygons are not drawn in the interval 10 − 20 km of panels (a) and (b) Similarly, good results were found from the comparison of the simulated volcanic clouds with those observed and estimated by the Toulouse VAAC during the three episodes.
These results confirm the reliability of the present workflow for the real-time evaluation of the volcanic hazard associated with tephra dispersal and fallout at Mt. Etna. The most accurate results were obtained when the sequence of VONAs issued for a specific episode was able to well constrain the duration of tephra emission and the height of the eruption column. On the contrary, the delayed release of the report or the lack of column height observation result in the overestimation of the mass load of tephra on the ground and the extension of the volcanic cloud in the atmosphere. Such eventuality is well justified in terms of volcanic risk assessment, since it prevents underestimating the hazard of poorly observed eruptions.
We set and tested the forecast procedure at Mt. Etna, but the same workflow can be easily applied to any volcano in the world if the eruptive parameters are supplied in real time by monitoring agencies, as well as the meteo data is available.

Appendix: Calculation of the umbrella area source
The use of the PLUME-MoM-TSM allows us to simulate the spreading of the umbrella cloud by solving a set of shallow-water equations describing the lateral spread of the mixture above the neutral buoyancy height. The umbrella model is fed by source terms derived from the neutral buoyancy level horizontal section of the column and the characteristics of the atmosphere at that height (de' Michieli Vitturi and Pardini 2021). At the neutral buoyancy height, the horizontal section of the column is described by the downwind distance from the vent of the center of the cross section, denoted with d 0 , and its radius, denoted with r 0 (see Fig. 19). It is worth noting that both d 0 and r 0 are computed by the column model in a few seconds. When the umbrella reaches a steady upwind spreading, we compute a new circular area which best fits the upwind part of the umbrella. This new circular area is the source from which the finest particles reaching the neutral buoyancy height are Fig. 19 Oblique, top and lateral view of the plume. In the top view, the small red circle highlights the cross-sectional area of the plume at the neutral buoyancy level, while the big circle defines the area source computed by the umbrella model. The radii of the two circular areas (r 0 and r steady ) and their downwind distances from the vent (d 0 and d steady ) are indicated injected into the atmosphere. The downwind distance of the center of this area from the vent is denoted with d steady , and its radius with r steady . We indicate with d = d steady − d 0 the change in the horizontal downwind distance (Fig. 19).
The calculation of d steady and r steady from the umbrella model is computationally expensive and not suited for real-time applications. For this reason, we searched for a functional relationship between the geometric characteristic of the umbrella source area and the characteristic of the eruption column and the atmospheric wind field at the neutral buoyancy level. In this way, the source area is computed directly from the outcomes of the column model (model which is solved in a few seconds), without running the umbrella model. To do this, we fit the results of 1000 simulations performed by varying the mass flow rate at the vent in the range [10 4 ; 10 8 ] kg s −1 , the wind at the tropopause in the range [15; 100] m s −1 , and the specific humidity at sea level in the range [0.004; 0.0018] kg kg −1 .
At the steady state (i.e., when the upwind spreading has reached the stationary), the fitting equation for d is: where v[ms −1 ] is the module of the horizontal wind velocity at the neutral buoyancy level and vf r[m 3 s −1 ] is the volumetric flow rate at the crossectional area of the plume at the neutral buoyancy level (i.e., the volumetric flow rate feeding the umbrella region). The coefficients a, b, and c were found through a least-squares minimization procedure and are equal to 4.02 · 10 −9 , −0.7456, and 14.07, respectively.
Knowing d from Eq. 4 and d 0 from the results of the column model, the new downwind distance of the plume center from the vent is: and it is evaluated in the direction of the atmospheric wind at the neutral buoyancy height. For the new radius (r steady ), a similar fitting equation based on v and vf r was constructed: r steady = Av B log 10 (vf r) C , where the least-squares minimization gives A = 1.89 · 10 −8 , B = −0.8856, and C = 13.9. Finally, we also fit the time necessary to reach the steady upwind spreading as follows: t steady = αv β log 10 (vf r) γ , where again the coefficients α, β, and γ derived from a least-squares minimization procedure and are equal to 1.9 · 10 −5 , −1.57, and 10.9, respectively. For eruptive episodes lasting a time t < t steady (i.e., the eruption column is not longer sustained at time t steady ), the definition of the circular area of the umbrella region should be based on the downwind distance from the vent and on the radius computed at time t (d(t) and r(t)) and not at time t steady . To do this, we defined the function varying in the range [0,1]: which well represents the variation of the radius of the umbrella cloud and the upwind spreading distance with time. Thus, to compute d(t) and r(t), we used Eq. 8 to lineally interpolate between time 0 (the beginning of the umbrella spreading for which d 0 and r 0 are known from the plume model) and t steady (for which d steady and r steady are known from Eqs. 4, 5, and 6): Supplementary information The online version contains supplementary material available at https://dx.doi.org/10.1007/s00445-022-01614-z.