Dual-Wavelength Polarimetric Lidar Observations of the Volcanic Ash Cloud Produced during the 2016 Etna Eruption

Lidar observations are very useful to analyse dispersed volcanic clouds in the troposphere mainly because of their high range resolution, providing morphological as well as microphysical (size and mass) properties. In this work, we analyse the volcanic cloud of 18 May 2016 at Mt. Etna, in Italy, retrieved by polarimetric dual-wavelength Lidar measurements. We use the AMPLE (Aerosol Multi-Wavelength Polarization Lidar Experiment) system, located in Catania, about 25 km from the Etna summit craters, pointing at a thin volcanic cloud layer, clearly visible and dispersed from the summit craters at the altitude between 2 and 4 km and 6 and 7 km above the sea level. Both the backscattering and linear depolarization profiles at 355 nm (UV, ultraviolet) and 532 nm (VIS, visible) wavelengths, respectively, were obtained using different angles at 20◦, 30◦, 40◦ and 90◦. The proposed approach inverts the Lidar measurements with a physically based inversion methodology named Volcanic Ash Lidar Retrieval (VALR), based on Maximum-Likelihood (ML). VALRML can provide estimates of volcanic ash mean size and mass concentration at a resolution of few tens of meters. We also compared those results with two methods: Single-variate Regression (SR) and Multi-variate Regression (MR). SR uses the backscattering coefficient or backscattering and depolarization coefficients of one wavelength (UV or VIS in our cases). The MR method uses the backscattering coefficient of both wavelengths (UV and VIS). In absence of in situ airborne validation data, the discrepancy among the different retrieval techniques is estimated with respect to the VALR ML algorithm. The VALR ML analysis provides ash concentrations between about 0.1 μg/m3 and 1 mg/m3 and particle mean sizes of 0.1 μm and 6 μm, respectively. Results show that, for the SR method differences are less than <10%, using the backscattering coefficient only and backscattering and depolarization coefficients. Moreover, we find differences of 20–30% respect to VALR ML, considering well-known parametric retrieval methods. VALR algorithms show how a physics-based inversion approaches can effectively exploit the spectral-polarimetric Lidar AMPLE capability.


Introduction
One of the major hazards associated with volcanic explosive eruptions is the injection of volcanic ash into the atmosphere and its subsequent dispersion and deposition. Volcanic ash mainly affects aviation safety, although the impact could be reduced using real time observations and characterization of eruptive activity [1]. A variety of ground impacts also exist that change with distance from the volcanic vent (e.g., [2,3]).

Dual Wavelength Polarimetric Lidar
The capability of Lidar systems to detect the smallest particles (from 0.1 to 100 μm) in volcanic plumes as well as to reliably estimate the ash concentration depends on the instrumental characteristics and style of eruptive activity. Lidar measurements, in terms of particle backscattering coefficient and cross-polarization ratio profiles, can be carried out in regions where the Lidar signal is not extinguished due to the volcanic plume optical thickness even though some path corrections can be applied [6]. Figure 2a shows an image obtained by MODIS spectrometers (Moderate Imaging Spectroradiometer) in which the dispersed ash cloud is observed in the morning of 18 May 2016, and the Mt. Etna area is shown in the cartographic data derived from satellite images, from Google Maps, Figure 2b.
The dual wavelength Lidar employed in this case study is a ground-based remote sensing system, permanently monitoring dispersed Etna ash plumes, with elevation and azimuth scanning abilities, as highlighted in the Figure 2c. This innovative Lidar system was developed in the framework of the VAMOS SEGURO project (http://www.vamosseguro.eu), with the aim of studying and forecasting volcanic ash plumes. The Lidar, shown in Figure 2d, named Aerosol Multi-Wavelength Polarization Lidar Experiment (AMPLE) and developed by the Consorzio Inter-Universitario per le Scienze Fisiche della Materia (CNISM), is a compact multi-wavelength elastic/Raman scanning system with cross-polarization capability. At present, the AMPLE system is part of the European Aerosol Research Lidar Network (EARLINET) network [4], and it is devoted to special measurement campaigns at Etna volcano. The Lidar system is operated

Dual Wavelength Polarimetric Lidar
The capability of Lidar systems to detect the smallest particles (from 0.1 to 100 µm) in volcanic plumes as well as to reliably estimate the ash concentration depends on the instrumental characteristics and style of eruptive activity. Lidar measurements, in terms of particle backscattering coefficient and cross-polarization ratio profiles, can be carried out in regions where the Lidar signal is not extinguished due to the volcanic plume optical thickness even though some path corrections can be applied [6]. Figure 2a shows an image obtained by MODIS spectrometers (Moderate Imaging Spectroradiometer) in which the dispersed ash cloud is observed in the morning of 18 May 2016, and the Mt. Etna area is shown in the cartographic data derived from satellite images, from Google Maps, Figure 2b.
The dual wavelength Lidar employed in this case study is a ground-based remote sensing system, permanently monitoring dispersed Etna ash plumes, with elevation and azimuth scanning abilities, as highlighted in the Figure 2c. This innovative Lidar system was developed in the framework of the VAMOS SEGURO project (http: //www.vamosseguro.eu), with the aim of studying and forecasting volcanic ash plumes. The Lidar, shown in Figure 2d, named Aerosol Multi-Wavelength Polarization Lidar Experiment (AMPLE) and developed by the Consorzio Inter-Universitario per le Scienze Fisiche della Materia (CNISM), is a compact multi-wavelength elastic/Raman scanning system with cross-polarization capability. At present, the AMPLE system is part of the European Aerosol Research Lidar Network (EARLINET) network [4], and it is devoted to special measurement campaigns at Etna volcano. The Lidar system is operated at the Istituto Nazionale di Astrofisica (INAF) in Catania about 25 km away from the Etna summit craters [1]. distribution, structure and chemical composition. In combination with other parameters, it can be used to distinguish between different aerosol typologies. This parameter is essential in the optical properties retrieval, and it is generally assumed to be known in the inversion procedure using the Klett-Fernald algorithm [1,6,18,19]. The Elastic/Raman technique, employed in the AMPLE system, allows for direct measurement of the LR in the volcanic plume, estimated to be 48 sr [1]. The detection system can measure both the polarimetric elastic Lidar returns at W0 and W1 (horizontal and vertical polarized signals) and the nitrogen (N2) Raman Lidar echoes at 386 nm. Each detected signal is acquired by a multi-channel scaler with a raw spatial resolution of 15 m [1].
The cross-polarization ratio is one of the primary parameters able to discriminate different aerosol components. The accuracy in the retrieval of cross-polarization ratio is the driving factor for assessing and improving this capability [16] and can be estimated using methods accepted by the scientific community [4,14].

Polarimetric Lidar Observations during the May 2016 Etna Eruption
Using the notation adopted in [6], Lidar elastic observations can be expressed in terms of the backscattering coefficient β xy (m −1 sr −1 ), hereinafter O1, i.e., first observable, where x = h, v (horizontal and vertical respectively) stands for the receiving mode and y = h, v for the transmitting mode polarization, as well as in terms of cross-polarization ratio coefficient δ cr = β vh /β hh (dimensionless), hereinafter O2, i.e., second observable. In this work, we focus on the processing of copolar horizontally polarized backscattering coefficient β hh and cross-polarization ratio β vh to derive δ cr .
The β hh coefficient from daytime measurements can be obtained by using the Klett-Fernald algorithm [1,17,18], known as the Lidar Ratio LR that is the ratio between the extinction coefficient and the backscatter coefficient. This parameter varies for aerosol types depending on the aerosol microphysical properties such as refractive index, size distribution, structure and chemical composition. In combination with other parameters, it can be used to distinguish between different aerosol typologies. This parameter is essential in the optical properties retrieval, and it is generally assumed to be known in the inversion Remote Sens. 2021, 13, 1728 5 of 19 procedure using the Klett-Fernald algorithm [1,6,18,19]. The Elastic/Raman technique, employed in the AMPLE system, allows for direct measurement of the LR in the volcanic plume, estimated to be 48 sr [1].
The δ cr value, generally referred to as δ a to show that it is an aerosol, is obtained from the elastic Lidar profiles measured in the horizontal and vertical-polarized channels at W0 and W1 according to inversion procedures in [20,21]. The cross-polarization ratio of the finest ash particles can be used to identify the presence of non-spherical particles with higher cross-polarization ratio in the atmospheric sample sounded by the Lidar [9]. Nevertheless, polarimetric measurements, whose uncertainty can affect the derived parameters accuracy, need a calibration procedure that takes into account the different efficiencies of the detection for the parallel and perpendicular polarization channels. In order to calibrate the depolarization channels of our system, we used a two-step procedure already described in [1].
Polarimetric Lidar observables in the W1 during the 2016 Etna case study are shown in the panels of Figure 3. In the upper panel, the altitude-time diagrams of the range-corrected Lidar backscattered signal (RCS). RCS is the Lidar signal multiplied by the square of the distance and is provided as arbitrary unit (a. u.). The cross-polarization ratio (DEP) allows to highlight the dispersed ash layer using four elevations: (1)  Note that the maximum distance observed from the Lidar along the beam direction is always 15,400 m. By varying the pointing angle, the same ash layer is seen at different distances. The different elevations were chosen in order to observe different portions of the ash cloud and to probe the various stratifications, taking into account the attenuation that the signal undergoes when it crosses the atmospheric particulate along the path between the Lidar and the volcanic cloud.
In Figure 3a,b, we also identify two ash layers as observed by Lidar at different elevation pointing angles and between 2-4 km and 6-7 km of altitude. The horizontal structure of the ash cloud layer is producing a signature whose range is indeed increasing as the elevation angle decreases. In fact, the two layers are fairly distinguishable for a zenith pointing. For lower pointing angles, we observed an increase in the optical thickness due to the greater amount of ash layer crossed. This fact does not allow distinguishing the highest ash layer between 6 and 7 km due to strong attenuation of the Lidar signal. Lidar returns before 2 km are mainly due to aerosols of anthropogenic origin and very fine particles of water that constitute the particulate of atmospheric boundary layer (within 2 km). This horizontal stratification generated a strong backscatter return before 14:15 UTC, probably associated with cirrocumulus clouds of ice. However, a thin layer between 6 and 7 km could be related to the volcanic clouds of the VOR, which reached 6.5 km as estimated by the satellite [15]. It is worth noting that for this event no dust cloud load was forecasted over Catania area (https://ess.bsc.es/bsc-dust-daily-forecast), confirming that the layering observed by Lidar is of volcanic origin.
In order to deepen the characterization of the detected ash layers, Figure 4 shows the vertical profiles of backscattering coefficient and cross-polarization ratio, analysed at 15:02 using an elevation angle of 90 • at W0 (red line) and W1 (blue line) [1]. The panels in the middle and in the right show the zoom of profiles including the signature between 2 and 4 km and between 6 and 7 km, respectively. Remote Sens. 2021, 13, x FOR PEER REVIEW 6 of 21 In order to deepen the characterization of the detected ash layers, Figure 4 shows the vertical profiles of backscattering coefficient and cross-polarization ratio, analysed at 15:02 using an elevation angle of 90° at W0 (red line) and W1 (blue line) [1]. The panels in the middle and in the right show the zoom of profiles including the signature between 2 and 4 km and between 6 and 7 km, respectively.
We note that the backscattering coefficient at W1 ranges between 0 and 2.5 m −1 sr −1 and at W0 between 0.5 and 3.3 m −1 sr −1 between 2 and 4 km (green dashed square). The cross-polarization ratio at W1 ranges between 0.05 and 0.22% and at W0 between 0.13 and 0.1%. Between 6 and 7 km (blue dashed square), the backscattering coefficient at W1 ranges between 0 and 0.6 m −1 sr −1 and at W0 between 0.8 and 1.6 m −1 sr −1 , whereas the crosspolarization ratio at W1 ranges between 0.1 and 0.5 % and at W0 between 0.08 and 1.6 %. We note that the backscattering coefficient at W1 ranges between 0 and 2.5 m −1 sr −1 and at W0 between 0.5 and 3.3 m −1 sr −1 between 2 and 4 km (green dashed square). The cross-polarization ratio at W1 ranges between 0.05 and 0.22% and at W0 between 0.13 and 0.1%. Between 6 and 7 km (blue dashed square), the backscattering coefficient at W1 ranges between 0 and 0.6 m −1 sr −1 and at W0 between 0.8 and 1.6 m −1 sr −1 , whereas the cross-polarization ratio at W1 ranges between 0.1 and 0.5 % and at W0 between 0.08 and 1.6%.
We can highlight how both cross-polarization ratios at W1 and W0 show their maximum value, indicating the presence of larger aspherical ash particles whose abundance is also manifested by larger backscattering values. We can highlight how both cross-polarization ratios at W1 and W0 show their maximum value, indicating the presence of larger aspherical ash particles whose abundance is also manifested by larger backscattering values. . The blue and green dashed squares identify the profile portions of ash layers. Lidar profiles at an altitude between 7 and 6 km, panels (c) and (d), and between 4 and 2 km using a zenith pointing, panels (e) and (f). Both cross-polarization ratios at W1 and W0 around the same altitude of 3.5 km show two maximum values associable with two ash layers. The purple hatched ellipses identify the peaks where the ash layers are present in the panels (e) and (f). . The blue and green dashed squares identify the profile portions of ash layers. Lidar profiles at an altitude between 7 and 6 km, panels (c,d), and between 4 and 2 km using a zenith pointing, panels (e,f). Both cross-polarization ratios at W1 and W0 around the same altitude of 3.5 km show two maximum values associable with two ash layers. The purple hatched ellipses identify the peaks where the ash layers are present in the panels (e,f).

Retrieval Methods from Lidar Data
A physically based inversion methodology, named Volcanic Ash Lidar Retrieval (VALR) and previously developed [6], is applied to estimate ash particle mean size and mass concentration in the investigated volcanic ash cloud. Several retrieval approaches, such as the VALR Maximum-Likelihood (ML), Single-Regressive (SR) and Multi-Regressive (MR) algorithms, are implemented in this analysis in order to explore potentialities and differences for retrieving ash profiles and compare their results. The goals are: (i) using the VALR ML forward model for dual-polarization and dual-wavelength Lidar observations and (ii) comparing VALR ML with SR and MR retrieval techniques and parametric relationships. For this analysis, we assume that the VALR ML is the reference one due to its capability to deal with non-linearity of the inversion processes.
In order to infer volcanic ash parameters from Lidar measurements, we need to solve an ill-posed inverse problem [22]. The latter can be approached by resorting to a physically based approach that translates into the solution of the related forward problem, i.e., the modelling of the Lidar multi-wavelength and polarimetric response from a microphysical description of the volcanic particle distribution within ash cloud layers [6]. Note that the pre-processing of the Lidar signal is also devoted to the elimination of possible sources of error (e.g., path attenuation, ground clutter, target blocking).
The physically based approach requests an ash particle microphysical model. Hence, we here adopt a scattering model based on the T-Matrix solution method for spheroidal particles [23]. The numerical implementation, previously developed and named Tephra Particle Ensemble Scattering Simulator (TPESS) numerical simulator [6], considers volcanic particles features (size, density, shape, refractivity), derived from in situ experimental data and the available literature. The rationale of TPESS forward modelling is that using a Monte Carlo approach with a uniform probability, we let all the microphysical and morphological driving parameters of ash dimensional spectra to vary within a constrained range of values. Assuming equivalent-spheroid shapes, we apply the T-matrix backscattering simulator to compare Lidar observations. For Lidar applications, we identify possible 180 classes as a combination of three different particle sizes (VA, FA, CA), five different orientations (TO.1, TO.2, TO.3, OO, PO), four ash mass concentrations (VC, SC, MC, IC) and three axis ratios r ax = l/w, length l over width w of particle [24] (RB, RR and SR) as shown in the Table 1 (see [6] for more details). Gaussian-PDF Gaussian-PDF Gaussian-PDF The most probable classes, that we can find above the Lidar station, are the VA and FA ones, having the smallest concentrations (VC and SC). Figure 5 shows that the classes with the best correspondence with the Lidar measures are the classes FA-SC (yellow points) and VA-VC (red points). Applying the VALR ML classifier, the most probable simulated size class is the VA, with an occurrence of 95% respect to the FA (about 5%). Given the low percentage of detected FA, we restrict the analysis only to the VA class characterized by Very small Concentration (VC) category with values ranging between 10 −6 and 10 −3 g/m 3 [6]. Figure 5 also shows the overlapping between the simulated (FA-SC in yellow dots, FA-VC in blue dots, VA-SC green dots and VA-VC in red dots) and AMPLE Lidar measurements measured (dark dots) backscattering coefficient β hhs and cross-polarization ratio δ crs (the subscript s stands for simulated) at W0 and W1. Note that we will prefer to express β hhs in dBβ, that is, a value in decibel equal to 10 log10(β hhs ) when β hhs is expressed in Remote Sens. 2021, 13, x FOR PEER REVIEW 10 of 21 10 −6 and 10 −3 g/m 3 [6]. Figure 5 also shows the overlapping between the simulated (FA-SC in yellow dots, FA-VC in blue dots, VA-SC green dots and VA-VC in red dots) and AMPLE Lidar measurements measured (dark dots) backscattering coefficient βhhs and cross-polarization ratio δcrs (the subscript s stands for simulated) at W0 and W1. Note that we will prefer to express βhhs in dBβ, that is, a value in decibel equal to 10 log10(βhhs) when βhhs is expressed in [m -1 sr -1 ] [6]. The synthetic and measured signatures of the Lidar are in good agreement, confirming that the selected VA-VC class well represents the ash particles observed by the Lidar between 2-7 km of altitude and 25 km distance from VOR. The distribution of the coloured points that are distributed along the axis of the backscatter coefficient of Figure  5 mainly identifies the presence of the smallest and spherical particles that show depolarization values close to zero. The VALR ML can treat a variety of input data: (i) a single observable such as the measured copolar backscattering coefficient βhhm (denoted hereinafter by O1); (ii) double observables such as measured copolar backscattering βhhm and cross-polarization ratio coefficients δcrm (denoted hereinafter by O2); (iii) a single or double observable at a single wavelength W0 or W1, respectively, and at both wavelengths W0 and W1 simultaneously (denoted hereinafter W2). Table 2    The synthetic and measured signatures of the Lidar are in good agreement, confirming that the selected VA-VC class well represents the ash particles observed by the Lidar between 2-7 km of altitude and 25 km distance from VOR. The distribution of the coloured points that are distributed along the axis of the backscatter coefficient of Figure 5 mainly identifies the presence of the smallest and spherical particles that show depolarization values close to zero.
The VALR ML can treat a variety of input data: (i) a single observable such as the measured copolar backscattering coefficient β hhm (denoted hereinafter by O1); (ii) double observables such as measured copolar backscattering β hhm and cross-polarization ratio coefficients δ crm (denoted hereinafter by O2); (iii) a single or double observable at a single wavelength W0 or W1, respectively, and at both wavelengths W0 and W1 simultaneously (denoted hereinafter W2). Table 2 summarizes the possible combinations of the retrieval techniques using single or double observables and wavelengths: (1) the Single Regressive (SR) technique, using power laws with one (O1) or both observables (O2) connected to the wavelength W0 or W1; (2) the Multiple Regressive (MR) technique, using power laws with one O1 observable simultaneously and both wavelengths W2; (3) the Maximum Likelihood (ML) estimation with different observable combination.

Maximum Likelihood (ML) Method
In the ML approach, the ash particle classification from a set of Lidar observables is performed using a Maximum a Posteriori (MAP) probability criterion.
The VALR ML approach basically reduces to a minimization process in order to assign the specific class to each available Lidar measurement. Under a Gaussian probability density function of the differences between simulated and measured observables, the VALR ML method minimizes a quadratic form. The estimated ash class c and the retrieved microphysical parameters are those that exhibit the minimum metric, computed as the average of first five minimum square distance d 2 between the Lidar measurement set (β hhm , δ crm ) and simulated set (β hhsWX , δ crsWX ) using the selected VA-VC class at the specific wavelengths WX, i.e., W0 or W1. We can write the quadratic distance for a single wavelength as Being each term of (1) weighted by the inverse of the variances σ 2 βhhs and σ 2 δcrs of the simulated ash class dataset, whereas C a and D n are the mass concentration (in mg/m 3 ) and mean diameter (in µm), respectively.
We can improve the VALR ML estimation, considering both observables measured at both wavelengths W0 and W1 at the same time, obtaining the following quadratic distance at W2: In order to retrieve the ash mass concentration and mean diameter in the selected VA-VC class, we can extract their value by minimizing the quadratic distance in (1) or (2), that is,Ĉ where Min is the minimum function. The symbol "|" in (3a) and (3b) indicates that the estimated concentration and mean diameter are those values minimizing the quadratic distance as a function of the possible parameters of the simulated class. It is worth noting that these retrievals are conditioned by microphysical-electromagnetic assumptions in the TPESS numerical model and their representativeness with respect to the observed scene.

Single and Multiple Regression (SR and MR) Methods
Starting from TPESS simulated data for the VA-VC class, we can derive the regressive power laws correlating ash concentration C a and mean diameter D n .
For the single-wavelength regression (SR), the Equations (4a) and (4b) use only the measured backscattering coefficient or both backscattering coefficient and cross-polarization ratio, respectively, for each wavelength that is considered separately: where the coefficients a, b, c, d, e, and f are the SR regression coefficients and WX can be W0 or W1. For the multiple-wavelength regression (MR), the following Equation (5) employs only the backscattering coefficient at both wavelength W0 and W1 together (MR) to estimate C a and D n :Ĉ where a , b , c , d , e , f are the MR regression coefficients. The choice to use only the backscattering coefficient in MR is due to the difficult use of the cross-polarization ratio with the regressive models.
To evaluate the discrepancy ε of C a and D n estimations, expressed in percentage and due to the different retrieval methods, we can compute the normalized absolute value of the difference between SR and VALR-ML as well as between MR and VALR-ML retrievals. We then calculate the average ε Ca% , ε Dn% along the vertical profile of all discrepancies within the detected ash layers (2, 4, 6 and 7 km). The dispersion of these values around the mean value ε is expressed in terms of standard deviation σ Ca% , σ Dn% . In formulas we have where Mean z and Std z are the mean and standard deviation function, respectively, along the vertical coordinate z. Note that in (6), in absence of in situ airborne data, we consider the VALR ML estimates as the 'reference approach' because it can deal with non-linearities of the inverse problems much better than statistical regression methods. In this respect, we expect that the VALR ML estimates are more realistic than SR and MR methods.

Parametric Inversion Methods at Visible Wavelength
We also apply two parametric models, available in literature, using only measured backscatter coefficient (O1) in the visible wavelength (W1) [6,7,11,13]. The first retrieval parametric model (hereinafter P1) is expressed bŷ where k c is the ash conversion factor, function of the PSD and mainly dependent on the ash effective radius r e k c = 2 3 r e [6,9,11,25] with r e equal to half value of D n from the VALR ML algorithm. In this case, we can use the point estimates of D n (hereinafter PML1) or the average of its values within the ash layer (hereinafter PML2). LR is the mean value of Lidar ratio (48) [7,26]; ρ a is the density of volcanic ash fixed to 2450 kg/m 3 [27], and β hhW1 is the measured volcanic ash backscatter coefficient at W1.
The second parametric method [6,9] used to derive C a (hereinafter P2) employs the measured backscattering coefficient at W1 and can be written aŝ The expression between square brackets is known as the mass-extinction conversion factor for volcanic ash concentration, depending on the particle effective radius in this case equal to half mean diameter previously retrieved. [6,9,28]. Both parametric P1 and P2 models need some a priori information.

Results
In this section, we summarize the ash mass concentration and/or mean diameter estimates obtained from the previously described retrieval algorithms, and we show the main differences among the various estimates. In the model-based approach, we have the ML, SR and MR techniques using only one or more combinations of backscattering and cross-polarization ratio (O1, O2) at different wavelengths (W0, W1 or W2; see Table 2). We finally apply four parametric approaches, named P1, P2, PML1 and PML2. Hence, we evaluate the discrepancies through the metrics defined in Equation (6). Figure 6 shows the different retrievals of C a and D n , respectively, observed at the zenith angle (90 • ) using all the retrieval approaches. Figure 6b,d,f,h,j,l shows a variability of the estimated values falling within the VA-VC simulated class boundaries. We compute the standard deviation of first five minimum square distance d 2 (Equation (2)) in order to show an uncertainty (±σ), closest to ML estimates, plotted as coloured area bounded by the blue dashed line. The peaks observable in the profiles of the mass concentration estimates (6.9, 6.7 and 6.3 km, between 3.8 and 3.2 km and between 2.8 and 2.2 km of altitude, panels (f) and (h), respectively) show a stratification of ash clouds in thin layers of ash.
The SR retrieval models (O1-purple line and O2-yellow line) show results with less variable values. The estimates are consistent with each other and with the ML estimates, using the various combinations. We consider the parametric models only at W1, Figure 6f,h. P1 model shows the same variability of the vertical profile typical of the ML method, whereas the other two models P2, PML1 and PML2, using a fixed r e , show a more uniform trend. In all cases, C a estimates vary between 1 and 400 µg/m 3 . The P1 estimates are greater than the ML and SR estimates, reaching values of 400 µg/m 3 , whereas the P2, PML1 and PML2 models do not exceed 200 µg/m 3 . At W0, the C a estimates, Figure 6b,d, are more consistent. The ML estimates are more variable, probably due to a stratification of the observed layer along the vertical profile, as highlighted by the ML method.
Combining the measurements at W0 and W1, that is W2, the ML-O2 and MR-O1 estimates are shown in Figure 6j,l. The ML estimates show a greater variability (with a minimum around 1 µg/m 3 and a maximum at 300 µg/m 3 ) than the MR estimates (with a minimum trend around 10 µg/m 3 and a maximum at 50 µg/m 3 ). For the mean diameter D n estimates, shown in the left side panels of Figure 6, we always consider the bands W0, W1 and W2. In W1, the ML estimates with two observables (O2), as a function of uncertainty σ, and the SR estimates with one or two observables (O1) and (O2) are shown in Figure 6e,g. For the ML method D n the values range between 1 µm and 6 µm, as seen in the blue area of each panel, whereas for the SR method between 1 µm and 2 µm.
Using W0, Figure 6a,c, the D n estimates show a narrower variation of size values. The retrieved D n values vary between 0.5 µm and 3.5 µm, whereas SR-O1 and SR-O2 between 0.2 µm and 1.3 µm. Figure 6i,k shows the ML-O2 and MR-O2 estimates with estimated D n range between 0.1 µm and 6 µm and oscillating between 0.5 µm and 1.5 µm, respectively. We note a reduction in the fidelity of ML, SR and MR estimates, both in the mass concentration and mean diameter, with increasing observation distance and with the degrading of the measured signal. Moreover, in this case, the D n estimates of the ML with their uncertainty allow to collect the regressive estimates in their dynamics.
The AMPLE Lidar can scan and observe along different azimuth and elevation, as shown in Figures 2 and 3. Figure 7 shows C a (panels (a), (c) and (d)) and D n (panels (b), (d) and (f)) retrievals, applying the various methodologies, previously discussed, to Lidar measurements performed only at W1, coinciding with the thin stratification of ash at low altitude (between 2 and 4 km with a zenith point). By increasing the pointing angle from 20 • to 30 • and 40 • , the C a and D n retrievals become more similar. This behaviour is strictly due to the increase of the atmospheric boundary layer crossed along the distance between the ash layer and Lidar and consequently of the optical thickness which attenuates the Lidar signal, not allowing to discriminate the ash stratification.  , panels (f) and (h) and combining both wavelengths (W2), panels (j) and (l), respectively, and related to layers between 2 and 4 km (panels (d), (h) and (l)) and between 6 and 7 km, (panels (b), f); ML uses both observables (O2, backscattering coefficient and cross-polarization ratio) and the coloured area bounded by the blue dashed line is obtained considering the discrepancy σ around the estimated value (line blue); SR uses only the backscattering coefficient (O1) and both Lidar observables (O2); MR uses only the backscattering coefficient (O1); PM1, PM2, PML1 and PML2 use only the backscattering coefficient. The six panels on the left show the retrievals performed at W1, panel (a) and (c), W0, panel (e) and (g) and combining both wavelengths W2, panel (i) and (k), respectively, and related to layers between 2 and 4 km, panels c), g) and (k)) and between 6 and 7 km, panels a), e) and (i). The employed retrieval methodologies are: ML using both observables (O2); SR-O1 and O2; MR-O1.
The AMPLE Lidar can scan and observe along different azimuth and elevation, as shown in Figures 2 and 3. Figure 7 shows (panels (a), (c) and (d)) and (panels (b), (d) and (f)) retrievals, applying the various methodologies, previously discussed, to Lidar measurements performed only at W1, coinciding with the thin stratification of ash at low altitude (between 2 and 4 km with a zenith point). By increasing the pointing angle from 20° to 30° and 40°, the and retrievals become more similar. This behaviour is strictly due to the increase of the atmospheric boundary layer crossed along the distance between the ash layer and Lidar and consequently of the optical thickness which attenuates the Lidar signal, not allowing to discriminate the ash stratification.
Regarding the ash mass concentration, at 20° elevation, the profiles show a distance between 6 and 10 km. Among the parametric models (P1, P2, PML1 and PML2, using backscattering coefficient O1), only the P1 method is confirmed as the one that reaches the greatest values (about 10 3 g/m 3 ), followed by PML1 and PML2. ML-O2, SR-O1 and O2 show a similar trend with a variability between 1 g/m 3 and peaks reaching 400 g/m 3 . At an altitude of 30°, the estimates between 4 and 8 km tend to be comparable between different methods, with the exception of the P1 method which also in this case is confirmed as the one with the highest values (about 10 3 g/m 3 ). The ML method ranges Figure 6. Vertical profile of mean diameter D n (panels (a,c,e,g,i,k)) and ash mass concentration C a (panels (b,d,f,h,j,l)) retrievals observed at 15:02 UTC on 18 May 2016. The panels on the right show C a retrievals performed at W1, panels (b,d); at W0, panels (f,h) and combining both wavelengths (W2), panels (j,l), respectively, and related to layers between 2 and 4 km (panels (d,h,l)) and between 6 and 7 km, (panels (b,f)); ML uses both observables (O2, backscattering coefficient and cross-polarization ratio) and the coloured area bounded by the blue dashed line is obtained considering the discrepancy ±σ around the estimated value (line blue); SR uses only the backscattering coefficient (O1) and both Lidar observables (O2); MR uses only the backscattering coefficient (O1); PM1, PM2, PML1 and PML2 use only the backscattering coefficient. The six panels on the left show the D n retrievals performed at W1, panel (a,c), W0, panel (e,g) and combining both wavelengths W2, panel (i,k), respectively, and related to layers between 2 and 4 km, panels (c,g,k)) and between 6 and 7 km, panels (a,e,i). The employed retrieval methodologies are: ML using both observables (O2); SR-O1 and O2; MR-O1.
Regarding the ash mass concentration, at 20 • elevation, the profiles show a distance between 6 and 10 km. Among the parametric models (P1, P2, PML1 and PML2, using backscattering coefficient O1), only the P1 method is confirmed as the one that reaches the greatest values (about 10 3 µg/m 3 ), followed by PML1 and PML2. ML-O2, SR-O1 and O2 show a similar trend with a variability between 1 µg/m 3 and peaks reaching 400 µg/m 3 . At an altitude of 30 • , the estimates between 4 and 8 km tend to be comparable between different methods, with the exception of the P1 method which also in this case is confirmed as the one with the highest values (about 10 3 µg/m 3 ). The ML method ranges include the other estimates (SR and P2, PML1 and PML2), with values ranging from 1 µg/m 3 to 300 µg/m 3 . At 40 • of elevation, we find the same dynamics observed in the other two elevations for a distance between 3 and 6 km, that is, the overlap between the various regressive and parametric methods with the ML one, probably related to the best observation geometry of the AMPLE Lidar. observation geometry of the AMPLE Lidar.
The detection approach discussed for the estimates can be basically proposed for the estimates. For a Lidar elevation angle of 20°, 30° and 40°, the ML method has significant oscillations. In contrast, the estimates with SR-O1 and O2 are quite similar. In both cases, the estimates show a variability between 10 −1 and 6 m. The retrievals at 40° are more consistent ranging between 10 −1 and 5 m. It should be noted that for the different elevation angles (20°, 30° and 40°), we took into account only the distance relative to the two lowest ash layers, because the lidar signal is not degraded. In Table 3 are listed the discrepancies and the standard deviations for and estimates in percentage, using indices defined in Equation (6). The addition of the second observable (depolarization or cross-polarization ratio, the latter used hereinafter) in the regressive SR and MR estimates increases the discrepancy for both and estimates. Regarding the estimates at zenith pointing, for W1 case, this effect is more evident, The detection approach discussed for the C a estimates can be basically proposed for the D n estimates. For a Lidar elevation angle of 20 • , 30 • and 40 • , the ML method has significant oscillations. In contrast, the D n estimates with SR-O1 and O2 are quite similar. In both cases, the estimates show a variability between 10 −1 and 6 µm. The D n retrievals at 40 • are more consistent ranging between 10 −1 and 5 µm. It should be noted that for the different elevation angles (20 • , 30 • and 40 • ), we took into account only the distance relative to the two lowest ash layers, because the lidar signal is not degraded.
In Table 3 are listed the discrepancies and the standard deviations for C a and D n estimates in percentage, using indices defined in Equation (6). The addition of the second observable (depolarization or cross-polarization ratio, the latter used hereinafter) in the regressive SR and MR estimates increases the discrepancy for both C a and D n estimates. Regarding the estimates at zenith pointing, for W1 case, this effect is more evident, whereas for W0 case, there is an opposite behaviour, i.e., a reduction in discrepancy when two observables are considered. Although the SR method shows low discrepancy, combining W0 and W1 and using a single observable (O1) in MR shows a greater discrepancy for both C a and D n estimates respect to SR methods. This trend is quite intuitive considering the different sensitivity of W0 and W1 radiation to particle size. Table 3. Overview of the mean discrepancies ε Ca% and ε Dn% and the relative dispersions σ Ca and σ Dn , expressed in percentage and related, respectively, to C a and D n , using the 4 different angles (20 • , 30 • , 40 • and 90 • ) performed during Lidar measurements and considering the possible combinations of observables (O1 and O2) and wavelengths (W0, W1 and W2) for the ML, SR, MR and parametric retrieval methods. Green colours indicate low discrepancy (ε Ca% , ε Dn% ≤ 10), blue moderate discrepancy (10 < ε Ca% , ε Dn% ≤ 30), yellow large discrepancy (30 < ε Ca% , ε Dn% ≤ 60), and finally, red extreme discrepancy (ε Ca% , ε Dn% > 60). Considering the parametric models, the estimates obtained using only the backscattering coefficient measured with a pointing of 90 • show variable discrepancy values. The greater C a values derived from the P1 are probably due to the value of r e associated a priori, compared to the other parametric models that use r e values deduced from the D n estimates by the ML method. If we look at the estimates and errors associated with the measurements in the W1 alone, if the pointing angle increases from 20 • to 40 • , the discrepancy terms relating to C a and D n estimates tend to decrease. This occurs, in our opinion, because the quality of the Lidar data can increase when the observed distance reduces affecting the regressive and parametric estimates. This confirms that, for single wavelength observations, the observation angle can influence this inter-comparison, being the elevation angle of 40 • the optimal case for a minimum discrepancy. In fact, at low elevation, the distance at which the observed ash layer is intercepted increases, increasing the mass of air crossed and consequently the measurement errors.

Altitude (km)/Elevation ( • )/
To better highlight the potential of the various methods, the boxes in Table 3 are shown with different colours according to the different discrepancy values. All the regressive methods (SR and MR) are those whose estimates match best with ML (ε Ca% , ε Dn% ≤ 10), as long as the aiming angle is greater than 30 • . The parametric methods are strongly affected by the a priori assumptions leading the C a estimates to differ more from those ML. The P1 method is confirmed as the one that differs most from the ML method (ε Ca%~1 00).

Conclusions
In this work, we have highlighted how the VALR ML algorithm, used in inversion of lidar data to derive particle size and concentration, can offer robust results. This result can be derived from the comparison of the results of the various regressive and parametric estimation methods. The analysis of the 2016 explosive event at Etna volcano highlights the potential of the Lidar dual-wavelength methodology compared to a single wavelength one to better describe both ash concentration and size. We have shown that a more detailed characterization of the dispersed ash layers at different altitude and distance from the volcanic vent is possible in terms of ash mass concentration C a and mean size D n retrievals. Reliable ash concentration and size retrievals from Lidars, such as AMPLE, can significantly help improve the capability to monitor and predict the airborne distribution of the volcanic ash clouds during transport and dispersion processes.
The VALR ML algorithm, applied to Lidar polarimetric observables, can recover the vertical profile of ash concentration and particle size in a physically coherent framework. The scanning capability of the dual wavelength polarimetric Lidar improves the detection of the concentration of finer ash of scattered plumes by combining the measurements at UV with those at VIS wavelengths. Given the flexibility of the VALR ML, it is possible to manage the non-linearity of the inverse problem, although the estimates may show more variability than regressive statistical approaches, such as SR and MR. On the contrary, the regressive approach can offer the advantage of more uniform estimates than ML ones, being less sensitive to variation, and of less intensive computational resources to execute operationally.
In this 2016 Etna case study, the ash concentration normally observed within thin ash layers ranges between 0.1 µg/m 3 and about 1 mg/m 3 , whereas the very fine particle sizes characterized by a mean diameter ranging between 0.1-6 µm is in agreement with particle size retrieved in [11]. This different concentration is expected considering the stratification at different levels of the finest dispersed ash. The different estimates of ash mean size and concentration at UV and VIS is mainly related to their different spectral sensitivity to the smaller particles present in the ash layers between 2 and 4 km and 6 and 7 km. The advantage of using both wavelengths in the estimation of C a and D n is to compensate for the retrievals between the estimates at the two bands with C a values between 0.1 µg/m 3 and about 1 mg/m 3 and of D n between 0.1 and 6 µm. A comparison between the concentration C a and average diameter D n estimates has been evaluated in terms of the discrepancy between the ML method and other regressive and parametric methods, obtaining values below 10% in most cases, in other cases values between 20% and 30% and only in the parametric approach P1 the discrepancy is about 100%.
Further work should be devoted to the verification of the ash cloud volcanic particle distribution, estimated from the proposed VALR approaches, possibly using airbornebased particle samplers flying within the ash clouds. In this respect, remotely controlled drones and unmanned vehicles can offer a unique opportunity for Lidar product validation campaigns.