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

Abstract

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/m³ and 1 mg/m³ 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.

1. 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]).

Light detection and ranging (Lidar) techniques represent the optimal tool to provide the range-resolved aerosol data [4]. Lidar systems are powerful techniques for monitoring dispersed particles in the troposphere and lower stratosphere because of their profiling capability at very high range resolution.

In order to mitigate from the impact associated with volcanic ash, Lidar observations allow to perform immediate and accurate detection of volcanic plumes, quantify volcanic ash concentration and characterize optical properties of volcanic particles [1], improving modelling of dispersed volcanic ash clouds [5,6]. Lidar systems have been widely used to monitor the vertical profile of the volcanic aerosol and study aerosol associated with volcanic clouds [7]. Lidar observations can also provide cloud morphological properties (such as top and bottom levels as well as thickness) and microphysical properties such as mass concentration and effective diameters [8]. Using the depolarization channel, ash particle shape can be in principle discriminated [6,9,10]. The capability of Lidar systems to detect the finest particles in volcanic plume and reliably estimate the ash concentration mainly depends on instrumental characteristics and the type of volcanic explosive activity.

This work aims at exploring the potential of the Lidar dual-wavelength methodology compared to a single wavelength methodology to better describe both ash concentration and size using the 18 May 2016 explosive event at Etna Volcano (Italy). Etna is one of the most active volcanoes in the world, and for this reason, it is regularly monitored mainly by several ground-based instruments of the Istituto Nazionale di Geofisica e Vulcanologia, Osservatorio Etneo (INGV-OE) [11]. The Lidar instrument, named Aerosol Multi-Wavelength Polarization Lidar Experiment (AMPLE), was upgraded in 2016 to perform dual-polarization and dual-wavelength measurements of the Etna explosive activity, extending the capability of the already existing Lidar [6,12,13]. The main limitation in using Lidar techniques is the decrease of the signal due to the presence of optically thick cloud layers. However, this limitation is mainly related to large explosive volcanic eruptions [14].

In this work, we analyse the medium-intensity explosive activity produced in the afternoon of 18 May 2016 at Mt. Etna, in Italy, from the Etna summit craters. Lidar measurements are performed at the INAF- Catania Astrophysical Observatory (37°31′43.8″ N, 15°4′18.5″ E, 196 m above sea level (asl) about 25 km from the Etna summit pointing at two thin volcanic layers clearly visible and dispersed at the altitude of 2–4 km (two layers) and 6–7 km (one layer) above sea level. Both ash particle backscattering and particle linear depolarization profiles (called also cross-polarization ratio [6]) at 355 nm (hereinafter W0) and 532 nm (hereinafter W1) wavelengths, respectively, are obtained using different pointing elevation angles at 20, 30, 40 and 90 degrees. For the first time, a combined dual-polarized dual-wavelength Lidar measurement at Etna volcano is here investigated using VALR ML. Uncertainties in the volcanic ash concentration and mean size estimations using different methods are also discussed.

The paper is organized as follows: Section 2 is devoted to describing the Etna case-study event as well as the dual wavelength polarimetric Lidar instrument and related measurements; Section 3 presents an overview of the retrieval methodologies developed to analyse the dispersed ash layers; Section 4 discusses the results, whereas the conclusions are drawn in Section 5.

2. Etna Case Study and Polarimetric Lidar Observations

2.1. Etna Case Study

Etna volcano generated an intense and complex eruptive activity between 16 and 25 May 2016, affecting several summit craters: North-Eastern Crater (NEC), Voragine (VOR), Bocca Nuova (BN), and Southeast Crater (SEC). In the night between 15 and 16 May, Strombolian activity was observed at the SEC and the day after intra-crater activity from the NEC was also visible by the INGV-OE monitoring network. In the morning of 18 May, ash emission from the NEC was mainly dispersed toward the south-east direction. At about 10:50 UTC, Strombolian activity from the VOR was clearly visible by the video-surveillance system. This activity anticipated by about less than 15 minutes the lava fountain that formed an ash plume directed toward ESE and having an altitude between of about 6.5 km above sea level [15]. The lava fountain activity ended after about 14:00 GMT and was replaced again by the Strombolian activity. This eruption was also characterized by the formation of lava flows, one of them formed on 18 May, and by two other lava flows on 19 May. The Etna explosive activity on May 18, 2016, at 11:46 UTC, observed by the EMOV and EMOT cameras of INGV-OE is shown in Figure 1a,b, respectively.

2.2. 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].

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].

2.3. 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⁻¹ sr⁻¹), 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 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) 90° in the time range between 13:45 and 15:16 and between 16:17 and 16:35 UTC, with the altitude limit of 15,400 m (all the heights are above sea level for 90°); (2) 20° between 15:20 and 15:55 UTC, with the distance limit of 15,400 m; (3) 30° between 15:55 and 16:07 UTC, with distance limit of 15,400 m and (4) 40° between 16:07 and 16:17 UTC, with distance limit of 15,400 m. 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.

We note that the backscattering coefficient at W1 ranges between 0 and 2.5 m⁻¹ sr⁻¹ and at W0 between 0.5 and 3.3 m⁻¹ sr⁻¹ 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⁻¹ sr⁻¹ and at W0 between 0.8 and 1.6 m⁻¹ sr⁻¹, 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.

3. 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. The microphysical particle features used in the simulations have been marked with an asterisk.

Ash Particle Ensemble Property	Very Fine Ash (VA)	Fine Ash (FA)	Coarse Ash (CA)
Ash diameter variability range DDn (mm)	Uniform PDF DDn 0.125–8	Uniform PDF DDn 8–64	Uniform PDF DDn 64–512
Ash particle concentration variability range (mg/m3) VC: Very small Conc. SC: Small Conc. MC: Medium Conc. IC: Intense Conc.	Uniform PDF VC: 10−3–100 (*) SC: 100–102 (*) MC: 102–103 IC: 103–104	Uniform PDF VC: 10−3–100 (*) SC: 100–102 (*) MC: 102–103 IC: 103–104	Uniform PDF VC: 10−3–100 SC: 100–102 MC: 102–103 IC: 103–104
Ash size distribution shape parameter mp (dimensionless)	Scaled Gamma PSD U-PDF mp = 1–2	Scaled Gamma PSD U-PDF mp = 1–2	Scaled Gamma PSD U-PDF mp = 1–2
Ash particle density rp (g/cm3)	Uniform PDF rp = 0.5–2.5	Uniform PDF rp = 0.5–2.5	Uniform PDF rp = 0.5–2.5
Ash particle canting angle mean and deviation mf (°) and sf (°) TO.1: Tumbling Orientation TO.2: Tumbling Orientation TO.3: Tumbling Orientation OO: Oblate Orientation PO: Prolate Orientation	Gaussian-PDF mf = 30°; sf = 30° mf = 45°; sf = 30° (*) mf = 60°; sf = 30° mf = 0°; sf = 10° (*) mf = 90°; sf = 10° (*)	Gaussian-PDF mf = 30°; sf = 30° mf = 45°; sf = 30° (*) mf = 60°; sf = 30° mf = 0°; sf = 10° (*) mf = 90°; sf = 10° (*)	Gaussian-PDF mf = 30°; sf = 30° mf = 45°; sf = 30° mf = 60°; sf = 30° mf = 0°; sf = 10° mf = 90°; sf = 10°
Non-spherical particle axial ratio rax: axis ratio (dimensionless) RB: basaltic ratio RR: rhyolitic ratio SR: spherical ratio	rax = AR RB: rax-b (*) RR: rax-r SR: rax-s (*)	rax = AR RB: rax-b (*) RR: rax-r SR: rax-s (*)	rax = AR RB: rax-b RR: rax-r SR: rax-s

(see [6] for more details).

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⁻⁶ and 10⁻³ g/m³ [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⁻¹sr⁻¹] [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. Overview of the different retrieval methodologies as function of observables and Lidar wavelength combinations employed in this work.

Method	Observable (O)		Wavelength (W)
	βhh	βhh&δcr	UV	VIS	UV&VIS
SR	O1	O2	W0	W1
MR	O1				W2
VALR ML	O1	O2	W0	W1	W2

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.

3.1. 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² between the Lidar measurement set $\left({\beta }_{hhm}, {\delta }_{crm}\right)$ and simulated set $\left({\beta }_{hhsWX},{\delta }_{crsWX}\right)$ 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

$$d^2\left(C_a^{\left(c\right)},D_n^{\left(c\right)}\right)=\frac{{\left[{\beta }_{hhmWX}-{\beta }_{hhsWX}\left(C_a^{\left(c\right)},D_n^{\left(c\right)}\right)\right]}^2}{{\sigma }_{\beta hhsWX}^2}+\frac{{\left[{\delta }_{crmWX}-{\delta }_{crsWX}\left(C_a^{\left(c\right)},D_n^{\left(c\right)}\right)\right]}^2}{{\sigma }_{\delta crsWX}^2}$$

(1)

Being each term of (1) weighted by the inverse of the variances σ²_βhhs and σ²_δcrs of the simulated ash class dataset, whereas $C_a$ and $D_n$ are the mass concentration (in mg/m³) 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:

$$\begin{gathered} d^2\left(C_a^{\left(c\right)},D_n^{\left(c\right)}\right)=\frac{{\left[{\beta }_{hhmW0}-{\beta }_{hhsW0}\left(C_a^{\left(c\right)},D_n^{\left(c\right)}\right)\right]}^2}{{\sigma }_{\beta hhsW0}^2}+\frac{{\left[{\delta }_{crmW0}-{\delta }_{crsW0}\left(C_a^{\left(c\right)},D_n^{\left(c\right)}\right)\right]}^2}{{\sigma }_{\delta crsW0}^2}+ \\ +\frac{{\left[{\beta }_{hhmW1}-{\beta }_{hhsW1}\left(C_a^{\left(c\right)},D_n^{\left(c\right)}\right)\right]}^2}{{\sigma }_{\beta hhsW1}^2}+\frac{{\left[{\delta }_{crmW1}-{\delta }_{crsW1}\left(C_a^{\left(c\right)},D_n^{\left(c\right)}\right)\right]}^2}{{\sigma }_{\delta crsW1}^2} \end{gathered}$$

(2)

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,

$${\widehat{C}}_a^{\left(C\right)}=C_a^{\left(C\right)}|Mi{n}_{\left(C_a^{\left(c\right)},D_n^{\left(c\right)}\right)}\left\{d^2\left(C_a^{\left(c\right)},D_n^{\left(c\right)}\right)\right\}$$

(3a)

$${\widehat{D}}_n^{\left(C\right)}=D_n^{\left(C\right)}|Mi{n}_{\left(C_a^{\left(c\right)},D_n^{\left(c\right)}\right)}\left\{d^2\left(C_a^{\left(c\right)},D_n^{\left(c\right)}\right)\right\}$$

(3b)

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.

3.2. 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:

$${\widehat{C}}_a^{\left(c\right)}=a{\beta }_{hhWX}^b{\widehat{D}}_n^{\left(c\right)}=d{\beta }_{hhWX}^e$$

(4a)

$${\widehat{C}}_a^{\left(c\right)}=a{\beta }_{hhWX}^b{\delta }_{crWX}^c{\widehat{D}}_n^{\left(c\right)}=d{\beta }_{hhWX}^e{\delta }_{crWX}^f$$

(4b)

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$:

$${\widehat{C}}_a^{\left(c\right)}=a^{\prime }{\beta }_{hhW0}^{b^{\prime }}{\beta }_{hhW1}^{c^{\prime }}{\widehat{D}}_n^{\left(c\right)}=d^{\prime }{\beta }_{hhW0}^{e^{\prime }}{\beta }_{{hhW1}^{\prime }}^{f^{\prime }}$$

(5)

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 ${\overline{\epsilon }}_{Ca\%},{\overline{\epsilon }}_{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 $\overline{\epsilon }$ is expressed in terms of standard deviation ${\sigma }_{Ca\%},{\sigma }_{Dn\%}$. In formulas we have

$${\epsilon }_{Ca}=100\frac{\left|C_{aXR}-C_{aML}\right|}{C_{aML}},{\overline{\epsilon }}_{Ca\%}=Mea{n}_z\left({\epsilon }_{Ca\%}\right),{\sigma }_{Ca\%}=St{d}_z\left({\epsilon }_{Ca\%}\right)$$

(6a)

$${\epsilon }_{Dn}=100\frac{\left|D_{nXR}-D_{nML}\right|}{D_{nML}},{\overline{\epsilon }}_{Dn\%}=Mea{n}_z\left({\epsilon }_{Ca\%}\right),{\sigma }_{Dn\%}=St{d}_z\left({\epsilon }_{Dn\%}\right)$$

(6b)

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.

3.3. 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 by

$${\widehat{C}}_{aP1}=k_c\overline{LR}{\rho }_a{\beta }_{hhW1}$$

(7)

where k_c is the ash conversion factor, function of the PSD and mainly dependent on the ash effective radius r_e $k_c=\frac{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). $\overline{LR}$ is the mean value of Lidar ratio (48) [7,26]; ρ_a is the density of volcanic ash fixed to 2450 kg/m³ [27], and ${\beta }_{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 as

$${\widehat{C}}_{aP2}=\left[1.346r_e-0.156\right]\overline{LR}{\beta }_{hhW1}$$

(8)

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.

4. 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² (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³. The P1 estimates are greater than the ML and SR estimates, reaching values of 400 μg/m³, whereas the P2, PML1 and PML2 models do not exceed 200 μg/m³. 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³ and a maximum at 300 μg/m³) than the MR estimates (with a minimum trend around 10 μg/m³ and a maximum at 50 μg/m³). 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 Figure 2 and Figure 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.

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³ μg/m³), followed by PML1 and PML2. ML-O2, SR-O1 and O2 show a similar trend with a variability between 1 μg/m³ and peaks reaching 400 μg/m³. 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³ μg/m³). The ML method ranges include the other estimates (SR and P2, PML1 and PML2), with values ranging from 1 μg/m³ to 300 μg/m³. 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.

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⁻¹ and 6 μm. The $D_n$ retrievals at 40° are more consistent ranging between 10⁻¹ 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. Overview of the mean discrepancies ${\overline{\epsilon }}_{Ca\%}$ and ${\overline{\epsilon }}_{Dn\%}$ and the relative dispersions ${\sigma }_{Ca}$ and ${\sigma }_{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 (${\overline{\epsilon }}_{Ca\%}$, ${\overline{\epsilon }}_{Dn\%}$ ≤ 10), blue moderate discrepancy (10 < ${\overline{\epsilon }}_{Ca\%}$, ${\overline{\epsilon }}_{Dn\%}$ ≤ 30), yellow large discrepancy (30 < ${\overline{\epsilon }}_{Ca\%}$, ${\overline{\epsilon }}_{Dn\%}$ ≤ 60), and finally, red extreme discrepancy (${\overline{\epsilon }}_{Ca\%}$, ${\overline{\epsilon }}_{Dn\%}$ > 60).

Altitude (km)/Elevation (°)/		Discrepancy (%)	SR O1 W0	SR O2 W0	SR O1 W1	SR 02 W1	MR O1 W2	P1 O1 W1	P2 O1 W1	PML1 O1 W1	PML2 O1 W1
2–4 km	90°	${\overline{\epsilon }}_{Ca\%}$	4.68	2.74	1.74	3.14	8.86	100	11.32	31.40	32.89
		${\sigma }_{Ca\%}$	25.90	17.30	8.13	14.20	40.74	100	58.62	100	76.51
		${\overline{\epsilon }}_{Dn\%}$	0.38	0.06	0.07	0.38	0.01
		${\sigma }_{Dn\%}$	1.18	0.89	1.10	1.46	1.24
6–7 km	90°	${\overline{\epsilon }}_{Ca\%}$	1.15	3.21	2.70	2.63	2.41	100	0.46	5.34	7.73
		${\sigma }_{Ca\%}$	2.17	20.99	20.21	18.40	15.14	100	12.68	38.29	37.51
		${\overline{\epsilon }}_{Dn\%}$	0.78	0.86	1.45	2.65	0.96
		${\sigma }_{Dn\%}$	1.27	1.48	1.11	1.69	1.36
2–4 km	40°	${\overline{\epsilon }}_{Ca\%}$			3.18	5.12		100	4.94	18.00	22.55
		${\sigma }_{Ca\%}$			16.78	21.76		100	36.90	86.38	59.82
		${\overline{\epsilon }}_{Dn\%}$			0.02	0.18
		${\sigma }_{Dn\%}$			0.88	1.19
2–4 km	30°	${\overline{\epsilon }}_{Ca\%}$			25.89	5.66		100	0.72	26.93	9.91
		${\sigma }_{Ca\%}$			54.00	5.52		100	2.09	77.83	14.86
		${\overline{\epsilon }}_{Dn\%}$			1.10	3.90
		${\sigma }_{Dn\%}$			1.15	12.30
2–4 km	20°	${\overline{\epsilon }}_{Ca\%}$			7.40	9.02		100	7.38	26.14	43.08
		${\sigma }_{Ca\%}$			30.01	31.25		100	33.34	80.89	100
		${\overline{\epsilon }}_{Dn\%}$			0.62	2.43
		${\sigma }_{Dn\%}$			1.18	10.32

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.

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.

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 (${\overline{\epsilon }}_{Ca\%}$, ${\overline{\epsilon }}_{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 (${\overline{\epsilon }}_{Ca\%}$ ~100).

5. 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³ and about 1 mg/m³, 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³ and about 1 mg/m³ 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 airborne-based 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.