Classifying Major Explosions and Paroxysms at Stromboli Volcano (Italy) from Space

: Stromboli volcano has a persistent activity that is almost exclusively explosive. Predomi-nated by low intensity events, this activity is occasionally interspersed with more powerful episodes, known as major explosions and paroxysms, which represent the main hazards for the inhabitants of the island. Here, we propose a machine learning approach to distinguish between paroxysms and major explosions by using satellite-derived measurements. We investigated the high energy explosive events occurring in the period January 2018–April 2021. Three distinguishing features are taken into account, namely (i) the temporal variations of surface temperature over the summit area, (ii) the magnitude of the explosive volcanic deposits emplaced during each explosion, and (iii) the height of the volcanic ash plume produced by the explosive events. We use optical satellite imagery to compute the land surface temperature (LST) and the ash plume height (PH). The magnitude of the explosive volcanic deposits (EVD) is estimated by using multi-temporal Synthetic Aperture Radar (SAR) intensity images. Once the input feature vectors were identiﬁed, we designed a k-means unsupervised classiﬁer to group the explosive events at Stromboli volcano based on their similarities in two clusters: (1) paroxysms and (2) major explosions. The major explosions are identiﬁed by low/medium thermal content, i.e., LSTI around 1.4 ◦ C, low plume height, i.e., PH around 420 m, and low production of explosive deposits, i.e., EVD around 2.5. The paroxysms are extreme events mainly characterized by medium/high thermal content, i.e., LSTI around 2.3 ◦ C, medium/high plume height, i.e., PH around 3330 m, and high production of explosive deposits, i.e., EVD around 10.17. The centroids with coordinates (PH, EVD, LSTI) are: C p (3330, 10.7, 2.3) for the paroxysms, and C me (420, 2.5, 1.4) for the major explosions.


Introduction
Stromboli volcano is a small island (less than 5 km wide) in the southern Tyrrhenian Sea (Italy), which has been continuously erupting for the past 2000 years [1]. Its activity is almost exclusively explosive, but lava flows do occur at times, most recently in 2014 and 2020 [2,3]. Its persistent but moderate explosive activity, termed "Strombolian", is occasionally interrupted by explosive events that are more violent and represent the main hazard for the inhabitants of the island (less than 500 in winter, but more than 5000 in summer). The products of the most powerful explosive eruptions can damage the buildings [4] and even cause extensive vegetation fires [5]. Nevertheless, fatalities and injuries may also occur, such as in the summer 2019 [6]. To assess the hazards faced by communities living on the island, it is fundamental to detect the explosions of different magnitude and intensity and understand how to classify them. The classification being a paramount issue for a correct hazard assessment [7], in the last decades several classifications have been proposed. Barberi et al. [8] distinguished three types of eruptive activity on a growing scale: (1) a persistent, ordinary explosive activity of low intensity, (2) major explosions and (3) paroxysmal events.
The persistent activity at Stromboli is characterized by several hundreds of short explosive transients (<10 s) per day erupting small ash columns (<100 m) that normally have no impact on the local population [7]. Conversely, the more intense and powerful explosions correspond to greater magnitude (i.e., erupted volume; [9]) and intensity (i.e., instantaneous effusion rate; [10]), resulting in longer explosive transients and higher eruptive columns [7,11].
Generally, more than one crater in the summit area of Stromboli is active during major explosions, and these produce eruptive columns a few hundred meters high, with fallout of bombs and ash up to a maximum distance of about 1 km from the crater area. In most cases, the proximal area is affected by ballistic fallout of lithic blocks (up to decametric sizes), impacting the Pizzo Sopra La Fossa location (Figure 1) where tourists gather to observe the activity [8,12,13]. The volume of products ejected range from 10 2 to 10 3 m 3 of tephra (fragmental pyroclastic materials that fall to the ground from eruption columns) [14,15].
Remote Sens. 2021, 13, x FOR PEER REVIEW 2 o paramount issue for a correct hazard assessment [7], in the last decades seve classifications have been proposed. Barberi et al. [8] distinguished three types of erupt activity on a growing scale: (1) a persistent, ordinary explosive activity of low intens (2) major explosions and (3) paroxysmal events. The persistent activity at Stromboli is characterized by several hundreds of sh explosive transients (<10 s) per day erupting small ash columns (<100 m) that norma have no impact on the local population [7]. Conversely, the more intense and powe explosions correspond to greater magnitude (i.e., erupted volume; [9]) and intensity ( instantaneous effusion rate; [10]), resulting in longer explosive transients and hig eruptive columns [7,11].
Generally, more than one crater in the summit area of Stromboli is active dur major explosions, and these produce eruptive columns a few hundred meters high, w fallout of bombs and ash up to a maximum distance of about 1 km from the crater area most cases, the proximal area is affected by ballistic fallout of lithic blocks (up decametric sizes), impacting the Pizzo Sopra La Fossa location (Figure 1) where tour gather to observe the activity [8,12,13]. The volume of products ejected range from 10 10 3 m 3 of tephra (fragmental pyroclastic materials that fall to the ground from erupt columns) [14,15]. Paroxysms are characterized by eruptive columns rising a few km above the cra and causing fires and damage to the populated villages on the lower flanks of the volca 2.0-2.5 km away from the summit craters [8,[16][17][18][19]. These large-scale events typic consist of a cluster of explosions involving all active vents and ejecting meter-sized sco bombs and blocks up to 2 km from the craters, and by greater volumes (from 10 4 to m 3 ) of emitted materials [8,20,21]. Usually, paroxysms cause considerable changes to crater area, leaving a wider and deeper crater depression [8]. Furthermore, lava flows m be emitted right before [6,11], after [17], or during a paroxysm [22]. Often, hot avalanc or pyroclastic density currents (PDCs) caused by the opening of flank fissures [23], by collapse of eruptive columns during paroxysms [6], by the collapse of small portion the summit cone [24,25], by flank failure [26] or by the brecciation of lava flow fronts [ spread along the steep and unstable Sciara del Fuoco (SdF; Figure 1) slope and over sea surface [3,25]. More rarely, PDCs may affect the other slopes of the island [8,27].
Phenomena observed during Strombolian paroxysms may vary case by case [7]. instance, large volumes (~10 5 m 3 ) of PDC are often ejected along the Sciara del Fu during paroxysms, whereas they are normally absent during major explosions, and thi Paroxysms are characterized by eruptive columns rising a few km above the craters and causing fires and damage to the populated villages on the lower flanks of the volcano, 2.0-2.5 km away from the summit craters [8,[16][17][18][19]. These large-scale events typically consist of a cluster of explosions involving all active vents and ejecting meter-sized scoria bombs and blocks up to 2 km from the craters, and by greater volumes (from 10 4 to 10 5 m 3 ) of emitted materials [8,20,21]. Usually, paroxysms cause considerable changes to the crater area, leaving a wider and deeper crater depression [8]. Furthermore, lava flows may be emitted right before [6,11], after [17], or during a paroxysm [22]. Often, hot avalanches or pyroclastic density currents (PDCs) caused by the opening of flank fissures [23], by the collapse of eruptive columns during paroxysms [6], by the collapse of small portions of the summit cone [24,25], by flank failure [26] or by the brecciation of lava flow fronts [23], spread along the steep and unstable Sciara del Fuoco (SdF; Figure 1) slope and over the sea surface [3,25]. More rarely, PDCs may affect the other slopes of the island [8,27].
Phenomena observed during Strombolian paroxysms may vary case by case [7]. For instance, large volumes (~10 5 m 3 ) of PDC are often ejected along the Sciara del Fuoco during paroxysms, whereas they are normally absent during major explosions, and this is related to the greater heights reached by the eruptive column during paroxysms [28].
However, the threshold level distinguishing these categories is not fixed and has not been addressed so far. Conversely, lava flows may or may not be emitted before or after a paroxysm, thus being a process not always associated with paroxysmal events. Although the island has a modern and robust instrumental monitoring system, it does not provide a robust way to identify automatically the different sizes of the most intense explosions. For example, the typology of the most powerful explosive eruptions occurring at Stromboli in the period 2018-2020 have been distinguished between major explosions and paroxysm events using as discriminant factors either the area affected by large ballistic projectiles [29] or the magnitude of seismicity combined with the ground deformation, muzzle velocity, and height of the ash plume [7]. It is worth noting that some of the events were not uniquely defined (such as the 19 July 2020 episode) and a generally accepted classification combining all measurements has not been established yet.
Nowadays, a variety of sensors on board satellites is widely used to monitor volcanic activity from space allowing to characterize different volcanic phenomena, from lava flow mapping [30] to ash cloud tracking [31]. Moreover, recent improvements in the frequency, type and availability of satellite images mean it is now feasible to use a machine learning multi-parametric approach by taking into account different satellite data to obtain a reliable classification of the more powerful explosions. Here, we seek to address this task by integrating infrared and radar satellite data. It is worth noting that the detection capability of the satellite measurements themselves only provides a piece of information about the monitored phenomena. Thus, we have investigated the most powerful events occurring from January 2018 to April 2021 at Stromboli volcano by examining the land surface temperature (LST) time series derived by the daily Moderate Resolution Imaging Spectroradiometer (MODIS) LST product, the variable magnitude of volcanic deposits retrieved by using Sentinel-1 SAR amplitude data, and the maximum height reached by the eruptive column estimated by Spinning Enhanced Visible and Infrared Imager (SEVIRI) and Visible Infrared Imaging Radiometer Suite (VIIRS) data. Furthermore, we show that the satellite derived LST information gives an insight into the pre-eruptive phase of some of the investigated paroxysmal events. We designed a k-means unsupervised classifier, exploiting the functionalities of Google Earth Engine (GEE), which is a cloud computing platform for environmental data analysis from local to planetary scales, with fast access and processing of satellite data from different missions [32].

Record of the Most Powerful Explosions
Stromboli, the northeasternmost island of the Aeolian Archipelago (Tyrrhenian Sea, southern Italy; Figure 1), is a stratovolcano that extends to a depth of about 2000 m below sea level. Only the summit of the volcanic edifice emerges above sea level, reaching an elevation of 924 m [33]. The persistent explosive activity originates from several vents nested within the Crater Terrace, a small depression located at about 750 m of elevation on top of the Sciara del Fuoco. The Crater Terrace comprises three areas where a variable number of eruptive vents open (Figure 1): the NE Crater zone (NEC), Central Crater zone (CC) and SW Crater zone (SWC). Typical Strombolian explosions eject pyroclastic fragments to a height of a few tens of meters, which fall at a short distance from the eruptive vent.
Between January 2018 and April 2021, the mild Strombolian activity was interspersed with nineteen high-intensity explosive events producing fallout of heavy materials (blocks and pyroclastic material) on the upper part of the volcano and, occasionally, on the villages of Ginostra and Stromboli (Figure 1). These short-lived, violent eruptions of variable scale, were identified as either major explosions or paroxysmal events depending on their size, covering a wide spectrum in terms of recurrence, intensity, and magnitude [7,29].
The historical catalog of the last 140 years includes 180 higher energy explosions, conventionally classified as major explosions and paroxysms, with an average frequency ranging from 0.7 to 2.1 events per year [29]. However, the current annual rate of higher energy explosions is 2.8 events per year based on the last 10 years, and increased up to Table 1. List of the powerful eruptive events from January 2018 to April 2021 with the classes derived from Calvari et al. [7], Bevilacqua et al. [29] and the INGV weekly bulletins. Major Explosion (ME), Paroxysm (PA), Not Available (NA) in the case no information is provided. Date and time in dd/mm/yyyy hh:mm:ss format. Time is expressed in UTC.

Satellite Data Sources
The satellite products used in this study are based on the outputs of the LST product from MODIS, radiometric measurements from SEVIRI and VIIRS, and microwave observations from Sentinel-1.
The Moderate Resolution Imaging Spectroradiometer (MODIS) is a key instrument on board of both satellites of Terra and Aqua missions of the NASA-centered international Earth Observing System (EOS). Terra and Aqua MODIS are viewing the entire Earth's surface every 1 to 2 days, acquiring images in 36 spectral bands. We processed the level 2 LST and Emissivity product in GEE. In particular, Emissivity in bands 31 and 32 is estimated from land cover types and in the day/night algorithm, daytime and nighttime LSTs and surface emissivity are retrieved from pairs of day and night MODIS observations in seven TIR bands.
The Spinning Enhanced Visible and Infrared Imager (SEVIRI) is the main payload of the current Meteosat Geostationary (MSG) satellite operated by EUMETSAT. SEVIRI provides radiances within 12 spectral channels covering the visible to the infrared [34] with a 3 × 3 km nominal spatial resolution at the subsatellite point. The primary location of Meteosat is at 0 • E above the equator, and SEVIRI has a temporal resolution of 15 min (four full disk images per hour).
The Visible Infrared Imaging Radiometer Suite (VIIRS) instrument is aboard the joint NASA/NOAA Suomi National Polar-orbiting Partnership (named Suomi NPP) satellite and the joint NASA/NOAA NOAA-20 satellite. The VIIRS Level 1 and Level 2 swath products are generated from the processing of 6 minutes of VIIRS data acquired during the NOAA-20 satellite overpass. The NOAA-20 product comprises 5 bands that are sensitive to visible, near-infrared (NIR), and thermal wavelengths. The spatial resolution of the instrument at viewing nadir is approximately 375 m for the I-bands, and 750 m for the DNB and the M-bands.
The Sentinel-1 mission of the Copernicus program of the European Space Agency (ESA) consists of two identical satellites equipped with C-band Synthetic Aperture Radar (SAR) instruments (frequency 5.33 GHz, wavelength 5.6 cm). The images are provided in dual or single polarizations, in four different operational modes, namely Stripmap (SM), Interferometric Wide swath (IW), Extra Wide swath (EW) and Wave (WV) modes. The spatial resolutions vary between 10 m and 40 m, depending on the selected mode. Sentinel-1 SAR images are made available at level-1C Ground Range Detected (GRD) in GEE, where they are already preprocessed based on the algorithms implemented by the Sentinel-1 Toolbox software [35].

Methods
Strombolian explosions of different magnitude and intensity do not fit easily into traditional styles of eruption [28] and the differentiation between paroxysms and major explosions on the basis of their intensity derives from measurements that are not always replicable (i.e., field surveys), so their classification is somewhat problematic and suffers from human subjectivity [7]. In recent years, machine learning (ML) methods have become common to solve classification problems by an automated processing in Earth and space sciences [36,37]. Machine learning methods can involve either supervised or unsupervised learning. Supervised learning techniques train machine learning algorithms using labeled data sets, which contain sample data that have been tagged with a target parameter. The algorithm uses a target to evaluate its accuracy in interpreting the training data. In unsupervised approaches, users feed unlabeled data to the algorithm, which tries to extract common features on its own based on data similarity.
The paroxysmal events do not occur often, i.e., about five events in the last 20 years [29], therefore the size of the training dataset would be insufficient to perform any supervised technique for their classification. Furthermore, some of the major events are not clearly classified by humans, thus, we designed a multivariate unsupervised classifier to automatically identify in an objective way an event as either major explosion or paroxysm. The justification behind the use of a multivariate approach lies in the fact that even individually less relevant features may become relevant when used in combination with other features and may allow a better separation among classes when used together [38], e.g., also a completely uncorrelated feature to the target can improve the separability between classes [38].
Since identifying distinct clusters through unsupervised methods depends closely on the input data, we use different independent data in order to obtain a straightforward classification that can be used any time an explosion occurs. We focus on the surface temperature time series over the volcanic summit area, the pre-and post-eruptive SAR intensity images, and satellite optical images acquired during each event and its corresponding atmospheric profile.
After collecting satellite data, we perform the feature extraction process [39], i.e., transforming raw data in informative features. In particular, we extract the magnitude of explosive volcanic deposits, a thermal index of the volcanic activity in the summit area, and the maximum plume height produced by the explosive events. Finally, the extracted features are fed to a k-means unsupervised classifier that learns to distinguish between the two classes.
The workflow of the proposed unsupervised classifier can be summarized in two main steps, as shown in Figure 2: • Feature extraction • K-means algorithm and the maximum plume height produced by the explosive events. Finally, the extracted features are fed to a k-means unsupervised classifier that learns to distinguish between the two classes. The workflow of the proposed unsupervised classifier can be summarized in two main steps, as shown in Figure

Feature Extraction
The selected features encapsulate the important aspects within data and synthesize them, providing an input for the unsupervised classifier [40]. Thus, the first step is to extract features from raw satellite data able to classify the most powerful explosive events.

Land Surface Temperature
The direct volcanic activity involved with the higher energy explosive events could lead to changes of the surface temperature over the summit area of Stromboli. These land surface temperature (LST) changes may be due to different phenomena, such as the persistent and vigorous gas emissions at the craters [41,42], or the eruption and emplacement of lava flows [22], or the multiplication and intensification of Strombolian explosive activ-

Feature Extraction
The selected features encapsulate the important aspects within data and synthesize them, providing an input for the unsupervised classifier [40]. Thus, the first step is to extract features from raw satellite data able to classify the most powerful explosive events.

Land Surface Temperature
The direct volcanic activity involved with the higher energy explosive events could lead to changes of the surface temperature over the summit area of Stromboli. These land surface temperature (LST) changes may be due to different phenomena, such as the persistent and vigorous gas emissions at the craters [41,42], or the eruption and emplacement of lava flows [22], or the multiplication and intensification of Strombolian explosive activity [43]. These phenomena are often associated with the most intense explosions detectable even before the eruption starts. The 2003, 2007, and 2019 paroxysms were heralded by intermittent phases of increases in ordinary Strombolian explosions starting weeks to months before the events and long-term anomalies in the geophysical and geochemical signals [43][44][45]. Thus, we retrieved satellite-based LST estimates to investigate surface temperature changes over the summit area of Stromboli.
The MODIS LST product (MOD11A and MYD11A) is exploited for this purpose because it provides two night-time temperature measurements per day with spatial resolution of 1 km. MODIS LST is estimated by using the generalized split-window LST algorithm [46] that utilizes two TIR bands, namely 31 and 32. This method corrects for atmospheric effects based on the differential absorption in adjacent thermal infrared bands being less sensitive to the uncertainties in optical properties of the atmosphere. A linear form for the LST algorithm is used [46]: where Ts is the effective radiometric surface temperature (i.e., land surface temperature), ε = 0.5 (ε 31 + ε 32 ) and ∆ε = ε 31 − ε 32 are the mean and the difference of surface emissivity in MODIS bands 31 and 32, T 31 and T 32 are brightness temperatures (BT) from bands 31 and 32. The coefficients C, A i and B i , i = 1, 2, 3 are given by interpolation on a set of multi-dimensional look-up tables (LUT) that can be found in [46]. We consider as an indicator of the thermal activity the LST magnitude calculated by subtracting the mean temperature from the maximum temperature measured over Stromboli Island [47]. For this purpose, we define a parameter that depends on LST in terms of both amplitude and duration. We have followed the approach proposed by Giudicepietro et al. [11] for the analysis of the Strombolian seismic signals related to volcanic dynamics. They obtained a robust parameter to represent the level of activity of Stromboli by using the amplitude and duration of the very long period seismic events. Thus, we firstly apply a smoothed 15-day moving average to the LST magnitude in order to get what we call Low Frequency Thermal Signal (LFTS), accounting for the long-term high thermal activity level before/during/after powerful explosions. The sliding window duration was chosen through a trial and error approach to obtain a good sensitivity to the LFTS signal variations. Then, we divide the entire time series in windows of 30 days and, for each time interval, we compute the max of the values; we call this parameter Land Surface Temperature Index (LSTI). This is done to investigate the persistent long-term thermal activity level; thus, considering the same time interval duration, i.e., 30 days, higher average LST magnitude values are expected for powerful explosions than ordinary sporadic explosions and consequently higher LSTI values. Therefore, when temperature changes are not observed in the selected time window, i.e., the maximum of LST is comparable with the average value over the island, it means that sporadic, weak or no activities take place; in these cases, LSTI is low. When effusive activities take place, LFTS increases from the start of the volcanic activity and stays high until new lava flows cool; this means that LST magnitude may remain high for days and thus characterized by medium/high LSTI depending on how long the effusive activity lasts. When violent explosions take place, LST magnitude may start increasing either before the eruption starts or after due to the occurrence of several concurrent phenomena raising the thermal level at the summit Crater Terrace in the selected time window, such as very high radiative power, lava flows emissions before/during/after the event, ejections of incandescent materials. In both the latter cases, the maximum LFTS values in this time window will be high implying high/extreme LSTI values. It should be clear that this parameter does not univocally identify events that are not extreme, i.e., LSTI values for small effusive events may be comparable with major explosive events. This is the reason why it is not the only parameter used to discriminate between major explosion and paroxysm but, together with other features, it can help characterize such powerful events. Therefore, we select the LSTI feature to account for the surface temperature changes over the Crater Terrace expressed in Celsius degrees.

Explosive Volcanic Deposit
The higher energy explosions of Stromboli erupt large quantities of products that fall over a widespread area far from the summit craters. Synthetic aperture radar (SAR) images have proved a valuable data source in detecting surface changes due to the emplacement of explosive volcanic deposits, including pyroclastic flows, ballistics, tephra fallout, short lava flows, as well as post-emplacement alteration [48][49][50]. Each pixel of an image corresponds to a resolution element on the ground (typically sized on the scale of the SAR wavelength ranging from several centimeters to tens of centimeters), which receives and scatters back an electromagnetic signal emitted by the radar. A pixel is characterized by two values: the amplitude and the phase. To make good use of the information content (both amplitude and phase) provided by SAR backscattering signals, detection methods are based either on the evolution of the reflectivity, namely, the amplitude of the radar images, or on variations in the correlation of the signal, which is a measure of the phase change.
The phase difference between two subsequent radar acquisitions is often used to investigate deposition of volcanic materials on the surface (e.g., [51]). However, in areas where the land surface changes rapidly, e.g., due to the emplacement of volcanic deposits, the phase difference between neighboring pixels observed by Interferometric Synthetic Aperture Radar (InSAR) will appear random (incoherent), thus little information can be retrieved. On the other hand, amplitude results from the sum of all scatters within a ground resolution element, hence slight variations in the satellite viewing geometry with time lead to the sum of a slightly different set of scatters at each time step, thus leading to pixel-by-pixel speckle pattern. Nevertheless, by comparing amplitude variations in Remote Sens. 2021, 13, 4080 8 of 21 time, changes in observed radar amplitude can be used to detect areas that are affected by volcanic eruptions, even when the SAR phase is incoherent.
Pyroclastic flows, tephra fallout and mudflow deposits produced during a volcanic activity result in changes in surface roughness [48]. Usually, deposits made by meter-sized blocks, smaller clasts and coarse tephra, appear "rough" to the radar resulting in an increase in backscattering, while fine-grained tephra appears "smooth" to the radar [48,53] resulting in a decrease in backscattering. Besides, the effects of the local morphology changes also produce changes in backscattering [48,52].
New deposits may be revealed by radar amplitude if their surfaces are sufficiently different from the underlying or surrounding deposits. Changes in surface reflectivity between two radar images can be visualized easier when plotted as a ratio image [53]. Computing ratio images allows reducing radiometric errors [54] and decreases the effect of speckle [53,55]. The image ratio depends on the relative average radar backscattering changes between two SAR images [54]. It is worth noting that for a given pixel in a SAR image, the backscatter is dominated by the largest objects contained within that pixel [56,57]. Even a small proportion of coarse-grained tephra deposits could be enough to increase the surface roughness at radar wavelengths. In these cases, the changes will appear greater in the pixels containing the largest objects, so that the maximum value of the change image will be related to the largest detected object.
We thus use Sentinel-1 SAR amplitude data to assess the magnitude of explosive volcanic deposits. As a change image, we use the amplitude ratio in order to reduce speckle noise while increasing changes in areas of high amplitude [48,54]. Therefore, we build a time series where each value at a given time t is obtained dividing the first available SAR intensity image recorded after the time t by the SAR intensity image recorded at the time t. In particular, in order to avoid misleading results due to changes not related to volcanic activities, e.g., vegetation changes occurring on the island and changes in the dielectric constant due to the presence of the sea, we investigated only the area of the Sciara del Fuoco. When a powerful explosion occurs, the amplitude ratio that is computed by using post-eruptive and pre-eruptive images will contain pixels with higher values since larger and coarser deposits are supposed to be emplaced. On the other hand, values smaller than one mean that erosion occurred and materials were removed. The latter information could thus be used to detect the material that is eroded and collapsed during an intense event, such as the inner crater removal during a paroxysm [6,8].
Sentinel-1 SAR amplitude data are gathered from the GEE collection "COPERNI-CUS/S1_GRD_FLOAT" that includes the Ground Range Detected (GRD) scenes, already processed for thermal noise removal, radiometric calibration and terrain correction. We used a local adaptive median filter that uses local statistics to detect SAR speckle noise and to replace it with a local median value. In comparison with established filters, it was found that the local adaptive median filter outperformed others in achieving the best balance between speckle suppression and image detail preservation [58]. Single polarization data ('VV') is used [48,49] with descending pass since descending-pass image changes are stronger [48]. This feature based on the Sentinel-1 SAR measurement is named Explosive Volcanic Deposit (EVD) from now on. It is an adimensional parameter being the maximum of the ratio of two SAR intensity images.

Plume Height
The height of the eruptive column (i.e., the height of the ash plume rising by buoyancy from the crater rim) is related to the magnitude of the explosion and thus to the volume of erupted products [59], and as such is a key discriminant between events of different magnitude and intensity. The plume produced during an explosion reaches greater heights during powerful events such as paroxysms [28]. Thus, we derived the volcanic cloud top altitude by means of satellite infrared images. For this purpose, we use the brightness temperature of the cloud surface at 10.8 µm [60,61]. Assuming a thermal equilibrium between the plume top height and the atmosphere, we acquire the brightness temperature T p in the Thermal InfraRed (TIR) wavelength of the plume pixels and compare it with the atmospheric temperature-altitude profile on Stromboli at the same time of satellite acquisition. We use the atmospheric temperature-altitude profile obtained from the National Centers for Environmental Prediction (NCEP)/National Center for Atmospheric Research (NCAR) Reanalysis Project (available at http://www.esrl.noaa.gov/psd/data/reanalysis/reanalysis.shtml, accessed on 20 June 2020). This technique is known as the Cloud Top Temperature [62] and relies on several assumptions: (i) the cloud must be optically thick (thickness τ > 5; where τ is a dimensionless quantity that indicates how transparent is a cloud) so that the brightness temperature is representative of the true cloud surface temperature; (ii) the ash particles must be in thermal equilibrium with the temperature of the ambient air; (iii) this method is only applicable to tropospheric emissions since the atmospheric temperature profile above the tropopause is not monotonic [63].
All the available SEVIRI satellite TIR images were taken into account. Although SEVIRI is available every 15 min, it has a low spatial resolution, thus in order to obtain a good level of detail in terms of spatial resolution, when available, the VIIRS data acquired at 375 m spatial resolution were used too. A linear interpolation was performed by locating T p within the temperature profile. Equation (2) shows how the maximum plume height PH was estimated using the interpolation technique [64]: where T 1 and T 2 are the temperatures within the profile that bound T p , with T 1 being the temperature at the highest (e.g., furthest from the ground) bounding level, while Z 1 and Z 2 are the corresponding heights of the bounding temperatures, T 1 and T 2 .

K-Means Clustering
Extracted feature vectors can become inputs for a variety of different techniques of machine learning. We use cluster analysis (CA) to classify the most powerful explosions of Stromboli. Cluster analysis is an unsupervised learning technique aimed at grouping similar data while separating different ones, where similarity is evaluated quantitatively using a distance function in the space of feature vectors. The clustering algorithms can be distinguished into hierarchical and non-hierarchical. In the former, a tree-like structure is built to represent the relations between clusters, while in the latter, new clusters are formed by merging or splitting existing ones without following a tree-like structure but simply grouping the data in order to maximize or minimize some evaluation criteria.
CA includes a vast class of algorithms. We have chosen the non-hierarchical k-means algorithm because it is simple and efficient [65][66][67]. It is an iterative procedure that partitions data into a predetermined number of clusters [68]. Firstly, inside the feature space, cluster centroids are randomly initialized and observations are assigned to the nearest centroid according to a Euclidean distance. Secondly, cluster centroids are updated with the mean location within their cluster. These steps are repeated until an optimal convergence is reached. We have used the k-means++ algorithm for cluster center initialization using a heuristic to find centroid seeds for k-means clustering. According to Arthur and Vassilvitskii [69], k-means++ improves the running time of Lloyd's algorithm [70] and the quality of the final solution. k-means++ algorithm achieves faster convergence to a lower sum of within-cluster, sum-of-squares point-to-cluster-centroid distances, than Lloyd's algorithm.
A challenge in applying the k-means algorithm is determining the number (k) of clusters to partition our dataset. When available, a priori knowledge is used to inform the choice of k. Our goal being to classify major explosions and paroxysms, we set k = 2 number of clusters, namely "major explosion" and "paroxysm" clusters.

Results
Once all the features are retrieved, a z-score normalization is performed to scale the data and reduce variability across arrays [71]. Then, an input vector made by the concatenation of the three normalized parameters PH, EVD, and LSTI is built for each violent explosion (listed in Table 1) occurring at Stromboli between January 2018 and April 2021, x = (PH, EVD, LSTI). The input vectors are then concatenated together getting the input matrix X, an n × m matrix where n is the number of event samples, i.e., n = 19 in this case, and m is the number of distinguishing features, i.e., m = 3 namely PH, EVD, and LSTI. Optical and Radar satellite images used to obtain EVD and PH values for each investigated volcanic explosion are listed in Table 2. The LST derived features for all the monitored events are shown in Figure 3, where both the LFTS (Figure 3a) and LSTI (Figure 3b) are reported. It is worth noting that these estimates are affected by any volcanic phenomena that leads to persistent increase in temperature over the summit area, particularly lava flow emissions. Thus, mid/high values of LSTI are found for time intervals when emission of lava occur, such as January-February 2020 (see INGV weekly bulletin).
The magnitude of explosive deposits retrieved for all the investigated events are shown in Figure 3c. Each value at a given time t is the EVD value derived dividing the first available SAR intensity image recorded after t by the SAR intensity image recorded at the time t. The explosive deposits associated with the investigated eruptive events listed in Table 1 are shown as red dots.
The maximum heights PH reached by the eruptive columns, derived from SEVIRI and VIIRS images, are reported in Table 3. In the case a satellite-based estimate cannot be retrieved due to sensor limits, e.g., in presence of low and transparent volcanic clouds, the values from the INGV weekly bulletins on the monitoring activity were considered. Note that the events with missing information on the plume due to the detectability limit of SEVIRI, are confirmed to be characterized by lower plume height compared with the estimated value of 3 July 2019 event (INGV weekly bulletins). Below, we focus on one event of each class, namely 3 July 2019 paroxysm and 10 November 2020 major explosion. The descriptions of the temporal evolution of the two events were derived from INGV bulletins and from Calvari et al. [7].

The 3 July 2019 Paroxysm
The 3 July 2019 event, one of the most violent explosions occurring at Stromboli volcano, was preceded by a significant increase in the intensity of explosive activity at all the summit vents. At 13:59 a strong explosion from the South West Crater (SWC) was followed by lava flow output along the upper part of the Sciara del Fuoco. A small vent opened at the base of North East Crater (NEC), widening and feeding a lava flow that started spreading along the Sciara del Fuoco from 14:02:40. Several lava flows from this vent overlapped the previous flows, spreading downslope. In the meantime, at 14:43:10 another small lava flow started from the Central Crater (CC), slowly spreading south within the Crater Terrace, followed by flows from the two NEC vents, and eventually joining with another flow erupted from the SWC vent. The maximum velocity of the jet,  Table 1. (c) EVD is computed for each Sentinel-1 acquisition (blue dots). EVD values corresponding to the investigated eruptive events are shown as red dots. Table 3. List of the input features and classification outcome between "paroxysm" and "major explosion" for each eruptive event at Stromboli volcano. Date and time in dd/mm/yyyy hh:mm:ss format. Time is expressed in UTC.  Below, we focus on one event of each class, namely 3 July 2019 paroxysm and 10 November 2020 major explosion. The descriptions of the temporal evolution of the two events were derived from INGV bulletins and from Calvari et al. [7].

The 3 July 2019 Paroxysm
The 3 July 2019 event, one of the most violent explosions occurring at Stromboli volcano, was preceded by a significant increase in the intensity of explosive activity at all the summit vents. At 13:59 a strong explosion from the South West Crater (SWC) was followed by lava flow output along the upper part of the Sciara del Fuoco. A small vent opened at the base of North East Crater (NEC), widening and feeding a lava flow that started spreading along the Sciara del Fuoco from 14:02:40. Several lava flows from this vent overlapped the previous flows, spreading downslope. In the meantime, at 14:43:10 another small lava flow started from the Central Crater (CC), slowly spreading south within the Crater Terrace, followed by flows from the two NEC vents, and eventually joining with another flow erupted from the SWC vent. The maximum velocity of the jet, estimated at 101.92 m s −1 , was recorded after 6 s of gradual increase due to the initial horizontal expansion of the jet [7]. The main phase of the paroxysm involved the whole crater area. The column collapse produced two pyroclastic flows that spread north-west along the Sciara del Fuoco and over the sea surface for several hundred meters, and caused a significant widening of the Crater Terrace [7]. In Figure 4a, a zoomed graph of the LFST is reported to highlight the increase begun before the 3 July 2019 paroxysm, for which LSTI = 2.63 • C.
ns. 2021, 13, x FOR PEER REVIEW 13 of 21 estimated at 101.92 m s −1 , was recorded after 6 s of gradual increase due to the initial horizontal expansion of the jet [7]. The main phase of the paroxysm involved the whole crater area. The column collapse produced two pyroclastic flows that spread north-west along the Sciara del Fuoco and over the sea surface for several hundred meters, and caused a significant widening of the Crater Terrace [7]. In Figure 4a, a zoomed graph of the LFST is reported to highlight the increase begun before the 3 July 2019 paroxysm, for which LSTI = 2.63 °C. The explosive deposits were assessed by using the amplitude ratio of two Sentinel-1 SAR descending intensity images. We used the 26 June 2019 and 08 July 2019 images as the pre-eruptive and post-eruptive respectively. The pixel by pixel ratio of the recorded images of the two dates is shown in Figure 4b. The red rectangle indicates the presence of The explosive deposits were assessed by using the amplitude ratio of two Sentinel-1 SAR descending intensity images. We used the 26 June 2019 and 08 July 2019 images as the pre-eruptive and post-eruptive respectively. The pixel by pixel ratio of the recorded images of the two dates is shown in Figure 4b. The red rectangle indicates the presence of new materials deposited after the eruption, while the green rectangle indicates the erosion of material after the eruption with respect to the pre-eruptive image. The gray color scale indicates the intensity of the backscattered signal from lower (black) to higher (white). The maximum value detected by Sentinel-1 SAR is EVD = 12.5.
The plume height has been retrieved by using the BT channel at 10.8 um of SEVIRI image acquired on 3 July 2019 at 15:00 UTC resulting in T p = 264 K, while from the atmospheric profiles we obtain Z 1 = 4450 m, T 1 = 274 m, Z 2 = 5888 m, T 2 = 263.8 K. The top plume height PH is estimated at 5.86 km in accordance with Equation (2) (Figure 5b).
July 2019 paroxysmal event. (b) Change detection image showing explosive deposits and eroded area detected by using the ratio of Sentinel-1 SAR descending images. The amplitude ratio is calculated between the pre-eruptive image of 26 June 2019 and the post-eruptive intensity image of 8 July 2019. The red rectangle indicates the brighter pixels that are related to the biggest materials deposited (widespread pyroclastic density current, ballistics and tephra fallout) during the eruption, whereas the green rectangle indicates the darker pixels related to the removal of materials due to the crater area erosion after the eruption.
The explosive deposits were assessed by using the amplitude ratio of two Sentinel-1 SAR descending intensity images. We used the 26 June 2019 and 08 July 2019 images as the pre-eruptive and post-eruptive respectively. The pixel by pixel ratio of the recorded images of the two dates is shown in Figure 4b. The red rectangle indicates the presence of new materials deposited after the eruption, while the green rectangle indicates the erosion of material after the eruption with respect to the pre-eruptive image. The gray color scale indicates the intensity of the backscattered signal from lower (black) to higher (white). The maximum value detected by Sentinel-1 SAR is EVD = 12.5.
The plume height has been retrieved by using the BT channel at 10.8 um of SEVIRI image acquired on 3 July 2019 at 15:00 UTC resulting in Tp = 264 K, while from the atmospheric profiles we obtain Z1 = 4450 m, T1 = 274 m, Z2 = 5888 m, T2 = 263.8 K. The top plume height PH is estimated at 5.86 km in accordance with Equation 2 (Figure 5b).

The 10 November 2020 Major Explosion
The 10 November 2020 episode started from the SWC at 20:04:21, forming an eruptive cloud reaching up to 600 m above the crater rim [7]. The explosive event then expanded to the NEC, forming a jet propagating horizontally and causing a wide spatter fallout on the upper Sciara del Fuoco, and eventually expanding to the Central Crater producing a low fountaining with little or no fallout outside the crater. The duration of the first pulse was 20 s, and the muzzle velocity of the ejecta was 54.50 m s −1 [7]. The index related to the temperature, as shown in Figure 3, is LSTI = 1.82 • C and a zoom of the LFTS is shown in Figure 6a.
The explosive deposits were assessed by using the ratio between the pre-eruptive Sentinel-1 SAR descending intensity image of 5 November 2020 and the post-eruptive intensity image of 11 November 2020. The pixel by pixel image ratio is shown in Figure 6b. The gray color scale indicates the intensity of the backscattered signal from lower (black) to higher (white). The maximum value detected by Sentinel-1 SAR is EVD = 2.28. The plume height was retrieved by using the INGV bulletin, as the plume was not detectable by using SEVIRI images, thus PH = 600 m.
Sentinel-1 SAR descending intensity image of 5 November 2020 and the post-eruptive intensity image of 11 November 2020. The pixel by pixel image ratio is shown in Figure 6b. The gray color scale indicates the intensity of the backscattered signal from lower (black) to higher (white). The maximum value detected by Sentinel-1 SAR is EVD = 2.28. The plume height was retrieved by using the INGV bulletin, as the plume was not detectable by using SEVIRI images, thus PH = 600 m.

Explosive Events Classification
Using the values of three distinguishing features PH, EVD and LSTI, the input matrix X was populated for all 19 event samples (Table 3). Then, the input matrix X of dimension [19 × 3] was given to the k-means classifier that grouped the investigated events as is shown in Figure 7. Three events are clustered together in one class, whereas the remaining sixteen events are grouped in the other one. The cluster identified by the centroid with higher coordinates is representative of the more intensive explosions and is named "Paroxysm", the other one characterized by the centroid with lower coordinates is named "Major Explosion".

Explosive Events Classification
Using the values of three distinguishing features PH, EVD and LSTI, the input matrix X was populated for all 19 event samples (Table 3). Then, the input matrix X of dimension [19 × 3] was given to the k-means classifier that grouped the investigated events as is shown in Figure 7. Three events are clustered together in one class, whereas the remaining sixteen events are grouped in the other one. The cluster identified by the centroid with higher coordinates is representative of the more intensive explosions and is named "Paroxysm", the other one characterized by the centroid with lower coordinates is named "Major Explosion".  Table 1.

Discussion
Unsupervised learning captures the eruptive patterns characterizing the most pow erful explosions at Stromboli. By combining the information provided from independen features for each event, a unique assignment to a specific representative class is automat

Discussion
Unsupervised learning captures the eruptive patterns characterizing the most powerful explosions at Stromboli. By combining the information provided from independent features for each event, a unique assignment to a specific representative class is automatically performed based on satellite data.
On one hand, results show that the major explosions are characterized by low/medium thermal content, i.e., LSTI around 1.4 • C, low plume height, i.e., PH around 420 m, and low production of explosive deposits, i.e., EVD around 2.5. However, the two events occurring on 13 and 15 July 2019 are distant (named eccentric) from the other major explosions in LSTI dimension. This happens because these events occur at a time very close to the 3 July 2019 paroxysm and all share similar LSTI values. During this paroxysmal phase, lava is erupted for days thus both LFTS and LSTI have high values for the entire period. Nevertheless, these events are classified as major explosions since the other features, i.e., PM and EVD, have characteristic values for this class.
On the other hand, results show that paroxysms are extreme events characterized mainly by medium/high thermal content, i.e., LSTI around 2.3 • C, medium/high plume height, i.e., PH around 3330 m, and high production of explosive deposits, i.e., EVD around 10.17. Some events may differ in the way they show up as happened for the powerful explosion occurring on 16 November 2020. It has weaker LSTI and PH than the 3 July and 28 August 2019 paroxysms, but very high values of EVD, leading the k-means algorithm to classify the 16 November 2020 event as a paroxysm (Figure 8). It is noteworthy that this event was preceded by significant ground deformation with the upward tilting of the northeast outer flank of the cone and characterized by the failure of the north portion of the Crater Terrace and by a lateral blast [7], thus EVD has a much larger value than those of major explosions.

FOR PEER REVIEW 16 of 21
On one hand, results show that the major explosions are characterized by low/medium thermal content, i.e., LSTI around 1.4 °C, low plume height, i.e., PH around 420 m, and low production of explosive deposits, i.e., EVD around 2.5. However, the two events occurring on 13 and 15 July 2019 are distant (named eccentric) from the other major explosions in LSTI dimension. This happens because these events occur at a time very close to the 3 July 2019 paroxysm and all share similar LSTI values. During this paroxysmal phase, lava is erupted for days thus both LFTS and LSTI have high values for the entire period. Nevertheless, these events are classified as major explosions since the other features, i.e., PM and EVD, have characteristic values for this class.
On the other hand, results show that paroxysms are extreme events characterized mainly by medium/high thermal content, i.e., LSTI around 2.3 °C, medium/high plume height, i.e., PH around 3330 m, and high production of explosive deposits, i.e., EVD around 10.17. Some events may differ in the way they show up as happened for the powerful explosion occurring on 16 November 2020. It has weaker LSTI and PH than the 3 July and 28 August 2019 paroxysms, but very high values of EVD, leading the k-means algorithm to classify the 16 November 2020 event as a paroxysm (Figure 8). It is noteworthy that this event was preceded by significant ground deformation with the upward tilting of the northeast outer flank of the cone and characterized by the failure of the north portion of the Crater Terrace and by a lateral blast [7], thus EVD has a much larger value than those of major explosions. The variable phenomenology of the most powerful explosions highlights the importance of using as many parameters as possible to characterize an event well rather than using a single measurement, enabling the classifier to provide an objective measure of the monitored explosion. It is worth noting that our results are in general agreement with the typologies shown in Table 1. However, the 19 July 2020 volcanic event that was not uniquely classified [7;29], is classified here as "Major Explosion"; whereas the 16 November 2020 event that was uniquely classified as "Major Explosion", is classified here as "Paroxysm". This observation is coherent with the fact that such event was violent and powerful ( [7]; INGV weekly bulletin) with the plume height reaching a higher value compared with major explosions, and which generated/caused a pyroclastic density current that reached the sea and a landslide that was clearly detected by SAR, i.e., high EVD value (see Table 3). Besides powerful explosive events, low/medium LSTI and medium/high EVD values may also be found in case of lava flow emissions as is shown in Figure 3, e.g., February-March 2020 events. That is why these features cannot be considered individually as the only distinguishing factors characterizing an explosive event.
Among the identified paroxysms, the most extreme one occurred on 3 July 2019, as indicated by all the investigated features. The LFTS appears to increase weeks before the paroxysm started (Figure 4a). Pre-eruptive long-term (years) variations in MODIS Band The variable phenomenology of the most powerful explosions highlights the importance of using as many parameters as possible to characterize an event well rather than using a single measurement, enabling the classifier to provide an objective measure of the monitored explosion. It is worth noting that our results are in general agreement with the typologies shown in Table 1. However, the 19 July 2020 volcanic event that was not uniquely classified [7;29], is classified here as "Major Explosion"; whereas the 16 November 2020 event that was uniquely classified as "Major Explosion", is classified here as "Paroxysm". This observation is coherent with the fact that such event was violent and powerful ( [7]; INGV weekly bulletin) with the plume height reaching a higher value compared with major explosions, and which generated/caused a pyroclastic density current that reached the sea and a landslide that was clearly detected by SAR, i.e., high EVD value (see Table 3). Besides powerful explosive events, low/medium LSTI and medium/high EVD values may also be found in case of lava flow emissions as is shown in Figure 3, e.g., February-March 2020 events. That is why these features cannot be considered individually as the only distinguishing factors characterizing an explosive event.
Among the identified paroxysms, the most extreme one occurred on 3 July 2019, as indicated by all the investigated features. The LFTS appears to increase weeks before the paroxysm started (Figure 4a). Pre-eruptive long-term (years) variations in MODIS Band 31 median radiant temperature have already been observed in different volcanoes before magmatic and phreatic eruptions [72]. In these cases, the large-scale (tens of square kilometers) changes of surface temperature provided an early indicator of renewed volcanic activity. Similarly, the significant long-term (from weeks to months) LFTS variations preceding the 3 July 2019 paroxysm can be an early detection of pre-eruptive conditions driving the most violent explosive eruptions.
This behavior is coherent with the changes observed in the parameters related to the dynamics of very long period (VLP) seismic events associated with Strombolian explosive activity, which were interpreted as due to the temporal evolution of the explosive source, indicating higher gas emissions at the summit craters, starting at least one month before the 3 July 2019 paroxysm [11,73]. This analysis pointed to a critical role of the gas as well as of the coexistence of two magma-types (with different gas content and porphyricity) for the eruptive dynamics of Stromboli [73].
It should be pointed out that LFTS remains high between the 3 July 2019 and 28 August 2019 events as a result of the on-going effusive eruption, and that the 28 August 2019 explosive event ended the July-August 2019 effusive phase [73]. Then the LFTS signal returned to lower values. This sudden drop is consistent with the interpretation that after the 3 July 2019 eruption the transitory conditions that enabled pressurization in the shallow conduit were not stabilized until the 28 August 2019 paroxysm, when a complete restoration of the gas flux recharge/discharge ratio in the volcanic conduit occurred together with a return to the ordinary volcanic activity [73]. This is also coherent with the high thermal content characterizing the events on 13 July 2019 and 15 July 2019. On the other hand, no evident long-term gradual changes in LFTS have been observed in other events, e.g., 10 November 2020 major explosion (Figure 6a).
Different values of the plume height of the 3 July 2019 event were estimated in the previous studies and reports. The maximum plume height was computed ranging: from 4 km [7], to 7 km [60], and about 8.4 km [6]. Here, we found a value of PH = 5.86 km, very close to the estimate of 6 km obtained with numerical simulations [43].
Besides the high detected EVD estimate, the change detection image in Figure 4b also shows the areas that were eroded during the 3 July 2019 paroxysm. In particular, the portion of the Crater Terrace that has been eroded appears as a dark area (green rectangle in Figure 4b). Thus, the bright areas include both the pyroclastic deposit emitted during the paroxysmal eruption and part of the Crater Terrace eroded that fell down on the slope during the explosion. The distribution of the deposits was well-documented through field surveys [43] and is in good agreement with the deposits detected in Figure 4b. The 3 July 2019 paroxysm deposited a 0.3-1 m (0.4 m average) thick continuous deposit on the Pizzo area upon which a discontinuous layer of 0.5-3 m-sized spatter bombs was laid. Tephra fallout deposits ranging from ash-sized to bomb-sized pyroclasts were dispersed towards south west and a continuous, coarse-grained (10-50 cm) pumice bombs deposit, up to 0.7 m thick in the proximal sector, covered the upper west flank down to 500 m a.s.l. [43].
Lastly, since the 3 July 2019 paroxysm was an exceptional and powerful event, fires broke out at some distance from the craters [5]. Thus, we have processed Sentinel-2 Multi-Spectral Instrument (MSI) images again to detect thermal anomalies. We have used the data collected on 4 July 2019 at 09:49 UTC to acquire the False Color image made by B12 (SWIR)-B11 (SWIR)-B4 (RED) bands shown in Figure 9. Although the satellite pass does not cover the whole of Stromboli Island, fires in the southern and eastern areas are clearly visible, thus confirming the violence of this paroxysm.
Lastly, since the 3 July 2019 paroxysm was an exceptional and powerful event, fires broke out at some distance from the craters [5]. Thus, we have processed Sentinel-2 Multi-Spectral Instrument (MSI) images again to detect thermal anomalies. We have used the data collected on 4 July 2019 at 09:49 UTC to acquire the False Color image made by B12 (SWIR)-B11 (SWIR)-B4 (RED) bands shown in Figure 9. Although the satellite pass does not cover the whole of Stromboli Island, fires in the southern and eastern areas are clearly visible, thus confirming the violence of this paroxysm.

Conclusions
The different sizes, covering a broad variability in intensity and magnitude, of the most intense explosive eruptions at Stromboli volcano from January 2018 to April 2021 were classified by combining satellite imagery and machine learning techniques. Modern satellites provide wide coverage with high resolution signals, generating large datasets. These remote sensing observations allowed us to obtain information on the monitored phenomena not retrievable by using only ground-based sensors, i.e., thermal cameras of the monitoring network. In particular, spatial mapping of the explosive volcanic deposits ejected during such events is derived from SAR images. Furthermore, the Sentinel-1 SAR amplitude signal has proved to be a valuable means not only to detect explosive volcanic deposits ejected during intense explosive events, but also to gain insights into the collapsed portion of the Crater Terrace.
We highlighted the importance of the information borne by the LST signal from two viewpoints. On one hand, it reveals changes due to the explosive source leading to intense paroxysmal explosions weeks before the event begins. This can be a first indicator of subsurface processes driving the most powerful explosions that can be tracked through satellite-based infrared observations with high temporal resolution. These findings give new insights into the importance of surface temperature variations to anticipate the most powerful eruptions that are not easy to forecast through traditional geophysical and geochemical measurements. On the other hand, the LST magnitude level can be used to characterize the most powerful events together with other satellite-derived features.
We demonstrated the capability of machine learning algorithms in distinguishing the most powerful explosions at Stromboli using a variety of satellite data. An unsupervised k-means classifier is proposed to automatically recognize explosive events as either paroxysms or major explosions based on the information related to the surface temperature of the summit area, the height reached by the plume and the magnitude of explosive deposits including pyroclastic flows, ballistics, and tephra fallout. Centroids are defined for both clusters, thus providing the main values of features that are characteristic for each class of explosive eruption. These results indicate that machine learning algorithms combined with satellite data have the potential to form an automated system to detect the intensity of explosive eruptions, also in remote and inaccessible volcanoes.

Data Availability Statement:
The European Space Agency and NASA provided high-quality openaccess data to the scientific community and Google Earth Engine platform facilitated the access to and the processing of the archive of publicly available satellite imagery.