Effects of Variable Eruption Source Parameters on Volcanic Plume Transport: Example of the 23 November 2013 Paroxysm of Etna

The purpose of the present paper is to investigate the effects of variable eruption source parameters on volcanic plume transport in the Mediterranean basin after the paroxysm of Mount Etna on 23 November 2013. This paroxysm was characterized by a north-east transport of ash and gas, caused by a low-pressure system in northern Italy. It is evaluated here in a joint approach considering the WRF-Chem model configured with eruption source parameters (ESPs) obtained elaborating the raw data from the VOLDORAD-2B (V2B) Doppler radar system. This allows the inclusion of the transient and fluctuating nature of the volcanic emissions to accurately model the atmospheric dispersion of ash and gas. Two model configurations were considered: the first with the climax values for the ESP and the second with the time-varying ESP according to the time profiles of the mass eruption rate recorded by the V2B radar. It is demonstrated that the second configuration produces a considerably better comparison with satellite retrievals from different sensors platforms (Ozone Mapping and Profiler Suite, Meteosat Second-Generation Spinning Enhanced Visible and Infrared Imager, and Visible Infrared Imaging Radiometer Suite). In the context of volcanic ash transport dispersion modeling, our results indicate the need for (i) the use of time-varying ESP, and (ii) a joint approach between an online coupled chemical transport model like WRF-Chem and direct near-source measurements, such as those carried out by the V2B Doppler radar system.


Introduction
Explosive volcanic eruptions inject large amounts of particles, aerosols, and gases into the atmosphere, where they are carried away for several thousand kilometers and remain suspended in the air for several days or months. Once in the atmosphere, these ash particles, typically ranging from a few millimeters to a few micrometers, are usually involved in physico-chemical processes, influencing weather and climate [1]. While the transport of ash particles mainly concerns airport security and flight safety and the sedimentation of ash particles affects ground infrastructures and human health, the oxidation of the volcanogenic SO 2 leads to the formation of sulfate particles, through chemical and Lagrangian particle trajectories models, we have utilized HYSPLIT [22], NAME [23], and Puff [24] in various contexts.
Likely, the Eyjafjallajökull eruption (15 April 2010) is the most studied recent event due to its strong impact on air traffic and economy in Europe. In this context, Plu et al. [25] simulated this eruption with the MOCAGE-CTM (Météo-France) and hourly changing MER from FPlume [6]. For the same eruption, Devenish et al. [26] and Webster et al. [27] both used the lagrangian model NAME while Folch et al. [28] used the FALL3D model. More recently, the CHIMERE model [29] coupled offline with WRF has been utilized to study the 18 March 2012 Mount Etna eruption [30].
The online coupling strategy is currently developed by using the WRF-Chem model [31][32][33] in which a specialized emission routine for the treatment of volcanic particles and gas (SO 2 ) has been implemented in its online chemistry-meteorology coupling framework [8]. This also permits inclusion of the direct and indirect feedback, together with other peculiar processes on the ash particles, like aggregation [33], microphysical nucleation [34], and the calculation of the optical properties (extinction coefficients and aerosol optical depth). Moreover, the "online CTM strategy" has been well suited for the utilization of time-dependent (transient) ESP like those from the V2B system. Nowadays, the WRF-Chem and radar combination, considering the transient and fluctuating nature of the volcanic emissions, represents an innovative advance to the current VATD modeling systems.
In this work, the source parameter evolution obtained from the measurements of the V2B Doppler radar system were assimilated inside the WRF-Chem input framework. The V2B raw data were pre-processed together with the other specific volcanic parameters to define a user-defined volcanic input file to allow multi-eruptions with a single run. The violent paroxysm of Mount Etna on 23 November 2013 was considered as a test case. Satellite retrievals by the Ozone Mapping and Profiler Suite (OMPS) sensor and the Visible Infrared Imaging Radiometer Suite (VIIRS) sensor both onboard the Suomi National Polar-orbiting Partnership (NPP) spacecraft, and the Meteosat Second Generation (MSG) Spinning Enhanced Visible and Infrared Imager (SEVIRI) images data were used to validate the outputs of the different setup configurations.
This work represents the preliminary step towards the forecast simulations of meteorology and ash dispersion of Mount Etna. This would be achievable in near real-time, assimilating into WRF-Chem model the input data from the VOLDORAD-2B monitoring system, aside from the end-member eruption scenarios used nowadays (one-weak long-lived emission similar to 2002-2003 and one strong short-lived eruption similar to 5 January 1990).
This new strategy (i.e., WRF-Chem and radar-based ESP) has three main objectives trying to respond to some of the new challenges of modern volcanology, namely (i) utilization of realistic time-dependent ESP and their assimilation in an online coupled CTM, (ii) provision of near-real time maps of volcanic aerosols that may help in identifying potential risk regions for flight security, and (iii) analysis of the effects of volcanic aerosols in atmospheric (longwave and shortwave) radiation.
The paper is organized as follows: Section 2 is dedicated to the description of the material and methods utilized in this work, in particular the suite of the remote sensing data and the configurations of the modeling system; in Section 3, we present all the simulations results and their comparison with experimental data; and in Section 4, we report the conclusions of the present study.

The 23 November 2013 Etna Paroxysm: Eruption Timeline and Synoptic Conditions
On the morning of 23 November 2013, Etna's New SE Crater produced the 42nd paroxysmal episode since January 2011. The ash plume top reached up to 12.6 km a.s.l. [35]. The strong wind spread the tephra plume toward NE, generating abundant fallout of bombs up to several kilometers and dispersing fine lapilli all along the Ionian coast of Sicily. Thick fallout caused serious trouble to infrastructures including the highway about 10 km east. Ash fallout was reported in southern Calabria and Puglia regions in southern Italy, and as far as 400 km from the volcano. From thermal images analysis, Bonaccorso et al. [36] inferred a time-averaged-discharge-rate of~360 m 3 s −1 , i.e.,~1.6 × 10 6 m 3 of dense-rock equivalent magma volume erupted in just 45 min.
The maps of the geopotential height (NOAA-Global Forecast System: https://www. ncei.noaa.gov/products/weather-climate-models/global-forecast (accessed on 18 February 2021) relating to the isobaric surface at 500 hPa of the day of 23 November 2013 show a large depression circulation that affects much of the central-western Mediterranean (not shown). The pressure minimum off-shore southern France tends to isolate itself in the middle Tyrrhenian Sea. This synoptic framework, in addition to bringing extensive conditions of atmospheric instability over much of Italy, with precipitation also of a stormy nature especially in the Tyrrhenian areas, favored the genesis of intense winds at high altitude. In particular, as shown by the wind speed maps at 300 hPa at the main synoptic hours (Figure 1a at 00:00 UTC, Figure 1b at 06:00 UTC, Figure 1c at 12:00 UTC, and Figure 1d at 18:00 UTC), thanks to the transit of a secondary branch of the Jet Current, Sicily was affected by winds at high altitude, locally even higher than 50 m s −1 (180 km h −1 ), coming mainly from the south-western quadrants. The presence of a low-pressure system in the middle of Tyrrhenian Sea, favored a flux of cold air at high altitude over the Mediterranean basin from the Scandinavian Peninsula (not shown in Figure 1). In particular, masses of cold air at the isobaric surface of 850 hPa (about 1500 m above sea level) mostly affected the western sector of Italy at 00:00 UTC. In the following hours of 23 November, the colder air masses affected the entire national territory, extending to a large part of the Adriatic Sea.  [37]). Using a fixed-pointing antenna beam, this 23.5 cm-wavelength pulsed radar sounds 13 atmospheric volumes, 120 m-deep, right above the summit craters. The echo power and radial velocities (level 1 data) are computed from the raw Doppler spectra recorded at intervals of 0.22 s. These unique near-source measurements of the emitted tephra allow maximum ejection velocities at the base of the eruptive column and mass eruption rates (level 2 data) to be determined in real-time [14][15][16]38]. Owing to its capacity for implementation in real-time, we used the methodology of Freret-Lorgeril et al. [15] to compute MER from the product between the echo power and radial velocities measured in the two beam volumes above the New SE Crater, from which the 23 November 2013 paroxysm originated. This product has been shown to be proportional to the MER of tephra passing through the probed volume (MER proxy) and correlates to the observed plume heights for many paroxysms at Etna. Therefore, the VOLDORAD-2B MER proxies could be calibrated with MER derived from the plume model of Degruyter and Bonadonna [39], matching the same plume heights and including wind vertical profiles measured during eruptions.

Suomi NPP (OMPS/VIIRS) Data
The Ozone Mapping and Profiler Suite (OMPS) is onboard the joint National Aeronautics and Space Administration/National Oceanic and Atmospheric Administration (NASA/NOAA) Suomi National Polar-orbiting Partnership (Suomi NPP) satellite, launched in October 2011. In particular, the OMPS-NPP L2 NM Aerosol Index swath orbital product provides aerosol index values from the OMPS Nadir-Mapper (NM) instrument. This is actually the official NASA aerosol index product, which replaces the aerosol index found in the OMPS-NPP L2 NM Total Ozone product. The aerosol index is derived from normalized radiances using two wavelength pairs at 340 and 378.5 nm in the ultra-violet spectrum. The dataset is accessible through the NASA Earthdata Search web portal (https://search.earthdata.nasa.gov/search (accessed on 5 May 2021).
The Visible Infrared Imaging Radiometer Suite (VIIRS) is onboard the Suomi NPP and the Joint Polar Satellite System (JPSS) spacecrafts. The VIIRS instrument has 22-spectral bands, covering the spectrum between 0.412 and 12.01 µm, including 16 moderate resolution bands (M-bands) with a spatial resolution of 750 m at nadir, and 5 imaging resolution bands (I-bands) with a 375 m resolution at nadir. The brightness temperature was downloaded from the NOAA-CLASS web portal at https://www.avl.class.noaa.gov/saa/ products/welcome (accessed on 6 July 2021) by selecting the VIIRS data records from the Suomi NPP products. The brightness temperature of the thermal infrared channels has been widely used in dust storm monitoring [40] but also to detect volcanic plumes [41].

MSG SEVIRI Images Data
The Spinning Enhanced Visible and InfraRed Imager (SEVIRI: https://www.eumetsat. int/seviri (accessed on 17 May 2021)) is a 12-channel sensor onboard Meteosat Second Generation (MSG) spacecraft. It is a geostationary, meteorological satellite developed by the European Space Agency in close co-operation with the European Organisation for the Exploitation of Meteorological Satellites. Product images are available at a frequency of 1 image/15 min with a pixel resolution of 3 × 3 km at nadir. The HOTVOLC observing system ( [42]; https://hotvolc.opgc.fr (accessed on 17 May 2021) is a web-based satellitedata-driven reporting system dedicated to the observation of volcanic products including in real-time. HOTVOLC is a Web-GIS (Geographic Information System) volcano moni-toring system developed at the Observatoire de Physique du Globe de Clermont-Ferrand (France). It uses on-site automated ingestion of MSG-SEVIRI geostationary data for direct dissemination on a full web-GIS interface of ready-to-use products characterizing volcanic ash, sulphur dioxide, and lava flow emissions from spectral band combinations. Satellite products are delivered in the form of geo-referenced images (geotiff) tiled on a background map and time series (csv) through an open-access database. Although several elaborated ash products (ash plume altitude, mass, contour) have been made available recently from combinations of thermal infrared channel radiance data, for the 23 November 2013 ash and SO 2 plume, we were able to use only ash RGB composite images from the HOTVOLC archive data combining IR8.7, IR10.8, and IR12.0 channels.

Multi-Satellite Volcanic Sulfur Dioxide Long-Term Global Database
The input parameters for SO 2

Numerical Setup
A new specialized module has been implemented inside the Weather Research and Forecasting (WRF) with Chemistry model (WRF-Chem, [31]) that allows simulation of the emission, transport, and settling of aerosols and gases released during volcanic eruptions [8]. Among the WRF-Chem specialized options that are available for the description of the transport of volcanic ash and gas (SO 2 ), we used the inert tracer options that utilize 10 ash variables and SO 2 .
The numerical grid, depicted in Figure   About the setting of the WRF physical parameters, the surface layer module corresponds to the Mellor-Yamada-Nakanishi and Niino (MYNN, [43]) scheme (sf_sfclay_physics = 5), the Rapid Update Cycle Land Surface Model [44] was used to represent the land surface interactions (sf_surface_physics = 3) and the MYNN 2.5 level turbulent kinetic energy parameterization was used to describe the planetary boundary layer parameterization (bl_pbl_physics = 5). Short-and long-wave radiation effects were modeled using the Rapid Radiative Transfer Model (RRTMG, [45]) for both short-and long-wave (ra_sw_physics = ra_lw_physics = 4); the two-moment cloud model microphysics scheme of Morrison et al. [46] was used for the treatment of the microphysical processes (mp_physics = 10).
This configuration has been already successfully tested by Rizza et al. [47] and well suited for the study of transport of natural aerosols in the Mediterranean basin.

Eruption Source Parameters
The definition of ESP together with the other specific volcanic parameters, namely the volcano location in the grid domain, start and end of each single paroxysm (UTC time), and TGSD, was calculated at the pre-processing level utilizing an NCAR Command Language script. The resulting ASCII file contained all the volcanic information needed at runtime by WRF-Chem. Future releases will consider a unique pre-processing environment also including the elaboration of the raw V2B radar data.  The value of MER [kg s −1 ] was computed from the level 1 data using the calibration method of Freret-Lorgeril et al. [15]. A total erupted mass (TEM) of 4.7 10 9 kg was computed from the integration of MER over the duration of the tephra emission totaling 193 min. The climax was found to release 83% of the TEM and to have an MER 8.5 times higher (3.4 × 10 6 kg s −1 ) than that of the whole paroxysm (4.0 × 10 5 kg s −1 ). It may be notified that the average peak values of the exit velocity (approximately 150 m s −1 ) and MER (approximately 10 7 ) in the climax agree well with an independent analysis performed by Montopoli [48] on the same case study.
The plume height data were measured from the ECV visible camera of Istituto Nazionale di Geofisica e Vulcanologia, Osservatorio Etneo [49] above the vent, so that plume altitude a.s.l. can be obtained by adding about 3.1 km [50].
The V2B raw data were pre-elaborated to be compatible with the WRF-Chem system in particular, they were averaged on a 40 s interval to be incorporated by the numerical model at each time step. The corresponding time series of the plume height ( Figure 4a) and MER (Figure 4b) are depicted in Figure 4. To establish the initial vertical distribution of the erupted mass, it was assumed that the volcanic plume is described by an idealized umbrella shape [51]. In this context, the erupted mass is distributed in the vertical according to the so-called umbrella cloud model. According to the default WRF-Chem configuration [8], it was assumed that 75% of the erupted mass is constrained in the umbrella cloud and 25% below, with a linear distribution from the umbrella base to the vent (Figure 1 of Stuefer et al. [8]).
The elaboration of plume height and V2B data reported in Figure 4a,b suggests that in the paroxysmal phase (09:55-10:14 UTC), the eruption reached an altitude of 11 km (a.s.l.) and then rapidly decreased to around 5 km (a.s.l.) after 10:30 UTC, causing different cloud layers being transported at different altitudes. These results are in agreement (within 1 km) with the analysis of Corradini et al. [35] on the same Etna eruption.
In WRF-Chem, as described by Stuefer et al. [8], the TGSD is distributed between 10 bins of aerosol particles with a diameter size range starting from 2 mm down to less than 3.9 µm. Table 1 gives the selected particle size bins, which were associated with the WRF-Chem variables named vash_1 to vash_10, and the corresponding mass fraction percentage.
The particle size diameter (mm) was computed according to the d = 2 −φ function and the corresponding E1 distribution among the 10 bins was recalculated elaborating the eruptive data obtained by Poret et al. [4] for the 23 November 2013 Etna eruption.  Table 2 shows the ESP and chemistry setup for the two configurations, namely run1 in which ESP are kept constant for the whole paroxysm event (07:13-10:30 UTC) and the run2 configuration that uses time-varying average ESP according to the elaboration of the V2B Doppler radar data reported in Figure 3. The setting chem_opt = 402 corresponds to the inert tracer options that utilizes 10 ash variables (cf. Table 1) and SO 2 . In this context, the model sensitivity to the use of transient vs. constant ESP may be evaluated. Table 2. Eruption source parameters and WRF-Chem chemistry option. ECV and V2B indicate time-varying MER retrieved from the VOLDORAD-2B radar ( Figure 3) and plume heights retrieved from a video camera (Figure 4)

Geopotential ERA5 vs. WRF
In order to analyze the synoptic conditions in the days following the eruption of Etna, we considered the ERA5 reanalysis [52], which have 137 vertical levels and about 30 km of horizontal resolution. Geopotential fields were obtained by downloading ERA5 hourly data on pressure levels (available from 1979 to present) in netCDF format (https://www.unidata.ucar.edu/software/netcdf/) from the Copernicus Climate datastore (https://cds.climate.copernicus.eu/#!/home (accessed on 13 February 2021).
To compare ERA5 with the WRF-Chem model output fields, we focused on the geopotential height at 500 and 200 hPa. These levels were chosen since they are representative of the upper tropospheric conditions, which determine the long-range transport of the volcanic plumes.
Panels (a, c, e) of Figures 5 and 6 represent ERA5 reanalysis of the geopotential at 200 and 500 hPa (in km), respectively, while panels (b, d, f) represent the WRF-Chem output of the same variable. In each figure, the top, middle, and bottom rows represent snapshots at 12:00 UTC on November 23, 24, and 25, respectively.  As depicted by Figure 6a-d the synoptic configuration at 500 hPa is quite similar; however, at these lower levels, the circulation in the region invested by the volcanic plume (36-40 • latitude north and 15-22 • longitude east) is slightly different and more oriented northward. These synoptic conditions highlight the different directions taken by the volcanic plume at approximately 12 and 5 km altitude (a.s.l.) as it will be depicted in the next paragraphs.
It is quite evident how the model (panels b, d, f) is able to reproduce the geopotential fields very well at both pressure levels, equally in terms of configuration and intensity.

Comparison for SO 2
The Goddard Space Flight Center (GSFC/NASA), through its Global Sulfur Dioxide Monitoring Home Page (https://so2.gsfc.nasa.gov (accessed on 5 May 2021)), makes data available for downloading, in particular SO 2 eruption alerts, SO 2 near real-time images, and finally archived daily OMI/OMPS/TROPOMI images.
The image reported in Figure 7 is relative to the eruption of 23 November 2013 retrieved from the OMPS sensor onboard the Suomi NPP satellite, which passed in its closest point to Sicily at 12:19 on the same day. This image is relative to the plume distribution about two hours after the eruption climax occurred around 10:00 UTC. As a consequence of the synoptic circulation, described in Section 3.1, the columnar SO 2 distribution after the main eruption is dislocated from Sicily to the Ionian Sea and moving in the north-east direction. The modeled hourly map at 12:00 UTC (same day) is reported in Figure 8a,b for run1 and run2, respectively. Two important points should be remarked: (i) the volcanic plume at 12:00 UTC for run1 is fragmented into two parts with the farthest portion already passed over the Albanian coast and (ii) the spatial pattern relative to run2 is very similar to OMPS retrievals at the same hour. This means that according to the run1 configuration, the injection height at 11 km (a.s.l.) for the whole duration of the paroxysm is a wrong assumption that leads to an unrealistic spatial distribution. This is caused by the presence of the ash cloud at the highest tropospheric levels with considerable high wind speed. Contrastingly, the transient ESP under the run2 configuration seems to produce a remarkably closer spatial pattern, both in terms of plume extension and centers of mass.
Moreover, if we consider the columnar density of SO 2 , expressed in Dobson Units (DU), according to the simulation outputs, the SO 2 concentration of run2 results in approximately 2.0 DU, which is close the observed value of about 2 DU, while run1 outputs provide a lower value around 0.3 DU. Unfortunately, on the following days (24 and 25 November 2013), there are no satellite retrievals that may be utilized for a further comparison. The closer match of the SO 2 concentration and spatial pattern highlights the need for specifying the appropriate emission height based on time-varying ESP.

Comparison for Volcanic Ash
To validate the model predictions with respect to ash, the simulation results were compared with the observed data characterized by the Aerosol Index (AI) that was explicitly developed to detect the presence of aerosols in the air. In particular, the OMPS Aerosol Index layer specifies the presence of ultraviolet (UV)-absorbing particles in the air. In this context, it is useful for identifying and tracking the long-range transport of volcanic ash from volcanic eruptions, but it also may be utilized for other aerosol types, like desert dust, smoke from wildfires, or biomass burning. It is a science parameter of OMPS/National Polar, and derived from normalized radiances using two wavelength pairs at 340 and 378.5 nm. In Figure 9, AI is visualized from Level-2 Swath and 50 × 50 km resolution, over 1 h time coverage around 12:00 UTC on 23 November 2013. A modeled quantity that could be qualitatively compared with the OMPS/AI is the columnar mass density, which is obtained by considering in each (x, y) grid point the following columnar quantity [53]: where (zmin = 0 km, zmax = 20 km) are the vertical integration limits, vashT = ∑ 10 1 vash(i) represents the sum over all volcanic bins expressed in mixing ratio units (µg kg −1 ), and ρ air is the air density (kg m −3 ). The resulting columnar quantity (Vash col ) is expressed in µg m −2 and reported in Figure 10a,b in Log 10 scale.
Considering that the comparison between Aerosol Index observations ( Figure 9) and the modeled columnar mass density ( Figure 10) is only qualitative, being performed between two quantities with different units, the very close reproduction of the spatial pattern by the run2 configuration (Figure 10b) at 12:00 UTC on 23 November 2013 is remarkable. In contrast, the extension at the same hour of the modeled plume over Albania and beyond in the run1 simulation (Figure 10a) reveals the unrealistic presence of volcanic ash at a high tropospheric level as far as 20 km beyond the downwind extension limit of the observed aerosols. This means that the use of averaged parameters to describe the eruptive column may lead to very substantial errors in the time and space dispersion of volcanic ash and aerosols. In particular, the optimal comparison of observations and run2 simulation is mainly the result of the appropriate specification of the time-varying MER and transient injection height. It is worth mentioning that the use of the FALL3D model as depicted by Poret et al. [4] on the same paroxysm provided comparable results, as depicted by their Figure 11 (panels a, b).

Comparison with MSG-SEVIRI Data
In this section, we compare the MSG-SEVIRI RGB composite images from HOTVOLC with the outputs of the WRF-Chem model under the run2 configuration. Figure 11 represents the plume snapshots at 11:00 (upper row), 12:00 (central row), and 13:00 UTC (lower row) on 23 November 2013. Panels (a, c, e) denote the RGB satellite images while panels (b, d, f) display the WRF-Chem outputs from the run2 configuration. In panels (a, c, e), the ash aerosols are colored in blue, whereas ice and water hydrometeors are colored in red and light green, respectively. The plume is traveling quite fast, being driven by north-east winds at almost 120 km h −1 , reaching Albania at 13:00 UTC.
The RGB images show a narrow volcanic plume denoted by the blue filament inside the black dotted oval that is "headed" by a red/green spot representing the presence of water/ice hydrometeors transported with the plume at 11:00 UTC (Figure 11a). This red/green spot is moving eastward in the following hours at 12:00 and 13:00 UTC (Figure 11c,e) while the rest of the plume, composed by ash aerosols, is moving in the north-east direction.
This probably indicates that a portion of the ash plume, depending on the local relative humidity and temperature, is ice nucleated and transformed into non-precipitating hydrometeors. This requires a microphysical-aerosol coupler [54] to be further investigated. Once more, it is important to remark on the optimal performance of the WRF-Chem model using time-varying ESP in input (run2 configuration), as depicted by panels b, d, and f of Figure 11.

Ash Transport at Vertical Layers
To analyze the plume transport at different heights, we have split the plume into two parts. In Equation (1), the lower portion is calculated between zmin = 0 km and zmax = 8 km (a.s.l.) and the upper portion between zmin = 8 km (a.s.l.) and zmax = 12 km (a.s.l.). In this discussion, we only consider WRF-Chem outputs under the run2 configuration as it demonstrates superior performances compared with the run1 setup.
The upper panels of Figure 12a,b are the snapshots at 12:00 UTC while the lower panels (c, d) refer to 13:00 UTC of 23 November. Panels (a, c) are the vertical integrals between 0 and 8 km and panels (b, d) between 8 and 12 km (a.s.l.). In the lower portion Figure 12a,c, we can observe typical thin and narrow volcanic plume stretching very rapidly in the north-east direction and reaching Albania at 13:00 UTC, while the upper plume portion (panels b, d) that contains the largest quantity of the ash load (~10 8 -10 9 µg m −2 ) is transported in a more easterly direction. This behavior of the volcanic plume is confirmed by Figure A1 (Appendix A) showing the plume splitting around 14:00 UTC ( Figure A1b). This is compatible with the geopotentials at 200 and 500 hPa, shown in Section 3.1. In particular, the 200 and 500 hPa geopotentials at 12:00 UTC (Figures 5b and 6b, respectively) show small differences in wind direction at these two pressure heights.
These results are in agreement with the satellite data (MSG-SEVIRI) analysis of Poret et al. [4] highlighting the presence of a north-eastern volcanic ash cloud at about 6 km (a.s.l.) and a volcanic ice cloud at 11 km (a.s.l.). Finally, Figure 13 reports the brightness temperature from the channel I5 (BTI5) at 11.5 µm of the NPP-VIIRS sensor. It reveals the location of the ash and ice plumes, respectively. The analysis of BTI5 shows that the ash plume was at around 260-270 K at approximately 3 km (a.s.l.), while the ice plume was at ∼220 K at 9 km (a.s.l.). This is confirmed by Figure A2 in Appendix B, obtained by analyzing the brightness temperature difference between the two thermal infrared VIIRS channels B15 and B16 at 11.763 and 12.03 µm, respectively. The dotted region highlights the presence of ash with thermal difference DT < 0 [41] and ice with DT > 0. These results likely indicate that a portion of the ash plume, depending on the local relative humidity and temperature, can act as ice-nucleating particles [34]. This requires a microphysical-aerosol coupler [54,55] to be further investigated with the WRF-Chem model.

Conclusions
In this work, we investigated the effects of variable eruption source parameters on volcanic plume transport of the Mount Etna paroxysm of 23 November 2013. It represents a further application of the WRF-Chem model for Mount Etna following Rizza et al. [32]. Concerning the model setup, the raw data obtained from the V2B Doppler radar system were elaborated at the WRF-Chem pre-processing level. This allowed the inclusion of the transient and fluctuating nature of the volcanic emissions to accurately model atmospheric dispersion of ash/gas, pointing to the exact definition of the most critical parameters, such as paroxysm onset and end, column height, total grain size distribution, and mass eruption rate (MER).
This joint approach between an online coupled chemical transport model (CTM) like WRF-Chem and the data from the V2B Doppler radar system represents an innovative development of current volcanic ash transport and dispersion models as it considers the two-way feedback for all the processes (meteo, chemistry, and physics) at the time step level, in a context of time-varying eruption source parameters (ESPs). The meteorology is part of the system and consequently there is no need to extrapolate/interpolate external fields for both passive and active tracers.
The system was evaluated considering a Mount Etna eruption that has been largely investigated by the scientific community. This paroxysm was characterized by a north-east rapid transport of ash and gas, caused by a low-pressure system in northern Italy.
Two WRF-Chem configurations were evaluated: the first one with the climax values for the ESP, like in Rizza et al. [32], and the second with the time-varying ESP according to the MER time profiles recorded by the V2B radar. It was shown that the second configuration produces a considerably better comparison with satellite retrievals (Ozone Mapping and Profiler Suite, Meteosat Second-Generation Spinning Enhanced Visible and Infrared Imager, and Visible Infrared Imaging Radiometer Suite). To verify the role of the initial column height, three additional runs were performed and discussed in Appendix C. The relative analysis indicated that it is fundamental to describe with "great realism" MER and column height parameters during the climax phase of paroxysm. The "old style" characterization of the ESP based on eruptive average values does not adequately reproduce the transport of ash particles after the eruption. These results suggest the "main avenue" for future development of coupled CTM-radar systems, namely their full integration, at the level of the model time step, to describe all the processes involved with reactive volcanic tracers.
In this work and in the work by Rizza et al. [32], we demonstrated that the model WRF-Chem may be utilized for forecast simulations of meteorology and ash/gas dispersion of Mount Etna emissions. In addition, it may also be particularly useful to better understand its impact on the meteorological phenomena over the whole Mediterranean area.
Specific future upgrade of the WRF-Chem model aiming to further increase the accuracy of the predictions are currently underway, namely (i) a new subroutine to describe the initial column height from V2B-derived MER and 1-D column models, and (ii) a new gravitational settling scheme to consider large ash particles.
The principal suggestion for future research in the area of volcanic ash dispersion is to include a high-resolution imbedded model of volcanic plume. It is likely that the dedicated volcanic plume model will allow explicit simulation of particle coagulation and the transition between the vertical transport mode in the cloud and the predominantly horizontal mode associated with the flow resolved by the WRF-Chem.
All these considerations indicate that simulations would be achievable in near realtime in the future using VOLDORAD-2B monitoring data as input. Concluding this initial setup phase, the new system is scheduled to run in near real-time at the University of Messina (Department of Mathematical and Informatics Sciences) and will be able to provide hazard evaluations to governmental agencies immediately following Mount Etna's explosive eruptions.  Data Availability Statement: The data used for this work can be made available upon request and on a case-by-case basis. radar monitoring. This study used the VOLDORAD open-access data base (http://voldorad.opgc.fr/ (accessed on 17 May 2021)) of OPGC-Université Clermont Auvergne, with support from the EU EPOS, EUROVOLC (#731070) and MED-SUV programs and from the French SNOV.

Conflicts of Interest:
The authors declare no conflict of interest. Figure A1 shows four snapshots of the whole model columnar mass density that may be considered as representative of the plume transport during the first day (23 November). At 10:00 UTC ( Figure A1a), when the eruptive phase reaches its climax, the plume is still near Etna; at 14:00 UTC ( Figure A1b), it already reaches Albania and begins to split in two sectors, the upper portion going eastward and the lower portion in the north-east direction; at 16:00 UTC ( Figure A1c), these different directions for the lower and upper plumes are emphasized and finally at 21:00 UTC ( Figure A1d), only the lower northward portion remains inside the numerical domain. These results are in agreement with the analysis made by Poret et al. [4] with the HYSPLIT forward trajectories [56].  Figure A2 shows the normalized brightness temperature difference of the VIIRS thermal infrared channels. It was calculated as in Prata [41]:

Appendix B
where BT15 and BT16 are the brightness temperature for the channels 15 (10.763 µm) and 16 (12.013 µm), respectively, of the NPP VIIRS image M (700 m spatial resolution) channels. The dotted region highlights the presence of (i) a narrow ash plume with DT < 0 and (ii) ice and water droplets clouds (blue shaded). In this context, as described by Prata [41], the temperature difference built in this way may be utilized to discriminate between ash and ice clouds.

Appendix C
To directly verify the influence of column height, a series of additional runs were performed with the setup showed in Table 2 of Section 2.3.2. The three simulations differ in the way the column height is specified, namely a variable column height and constant MER (run3), a fixed MER and column height based on the climax values (run4), and finally a fixed MER and column height but based on average values (run5) during the whole duration of the paroxysm. Figure A3 shows that only a proper combination of MER and column height ( Figure A3b) reproduces the experimental spatial distribution depicted in Figure 9.
This series of additional runs indicate that it is fundamental to describe with "great realism" MER and column height parameters during the climax phase of paroxysm. The "old style" characterization of the ESP based on eruptive average values does not adequately reproduce the spatial pattern of ash after the eruption.