Volcanic Cloud Detection and Retrieval Using Satellite Multisensor Observations

: Satellite microwave (MW) and millimetre-wave (MMW) passive sensors can be used to detect volcanic clouds because of their sensitivity to larger volcanic particles (i.e., size bigger than 20 µm). In this work, we combine the MW-MMW observations with thermal-infrared (TIR) radio-metric data from the Low Earth Orbit (LEO) spectroradiometer to have a complete characterisation of volcanic plumes. We describe new physical-statistical methods, which combine machine learning techniques, aimed at detecting and retrieving volcanic clouds of two highly explosive eruptions: the 2014 Kelud and 2015 Calbuco test cases. For the detection procedure, we compare the well-known split-window methods with a machine learning algorithm named random forest (RF). Our work highlights how the machine learning method is suitable to detect volcanic clouds using different spectral signatures without fixing a threshold. Moreover, the RF model allows images to be automatically processed with promising results (90% of the area correctly identified). For the retrieval procedure of the mass of volcanic particles, we consider two methods, one based on the maximum likelihood estimation (MLE) and one using the neural network (NN) architecture. Results show a good comparison of the mass obtained using the MLE and NN methods for all the analysed bands. Summing the MW-MMW and TIR estimates, we obtain the following masses: 1.11 ± 0.40 10 11 kg (MLE method) and 1.32 ± 0.47 10 11 kg (NN method) for Kelud; 4.48 ± 1.61 10 10 kg (MLE method) and 4.32 ± 1.56 10 10 kg (NN method) for Calbuco. This work shows how machine learning techniques can be an effective tool for volcanic cloud detection and how the synergic use of the TIR and MW-MMW observations can give more accurate estimates of the near-source volcanic clouds.

time-space dispersion of volcanic particles can ensure an improvement in defining the golden standard for airline flights and prevent damages [1]. Orbiting satellite observations can provide a large amount of daily data. The global perspective offered by Geosynchronous Earth Orbit (GEO) and the Low Earth Orbit (LEO) satellite systems is of vital importance for monitoring volcanoes [4,5], especially those in remote and inaccessible areas. Data from the LEO satellite microwave and millimetre-wave (MW-MMW) radiometers, such as Advanced Technology Microwave Sounder (ATMS) and Microwave Humidity Sounder (MHS), can be combined with thermal-infrared (TIR) spectroradiometers, such as Visible Infrared Imaging Radiometer Suite (VIIRS) and Advanced Very High Resolution Radiometer (AVHRR), to have more accurate estimates of the near-source volcanic cloud [1]. Whereas the LEO infrared data analysis represents the classic approach in the study of explosive activity by satellite, given their remarkable spatial resolution and sensitivity to ash clouds, their brightness temperature (BT) difference signatures can saturate because of large amounts of tephra mass mainly near the volcanic source. Moreover, the presence of tephra with a radius larger than 15 µm cannot be detected using these sensors, due to their detection limits [1,6], leading to inaccurate detections of some properties in the proximal or near-source volcanic cloud. In support of the TIR sensors, the MW-MMW sensors allow larger volcanic particles (⪆ 50 µm) to be observed. Actually, the near-source volcanic plume does not typically extinguish the MW-MMW signal, especially in the first hours of the eruptive event [1], probing the entire eruption column [7]. However, one of the drawbacks of MW-MMW sensors is the poor spatial resolution (up to 32 km). In the past, the MW-MMW bands have also been exploited to study volcanic plumes and their electrification processes [8][9][10].
The goal of this paper is to show how the analysis of both MW-MMW and TIR spectral signatures improves the characterisation of the volcanic plume in terms of detection and mass retrieval. For the detection, we worked in a statistical learning framework to develop a model capable of automatically processing images without the arbitrary threshold choice. In terms of retrieval, the new developed radiative transfer model algorithm ( ) is used to estimate the total columnar content (TCC) and, in turn, the mass for both MW-MMW and TIR. In this respect, two minimisation techniques, the MLE [11][12][13] and the NN [14][15][16], are also compared and discussed. The paper is structured as follows: Section 2 introduces the satellite sensors, the detection methods, the and the retrieval methods; Section 3 illustrates the 2014 Kelud and the 2015 Calbuco eruptions with an explorative data analysis of the satellite products; Section 4 shows the detection and retrieval results in two subsections, one per each eruption; Section 5 discusses the results and addresses future developments.

Methods
This section is divided into four subsections; the first describes the considered sensors for the analysis, the second discusses the split-window and the RF detection methods, the third describes the and the fourth focuses on the empirical approach.

Satellite Data
For the purpose of this work, we observed volcanic plumes using data coming from the ATMS and VIIRS, both on board the Suomi-National Polar-orbiting Partnership (S-NPP) LEO weather satellite provided by National Oceanic and Atmospheric Administration (NOAA) (all the acronyms are in Appendix A). The ATMS sensor is a passive crosstrack total power MW radiometer that measures microwave energy emitted and scattered by the atmosphere. It provides daily global atmospheric temperature, moisture and pressure profiles. It collects information in the frequency window from 23.80 to 183. 31 GHz. This window is divided into 22 channels, with channels 1-2-16 detecting quasi-vertical polarisation, while the other 19 quasi-horizontal polarisation [17][18][19]. The spatial resolution is 16 km for the channels from 165.50 to 183.31 GHz, 32 km for the channels from 50.30 to 88.20 GHz and 75 km for the channels from 23.80 and 31.40 GHz. The VIIRS instrument is a whiskbroom scanning radiometer with an angular field of view of 112.56° in the cross-track direction and a swath width of 3000 km. It acquires information in 22 spectral bands from 0.412 to 12.01 µm. The measured spectrum range is divided into 16 moderate resolution bands with 750 m spatial resolution, 5 imaging resolution bands with 375 m spatial resolution, and 1 panchromatic day/night band with a spatial resolution of 375 m [20,21]. The level 1B products used in this work were downloaded from https://www.avl.class.noaa.gov/saa/products/welcome (accessed on 4 February 2023) for ATMS and VIIRS sensors and from https://eoportal.eumetsat.int/userMgmt/confirmed.faces (accessed on 4 February 2023) for MHS and AVHRR sensors [22,23]. The MHS and the AVHRR sensors were used only to train the RF model. The scripts used for the analyses were coded in Matlab and Python. Finally, the colour map plots were generated using the "colorcet" library, which is a collection of perceptually uniform colour maps (https://colorcet.com/index.html, accessed on 4 February 2023).

Volcanic Cloud Detection
In this subsection, we focus on the detection of volcanic clouds. We first present the multi-spectral brightness temperature difference (or split-window) method (see Section 2.2.1), and then the RF method (see Section 2.2.2).

Multi-Spectral Brightness Temperature Difference
Volcanic clouds in the infrared spectrum can be identified by selecting pixels with a brightness temperature difference (BTD) below a fixed threshold ( ) [24][25][26]: A good approach is to apply different thresholds and see which one is most likely to represent the eruption [26,27]. The style of the eruption and the different weather conditions may have an influence on the threshold values. The higher the threshold, the higher the chance of having false alarms (i.e., ~ 0 K). In contrast, a lower threshold (i.e., < 0 K) increases the chance of underestimating the ash-contaminated areas. In terms of BTD thresholds, we use −0.2 K and −1.5 K for Kelud and Calbuco eruptions, respectively.
Since the clouds are located at high altitudes, they appear cold while observed at frequencies above 100 GHz [28]. When working at higher frequencies (>100 GHz), bare soil and oceans present higher BTs compared to the clouds' BTs. The identification of the volcanic cloud in the MW-MMW is a two-step identification method called microwave spectral difference (MSD). The first step is called the microwave spectral difference window ( ) [29]: where and are two frequencies. The former is in the range 155-165 GHz, the latter in the range 85-95 GHz. For the ATMS sensor, = 165.50 GHz and = 88.20 GHz. All pixels, which are below a given threshold, are detected as clouds. By convention, this threshold is set equal to 0 K [29]. Once all the cloud pixels are detected, the second step separates meteorological from volcanic clouds [1,29]. Meteorological clouds are characterised by a marked amount of water, which makes the channel around 183.31 GHz more sensitive due to the presence of the water vapour absorption peak. Volcanic clouds, on the other hand, respond better at lower frequencies. The formula below implements the microwave spectral difference absorption ( ), the second step of the method: where and are two different frequencies. For ATMS sensor = 183.31 ± 3 GHz and = 165.50 GHz. Volcanic plumes are separated from meteorological clouds by selecting the values below an arbitrary threshold, which is, by convention, set equal to 0 K [29]. False alarms can still occur also in this case. If the scene is mostly cloudy, fixing a threshold close to zero can increase the chance of having false alarms. Moreover, if scattered pixels are identified in the scene, these are excluded from the identification, since most aggregated pixels are more likely to belong to a proximal volcanic cloud [29]. In this way, only the most aggregated pixels are considered. In our analysis: the threshold is set to −9 K for Kelud eruption and to 0 K for Calbuco eruption; the threshold is set to 0 K for both eruptions. To summarise, it is not possible to assess that a given threshold can be useful to detect different volcanic clouds or even the same cloud at different instants mainly due to a change in the particle composition or different weather conditions, different times and different sensor settings [29,30].

Random Forest Classification Technique
In the last decade, the use of machine learning algorithms became popular also in atmospheric science. Most of the works present in the literature use the collected radiances only in the infrared to make the detection [31][32][33]. In addition, most of these models are designed to work with one specific sensor, even if some new studies are developing sensor-independent schemes [34]. The presented model instead combines different spectral information (i.e., MW-MMW and TIR) to attempt a general rule applicable on different sensors onboard of different platforms. The goal is to find a model that can detect the volcanic cloud automatically (i.e., without user action) and independently from the used sensor. The learned model is obtained by the following procedure: (1) definition of training images; (2) identification of the most relevant features, i.e., those that can better reflect the variability of the dependent variable; (3) model training; (4) model performance and selection; (5) classification of unseen images. The training data sample uses data collected by the MW-MMW and TIR sensors for the 2014 Kelud and 2015 Calbuco eruptions. Table  1 summarises the data and their application in the RF classification. The considered features are the common channels between the MHS-ATMS and AVHRR-VIIRS sensors plus the rectilinear distance (the sum of the absolute differences of the pixels' Cartesian co-ordinates) of each pixel from the volcano. More specifically, we use the channels 88.20 GHz (89.00 GHz for MHS), 165.50 GHz (157.00 GHz for MHS), 183.31 ± 1 GHz, 183.31 ± 3 GHz, 12.01 μm, 10.80 μm and 3.70 μm (3.74 μm for AVHRR). Each pixel is then treated as an independent sample during the training process. Due to the sensors' different pixel sizes, the MW-MMW observations were resampled on the TIR grid. In this way, each area is represented by the same number of pixels. The collected data are divided into sets, henceforth, train (80% of the pixels) and test (20% of the pixels) sets. The former is used to fit the model. The test is instead used to provide an unbiased evaluation of the model. The stratified cross-validation (CV) method was used to evaluate the actual model performance. With this method, we divide the training dataset into Kfolds (i.e., different k datasets of equal sizes where the percentage of each class does not change in each fold). On each iteration, a k th dataset is used as a validation set and the remaining k-1 folds are used as training sets. The error is calculated on the held-out fold (i.e., the k th set). This process continues up until all the k folds (10 in total) are used as validation sets. The k-fold CV result is computed by averaging all the validation errors [35]. During the training process, different bagging and boosting tree-based algorithms with different configurations are considered. All the different learned models are compared with each other by using the averaged CV score to find the one that better captures the underlying distribution. The F1 score is used as the evaluation metric: The F1 score allows false alarms and false negatives (FN) to be balanced and can be considered as the harmonic mean of the recall and the precision scores, assuming values between 0 and 1 [35]. The precision is computed as TP/(TP + FP) and the recall as TP/(TP + FN), where TP is the true positives and FP is the false positives. The precision measures the model capability of mislabelling as positive a negative pixel; the recall measures the model capability to correctly classify the positive pixels. The best considered model is the one learned by the RF algorithm with an F1 CV score of 0.8900 and of 0.9000 on our test. The RF configuration is summarised as follows: N° estimators, that is the number of trees in the forest, is set equal to 25; Criterion, which controls the function used to measure the quality of each split in each tree, uses the Gini index; Max depth, which controls the depth of the trees (i.e., how many splits in total), is set equal to 6; Max features defines how many features are considered when defining the best split and it is set equal to log N° features ; Class weight is used to associate weights to each class considering the class frequency and it is set equal to balanced; Bootstrap is set equal to true. The output of the model is a 0-1 binary mask, 0 for clear pixels (coloured in purple) and 1 for contaminated tephra pixels (coloured in orange). Once the mask is generated, we apply a nearest neighbourhood approach to identify and then exclude small clusters far from the volcanic crater. For the RF algorithm, we use the scikit-learn library version 1.0.2 (https://scikit-learn.org/stable/, accessed on 4 February 2023) and Python 3.7 version [36][37][38]. A schematic representation of the classification methodology is displayed in Figure 1. The output binary mask, where values of 1 (orange pixels) represent tephra pixels and values of 0 (purple pixels) represent clear (no tephra-contaminated) pixels.

Radiative Transfer Modelling
The radiative transfer model (RTM) scheme is designed to simulate the BTs observed by the TIR and MW-MMW sensors ( Figure 2). The contributions coming to the sensor are the background radiance (land or water), the atmosphere radiance and the space radiance (which is neglected during cloud physical properties retrieval [30]). We solve the radiative transfer differential equation by considering one-and two-layer approximations, with a normal pointing to the surface (i.e., we do not consider the inclination angle). We ignore the multiple scattering effect and the additional extension due to the electrostatic charge within highly explosive activity [39]. The one-layer model ( ) is applied to the TIR observations, as proposed in [24], whereas the two-layer model ( ) is applied to the MW-MMW observations. The ( Figure 2 left panel) is the solution to the radiative transfer differential equation considering the one-layer approximation, where, and are, respectively, the surface and the top cloud temperature and is the transmittance, expressed as a function of the optical depth τ, which is computed applying the following approximated formula: where is the geometric thickness of the cloud [24] and is the extinction coefficient, which is a function of the particle radius. The ( Figure 2 right panel) represents the solution to the radiative transfer differential equation considering the two-layer approximation. The contributions coming to the sensor are the , which is the surface contribution that is irradiated upward to the satellite and takes into account the surface temperature, the ground emissivity and the transmittance of the two layers; the , which is the first-layer brightness temperature that is irradiated upward; , which is the first-layer brightness temperature that is first reflected and then irradiated upward to the sensor; the , which is the second-layer brightness temperature that is first pushed downward and then it is irradiated upward to the sensor; and the is the secondlayer brightness temperature that is irradiated upward.
is the average temperature between and . and are the temperature of first and second layers, respectively, and , and are the respective layer heights. The analysis in the remote sensing framework is solved by inverting the forward problem. The forward problem starts from volcanic cloud properties based on the literature (e.g., particle size, density and concentration [2]), and then computes synthetic BT measurements. The synthetic BTs are then compared with the observed BTs to retrieve the main features of volcanic clouds. Our can be summarised into four main blocks; in the first block, we obtain the particle size distribution, in our case, a scaled Gamma. The second block estimates the cloud physical properties and the third block simulates the BTs considering the previously estimated layer features. Finally, the fourth block uses both the MLE and NN methods to estimate the total columnar content (TCC). We consider volcanic plumes composed by particles having particle density (ρ), effective radius ( ) and concentration . Effective radius ( ) and concentration ( ) values vary depending on particle classes. In this case, a positive relation is considered, i.e., larger particles correspond to larger concentration, as expressed in [2]. For the MW-MMW, we use the following classification: small lapilli (SL) having a particle density of 1200 kg/m 3 , in the interval 256-2048 µm and in the interval 10 3 -2 10 3 mg/m 3 [2,40]. For the TIR, instead, we consider the class fine ash (FA) with the following characteristics: particle density of 2600 kg/m 3 , in the interval 0.07-10.00 μm and in the interval 10 0 -10 1.5 mg/m 3 [2]. The complex refractive index (n) is called by the second block when cloud properties are retrieved. For MW-MMW, we use = 2.48 − i0.016 [40]. In the case of TIR, the refractive index and the extinction coefficient are derived by analysing results contained at Aerosol Refractive Index Archive (ARIA, University of Oxford, http://eodg.atm.ox.ac.uk/ARIA/, accessed on 4 February 2023). The used refractive indexes at 10.80 and 12.00 μm are, respectively, 2.10 + i0.41 and 1.79 + i0.19 for Kelud and 2.435 + i1.079 and 2.084 + i0.197 for Calbuco. The scattering of these particles is well described by Mie scattering theory. The volumetric extinction cross-section is calculated by considering the extinction cross-section ( ) and the particle distribution N(r) [41,42]: where the extinction cross-section is calculated as: with as the total extinction power, the incoming incident power density, the absorption cross-section coefficient and the scattering cross-section coefficient [41,42]. The particle size distribution N(r) is computed as: where is the intercept parameter, 2 is the volume-weighted median diameter and is the shape parameter of the gamma distribution [40]. The surface and the cloud temperatures, called by the fourth block, are used to compute the synthetic BTs with the surface emissivity set to 0.90 for MW-MMW [29]. The well-known radiative transfer differential equation is solved by considering the one-and two-layer approximations [41,42]: where is the absorption coefficient, is the albedo, is the emissivity, 0, and , are, respectively, the entire atmosphere optical thickness and the optical thickness from layer to the top of the atmosphere, is the surface temperature and is the physical temperature of the medium. Figure 3 shows the block diagram of the developed . . Block diagram of MW-MMW and TIR retrieval procedure. The first orange block computes the particle size distribution (PSD) using the particle density, effective radius and concentration. The second block is involved to retrieve cloud physical properties. It first computes extinction crosssection (based on Mie Theory) considering refractive index and wavelength, then the volumetric extinction cross-section given PSD and . Finally, the optical depth is estimated by considering and geometric thickness (l). The third block simulates the BTs considering the surface temperature, the layer temperature ( ), the PSD, the optical depth and the emissivity. The is the simulated BT at a given wavelength. The fourth block performs the retrieval using two methods: the MLE and NN. Figure 4 shows the solution of the plotted in arch curves. For both MW-MMW and TIR curves, we return the associated effective radius. For MW-MMW we show the effective radius interval that goes from 0.28 to 2.00 mm. Instead, for TIR measurements, we show the effective radius interval from 1.50 to 10.00 μm. The TCC or mass loading identifies the superficial density of the volcanic cloud. For that purpose, the cloud is treated as a cylinder and the TCC is computed by the integral of the ash content into a cylinder of unit area as [24,25]: where is the ash cloud density, is the particle radius and is the height. The total mass is then estimated by multiplying (10) for the area of each pixel and aggregating them. The pixel area is calculated from the image itself by looking at the latitude and longitude information. This allows a higher level of flexibility since the estimate is independent from the sensor and the scanning angle. The uncertainty related to the TCC estimate is given by: where we assume the following values for each independent parameter: δr = 0.20r and δH = 0.30H. The relative percentage error is equal to 36%. For each simulated BT, different information is known, such as , and . For observed pixels these quantities are retrieved by associating each of them to their closest curves. In this work, we compare two methods to pair synthetic BTs with observed BTs: the maximum a posteriori probability (MAP) and an NN architecture. Assuming uniform prior information, the MAP and the MLE converge to the same solution. Moreover, assuming Gaussian-likelihood statistics of the difference between the synthetic BTs and observed BTs and uncorrelation between the deviations (i.e., errors) of the same observables, the MLE reduces to the minimisation of a quadratic form [2]. The MLE method is used to retrieve information about and TCC, computed by the , for the collected BTs.
The NN architectures are particularly powerful to solve nonlinear problems. In this work, we consider an architecture composed by one hidden layer. This configuration returns good results to solve this multiple-nonlinear regression problem. The NN fitting procedure is involved in defining the unknown parameters that fit the training data well. In NN, the unknown parameters are the weights θ. Since it is a regression problem, the root mean squared error (RMSE) is used as a measure of fit. The stochastic gradient descent optimisation method was used to minimise the mean squared error (MSE) loss function [14,35]. This speeds up the computation since the gradient is not computed on the whole dataset but on a randomly selected subset. More specifically, the adaptive moment estimation (ADAM) method is used [43]. The learning rate for the batch size is optimised by a linear search that minimises the error function at each iteration. This means that goes to zero while increasing the number of iterations r. ADAM uses the first and the second moments of the gradients to adaptively adjust the learning rate. The weights are randomly initialised with values near to zero. In this way, the model becomes nonlinear as soon as the weights are updated [14]. The response variables TCC and are dependent upon the BTs and upon each other. This requires a model which is able to predict both quantities at once. NN allows a multi-output model to be defined. For that reason, the NN architecture presents two nodes in the output layer, one per each response variable. The input layer, instead, is composed by n numbers of nodes based on the number of simulated BTs. Mathematically, the used architecture can be expressed as follows [14]: where is the output function, is the weight matrix , a is the activation function, is the input (simulated BTs) and the bias matrix . The dataset used to train and evaluate the model is composed by the BTs simulated with the . In particular, for each simulated channel, there is a 500 × 500 matrix where each entrance is a simulated BT expressed in K. Then, the dataset is divided into train (80%) and test (20%) sets. The CV method, applied on the train set, was used to evaluate the actual NN performance by using the RMSE metric. The NN architecture is summarised in Figure 5. Since we are working with both the MW-MMW and TIR sensors, we developed two NN models. The NN model for the MW-MMW observations has three input nodes represented by the synthetic BTs at frequencies 88.20, 165.50 and 183.31 GHz. The hidden layer is composed of 256 nodes. The batch size is set to 512. The number of epochs is set to a large number with the early stopping method. In this way, the training stops when the model performance stops, improving after a given number of iterations. The rectified linear unit (ReLu) activation function is used in the hidden layer [14,44,45]. The NN architecture for the TIR has instead two input nodes: BTs at 10.80 and 12.01 µm. The hidden layer is composed of 128 nodes. The batch size and number of epochs are treated as for the MW-MMW network. Moreover, here, the ReLu activation function is used in the hidden layer. For the NN implementation, we use the TensorFlow 2.8.2 open-source platform (https://github.com/tensorflow/tensorflow, accessed on 4 February 2023) and the Keras 2.8 deep learning API (https://github.com/keras-team/keras, accessed on 4 February 2023) written in Python [46].

The Empirical Parametric Retrieval (EPR) Method
The MW-MMW mass loading estimates, obtained with the MLE and NN methods, are compared with the parametric regressive formula named empirical parametric retrieval (EPR) [1]. The EPR is used as a quality measure of the MLE and NN estimates. The EPR comes from a regression analysis applied on a numerical model simulation in terms of MW BTs of volcanic plumes over land [1]: where = 63.84, = −0.2564 and = = 2500 kg/m 3 . The is the brightness temperature around the 183 GHz channel, the frequency near the water vapour absorption peak. The tephra mass loading estimate ( ) is almost not dependent on the surface land emissivity variation due to the use of the absorption channel around 183 GHz [1].

Test Cases
In this section we describe the Kelud and Calbuco eruptions conjointly with an explorative data analysis of the ATMS and VIIRS products.

The Kelud Eruption in 2014
Kelud volcano is one of the most dangerous stratovolcanoes in Indonesia. It is located at −7°55′48.00″S 112°18′28.80″E, East Java. At the beginning of 2014, the earthquake activity started increasing, alerting the local community of an awakening of the volcano. On 13th February 2014 at 15:46 UTC, the seismic activity stopped, indicating the onset of the eruption. The first recorded explosion formed a high eruption column. The top of the umbrella cloud's height was estimated at approximately 20 km, with a diameter of more than 200 km [47]. The plume drifted south-west across the Java Island and the Indian Ocean [48]. The eruption caused many problems to air traffic; indeed, many local flights were cancelled or rerouted. The tephra ashfall damaged schools, homes and local businesses. More than 76,000 people were evacuated from their homes until the eruption activity stopped on 17th February [47,48]. Figure 6 shows the explorative analysis from the ATMS and VIIRS sensors. The first row of Figure 6 shows the channels 88.20, 165.50, 183.31 ± 4.50 and 183.31 ± 1 GHz. Over the land, the use of frequencies above 80 GHz displays how the scattering is stronger if compared with the emission-based process. The background BT signatures at frequencies above 90 GHz appear to be hotter (yellow pixels) than the cloud ones (blue and dark blue pixels). This allows the visual detection of volcanic cloud pixels. The second row represents the channels 12.01, 10.80, 8.55 and 3.70 µm of the VIIRS sensor. The presence of meteorological clouds coupled with shorter wavelengths (compared to MW-MMW) does not allow full discrimination between the volcanic cloud and the background. Indeed, by looking only at channel 3.70 µm, both background and the contour of the cloud appear blue. The other TIR channels better distinguish the cloud and the background (yellow areas). However, the borders of the clouds are still not well separated by the meteorological clouds (light blue pixels around 240 K).

The Calbuco Eruption in 2015
Calbuco is an active and hazardous stratovolcano located at 41°20′S-72°37′W in southern Chile. It covers an area of 150 km 2 , with its summit at 2003 m above sea level. After 54 years of calm, it reawakened on 22nd April 2015 at 21:08 UTC. Alarms of the volcano's awakening were signalled by an intense seismic activity [49]. On 23 April 2015 at 04:00 UTC, a strong explosion occurred, which lasted about six hours. The ash plume, with a maximum estimated altitude of about 15-20 km, moved in the north-east direction [1]. On 24th April 2015 at 02:30 UTC, a smaller explosion occurred. During those events, different cities registered inconveniences and mobility problems. The air traffic in the airport of Puerto Montt was also suspended. Different cities and areas were affected by tephra fallout, such as Los Lagos, Los Rios and Araucani [1]. Figure 7 is an explorative analysis from the ATMS and VIIRS sensors. Moreover, in this case, the first row of the figure shows the channels 88.20, 165.50, 183.31 ± 4.50 and 183.31 ± 1 GHz. The volcanic cloud pixels appear colder while moving at frequencies above 90 GHz, due to the increasing ratio between particle sizes and wavelengths. In particular, frequencies around 165.50 GHz and 183.31 GHz allow a clear discrimination between the background and the volcanic cloud. The coldest pixels (i.e., darkest blue and black) highlight the area covered by the volcanic cloud, whereas yellow pixels the background. It is interesting to notice how, using the ATMS sensor, the distal cloud is not visible in any of these channels. This is because these frequencies are not sensitive to smaller dispersed ashes. The channels 12.01, 10.80 and 8.55 µm show in dark blue and black the near-source cloud, in light blue the dispersed cloud and in yellow the background. The channel 3.70 µm allows a good visual representation of the volcanic cloud at the source (darkest blue and black pixels) but, due to its sensitivity to finest particles, the volcanic distal cloud may be partially masked by the presence of the other background atmospheric effects.

Results
In this chapter, we show the detection and the retrieval using MW-MMW and TIR spectral signatures. We first describe the detection results obtained using the split window methods and the RF model, then the estimates obtained using the MLE and the NN methods. For the detection, the results are displayed with three images, where images A and B show the detection using the MSD and BTD methods, respectively; image C shows the detection using the RF model. For the RF outputs, the results are plotted using three different colours. Regarding the retrieval, we provide simulated arch curves, representing the simulated BTs, for both spectral signatures. The TCC maps are grouped into two figures composed of five panels each. Section 4.1 is dedicated to the Kelud event and Section 4.2 to the Calbuco event.

Kelud: Detection and Retrieval
For the Kelud eruption, the MSD cloud pixel identification is visible in Figure 8A. The is higher for pixels which are proximal to the volcanic crater. This effect is shown in this image by the pixels with the yellow shade colour (higher ) and dark blue pixels (smaller ). The BTD detection is displayed in Figure 8B. The reported values are already adjusted with the water vapour correction explained in Appendix B. The edges of the cloud are characterised by the presence of pixels with a BTD smaller than zero (pale-blue and light-yellow colours), whereas the pixels near the volcano are characterised by a BTD slightly above zero (yellow pixels), where the abundance of water particles is relevant. Figure 8C shows the output of the RF model. As a measure of comparison, the images A and B are overlapped and used as ground truth. In this way, every RF-classified pixel is compared with its symmetrical available in the A + B combined image. The pixel is coloured in purple if it is identified by both the methods, the MSD-BTD and the RF; it is magenta if it is identified by the RF but not by the MSD-BTD; it is orange if it is identified by the MSD-BTD but not by the RF. Almost 90% of the area is identified by the RF model (purple pixels). A small fraction is missing (i.e., orange pixels) and another small part disagrees with the MSD-BTD detection (magenta pixels). The overall F1 score on this image is 0.9049, which is a very high score and higher than the one obtained during the test phase, meaning that the model is performing correctly and it is not overfitting. Regarding the retrieval, Figure 9A details the simulated curves for MW-MMW volcanic particles. The black points represent the volcanic particle pixels previously detected using the MSD method. These points are associated with particles of effective radius that vary from 0.28 to 0.60 mm. The points that lie outside the simulated curves are associated with their closest curve [24]. For the TIR curves, the detected ash particles have a particle radius that ranges between 5.00 and 10.00 μm. These points are represented in Figure 9B as orange points. The black points represent the MW-MMW-detected pixels, which were resampled on the TIR grid (this is the reason why there are more than those displayed in Figure 9A). As expected, many points lie outside the simulated curves and have a positive BTD. These are the points that represent the areas contaminated by particles not visible at the TIR wavelengths due to sensor detection limits. The simulated arch curves for TIR observations. The black points represent ash pixels identified using the MSD method. The orange pixels highlight the pixels detected using the BTD method.
Starting from these simulated curves, we perform the mass retrieval. Per each of the three methods, the MW-MMW masses are, respectively, EPR 1.54 ± 0.55 10 11 kg, MLE 1.11 ± 0.40 10 11 kg and NN 1.30 ± 0.47 10 11 kg. All the three methods highlight how the higher masses are estimated to be in the centre of the cloud (yellow pixels of Figure 10A-C). The smaller masses are associated with the areas far away from the volcanic crater (blue pixels). The instantaneous total masses and the TCC maps, for both the MLE and the NN, are in good agreement with the one obtained using the parametric formula (EPR). Figure 10 D, E show the estimated TCC maps for both MLE and NN methods considering TIR observations. The estimated masses are, respectively, MLE 2.25 ± 0.81 10 9 kg and NN 2.16 ± 0.78 10 9 kg. The MLE and the NN TIR results show the highest mass loading in the centre of the cloud (yellow pixels). However, the NN method shows a more homogeneous path in the map, with high mass loading areas also outside the centre of the cloud. This is more in line with the MW-MMW retrievals, in particular, the EPR ( Figure 10A) and the NN ( Figure 10C). Indeed, Figure 10 A, C, E are in agreement with each other (see yellow and blue pixels) in terms of mapping. Then, the NN results (MW-MMW and TIR) better characterise the volcanic cloud. Moreover, our NN estimate of the mass, given by the sum of the MW-MMW and TIR estimates, is 1.32 ± 0.47 10 11 kg, which is in line (in terms of order of magnitude) with volcanological studies. Following [48,50], using field data collected on ground, a total erupted mass was estimated to be 3.30-6.60 10 11 kg. It is noteworthy that we are observing the eruption only at the time of the satellite overpass (17:28 UTC) when the eruption was still ongoing. This implies that we are not considering in the comparison the mass emitted after the satellite acquisition time until the end of the eruptive event.

Calbuco: Detection and Retrieval
For the Calbuco eruption, the MSD cloud pixel identification is displayed in Figure  11A. The assumes high values in almost the whole detected area, except for one pixel with a value of ≈ −11 K. When the plume is spreading, it has a greater amount of fine particles, not detected by the MSD. The presence of smaller particles is, instead, well visible in Figure 11B (see outlines coloured in pale blue). As for Kelud, the reported BTDs values are already adjusted with the water vapour correction described in Appendix B. Figure 11C shows the output of the RF model. Moreover, here, the images A and B are overlapped and used as ground truth. Most of the pixels are coloured in purple, with an area coverage of about 90%. A small fraction is missing (i.e., orange pixels) and a very small fraction is wrongly classified (magenta pixels). The overall F1 score is 0.9271, which is even higher than the one obtained for the Kelud event ( Figure 8C). This is because the fraction of false alarms is very small compared to the fraction visible in Figure 8C.   Figure  12A displays the simulated arch curves for MW-MMW observations. The black points identify the detected MW-MMW pixels. The particle radius ranges in between 0.28 and 0.85 mm. For the TIR observations we show the results in Figure 12B. Moreover, here, the black points highlight the MW-MMW detected pixels which were resampled on the TIR grid. As expected, many of these points lie outside the simulated TIR curves. At the time of the pass of the Soumi-NPP satellite, the volcano was erupting for the second time. Indeed, the VIIRS sensor also detected a first event (22 April 2015 at 21:08 UTC) represented in Figure 12B by the yellow points. The orange points, instead, represent the pixels belonging to the near-source cloud (i.e., the second event at 05:09 UTC), the one detected also by the ATMS sensor. It is visible how the two events, observed in the TIR region, have particles with different sizes (from 2.00 to 4.00 μm for the first event and from 2.50 to 5.00 μm for the second event). The reason is that bigger particles remain suspended in atmosphere for less time due to their bigger mass. Associating each observation to its closest curve, the retrieval is, hence, performed. The MW-MMW masses obtained with the considered methods are EPR 4.21 ± 1.51 10 10 kg; MLE 4.33 ± 1.56 10 10 kg; and NN 4.17 ± 1.50 10 10 kg, respectively. For all three methods, the higher mass is estimated on the west side of the volcano (yellow pixels), in proximity to the top of the umbrella at an altitude of ≈ 19 km [1], see Figure 13A-C. The smaller mass is associated with the areas far from the volcanic crater (bright blue pixels). The estimated mass and the TCC maps, for both the MLE and the NN, are in good agreement with the one obtained using the parametric formula (EPR). The small changes between the estimates can be explained with the different polarisation corrections and the ignored multiple diffusion effects. For the TIR observations, the estimated masses are, respectively, MLE 1.48 ± 0.53 10 9 kg and NN 1.54 ± 0.55 10 9 kg. For consistency, we also report the estimates for the distal ash cloud, which are 2.37 ± 0.85 10 9 kg for MLE and 2.17 ± 0.78 10 9 kg for NN. The TIR proximal and distal masses are reported in Figure 13 D, E with the acronyms and . For both MLE and NN TIR results, the highest mass loading is associated with the centre of the cloud (yellow pixels) with values in the range of 0.06-0.08 kg/m 2 , see Figure 13 D, E. However, the NN TCC map ( Figure 13E) is more in line with the EPR map ( Figure 13A). Indeed, Figure 13 A, C, E are in good agreement with each other (see yellow and blue pixels) in terms of mapping. Moreover, in this case, we conclude that the NN results are more in line with the EPR results and past studies [1,49]. If we sum the NN MW-MMW and TIR mass estimates, we obtain a mass of 4.32 ± 1.56 10 10 kg that is in line with previous studies on the deposit (3.33-4.71 10 10 kg) [49].

Conclusions
The aim of the work was to show the synergy given by the combination of the MW-MMW and TIR observations to study volcanic clouds. The application on the Kelud and Calbuco case studies highlighted how the combination of these two spectral signatures can increase the quality of the detection and the mass retrieval. The proposed methodologies allow the study of the time evolution in terms of scene coverage and mass of the volcanic cloud.
The RF model allows images to be automatically processed without the arbitrary threshold choice and has, in our opinion, promising results. Moreover, our results underline how it is possible to combine different spectral information, coming from different platforms, to have a better characterisation of the volcanic cloud. The RF model is able to balance false alarms and the FN while classifying new images (see Figures 8C and 11C). Indeed, for the Kelud eruption the F1 score is 0.9049 and for the Calbuco is 0.9271. It is well known that the MW-MMW signatures between 90 and 183 GHz are more effective in identifying the near-source plume, even though they have a poorer spatial resolution, up to 32 km compared with 750 m of TIR observations. Instead, the TIR signatures perform better in identifying the edges of the proximal cloud, since, close to the crater, the signal tends to saturate due to the high optical extinction of the near-source volcanic cloud and the strong presence of water particles. The areas detected by both sensors (ATMS and VIIRS) highlight how some areas are contaminated by both bigger and finer volcanic particles. Further improvements in terms of detection should be considered, especially for the machine learning method. The RF results are promising but collecting more heterogeneous data should make the training dataset more statistically representative for this task. We know that using only a few images can affect the performance of the model. However, it is not easy to find cases where both MW-MMW and TIR signatures are available with a small difference in the time span. The absence of GEO MW-MMW sensors makes the acquisition even harder. It is also worth mentioning that this was a first attempt to prove the possibility of combining different information without complex data manipulation (polarisation and scanning angle corrections). Furthermore, the RF model can be improved by moving from a binary to a multiclass classification problem by considering, for example, "ash on meteorological cloud", "ash on land" and "ash on water" as already tested for infrared signatures [31,34]. This could increase the performance of the model, in particular when the scene is significantly contaminated by meteorological clouds, as it is in the Kelud event.
Concerning the retrieval, further improvements could also be achieved by making the radiative transfer model more realistic, for example, by considering the cloud as a multiplanar nonhomogeneous layer [51]. The mass loading and the mass estimates are valid until the particle spherical assumption holds. Indeed, volcanic ash particles are usually non-spherical [24]. Further works should also start to consider the particles of irregular shapes, since different shapes give different types of scattering [52]. Even so, the nonspherical assumption is still under investigation, making this assumption a fair and conservative compromise to date in real-case applications [53,54]. In addition, the electrical charging can influence the spectral signature during the first instants of an eruption [10,39,55,56]. In this work, we are not considering the additional extinction caused by electrostatic charging, but we plan to study how this component attenuates the signals in future studies. Moreover, a single NN architecture could be implemented in the future to return estimates for both MW-MMW and TIR observations simultaneously. Lastly, we are also considering combining satellite estimates with ground-based radar measurements and lidar observations to improve the retrieval in terms of mapping and mass estimates. Funding: This project has received funding from the PON-RAFAEL project of CETEMPS, University de L'Aquila. This work has also benefited from the e-shape project that receives funding from the European Union's Horizon 2020 research and innovation programme under grant agreement 820852.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author and first author.

Acknowledgments:
The authors would like to thank NOAA and EUMETSAT for providing MW-MMW and TIR satellite data and the ERA5 atmospheric data provider. The authors also thank Stephen Conway, an English native speaker reading the paper for INGV-OE.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript or in the decision to publish the results.

Appendix A
This appendix has the scope of summarising all the acronyms (see Table A1).

Appendix B
This appendix deals with the water vapour correction for the TIR observation applied in our estimations. In the 8-14 µm range, the silicate particle absorption is greater around 8-11 µm, whereas the water particle absorption is greater at higher wavelengths (12-14 µm) [25,57]. The volcanic clouds are mainly composed of silicate particles but, if the presence of water particles is high in the atmosphere, the BT of the top atmosphere is higher at wavelength 11 µm than at 12 µm. This has the side effect of making the BTD positive. Prata et al. [25,30] proposed a semi-empirical water vapour correction: where BT* = BT10.8/BTmax. The denominator is a normalisation constant set at 320 K. The free parameter b governs the water vapour effect on the BTD at the maximum value of . . The free parameter is estimated as the difference of the maximum value of . and its corresponding value at [26]. The water vapour correction for the TIR observations is estimated for the Kelud and the Calbuco eruptions ( Figure A1). The images A and C show the count corrections, whereas the images B and D show the lower and upper water vapour limit curves, as expressed in [26,30], where the points represent all the pixels associated with the volcanic clouds. The water vapour correction is applied for all the points that lie inside the two curves. The water vapour correction allows to detect, by using the BTD, 8982 more ash pixels for the Kelud near-source cloud (i.e. 30.57% more of the area) and of 5104 more pixels for the Calbuco near-source cloud (i.e. 33.72% more of the area).