Volcanic Hot-Spot Detection Using SENTINEL-2: A Comparison with MODIS–MIROVA Thermal Data Series

: In the satellite thermal remote sensing, the new generation of sensors with high-spatial resolution SWIR data open the door to an improved constraining of thermal phenomena related to volcanic processes, with strong implications for monitoring applications. In this paper, we describe a new hot-spot detection algorithm developed for SENTINEL-2 / MSI data that combines spectral indices on the SWIR bands 8a-11-12 (with a 20-meter resolution) with a spatial and statistical analysis on clusters of alerted pixels. The algorithm is able to detect hot-spot-contaminated pixels (S2Pix) in a wide range of environments and for several types of volcanic activities, showing high accuracy performances of about 1% and 94% in averaged omission and commission rates, respectively, underlining a strong reliability on a global scale. The S2-derived thermal trends, retrieved at eight key-case volcanoes, are then compared with the Volcanic Radiative Power (VRP) derived from MODIS (Moderate Resolution Imaging Spectroradiometer) and processed by the MIROVA (Middle InfraRed Observation of Volcanic Activity) system during an almost four-year-long period, January 2016 to October 2019. The presented data indicate an overall excellent correlation between the two thermal signals, enhancing the higher sensitivity of SENTINEL-2 to detect subtle, low-temperature thermal signals. Moreover, for each case we explore the speciﬁc relationship between S2Pix and VRP showing how di ﬀ erent volcanic processes (i.e., lava ﬂows, domes, lakes and open-vent activity) produce a distinct pattern in terms of size and intensity of the thermal anomaly. These promising results indicate how the algorithm here presented could be applicable for volcanic monitoring purposes and integrated into operational systems. Moreover, the combination of high-resolution (S2 / MSI) and moderate-resolution (MODIS) thermal timeseries constitutes a breakthrough for future multi-sensor hot-spot detection systems, with increased monitoring capabilities that are useful for communities which interact with active volcanoes. of lava lava-lake pulsations, fumarolic activity, multiple active craters thermal

and location of volcanic hot thermal anomalies with high accuracy, in order to provide useful and reliable information about thermal signals related to ongoing volcanic activity. In this regard, we preferred to avoid calculation of pixel-integrated temperature or radiant fluxes, because the SWIR bands at a high spatial resolution (such as ASTER and OLI/LANDSAT-8 sensors) offer smaller pixels and a more sensitivity to detect hot magmatic features (200-1200 • C), and thus are often saturated over high-temperature targets (fires, volcanoes), since a hot surface constitutes a greater part of pixel footprint [21,48]. In fact, some previous works on fire detection, using SWIR wavelengths, reported common saturation issues, leading to counterfeit DN pixel values and artificial radiances conversions, thereby hindering accurate quantitative analyses [43][44][45]. This is one of the first works, even if embryonic, to specifically design a volcano thermaldetection algorithm. using S2 high-resolution imagery, tested on several different volcanic cases through a multi-year-long analysis. Moreover, it is the first attempt to compare the two different thermal signals derived from SENTINEL-2 MSI and MODIS sensors on volcanic environments. The algorithm was developed in order to be applicable for volcanic-monitoring purposes at a global scale and to be integrated into thermal space-based systems, such as MIROVA, and in multiparametric applications. In this respect, the algorithm here presented works operationally nowadays on MOUNTS (Monitoring Unrest from Space, http://www.mounts-project.com/home; [30]), a multiparametric monitoring satellite system combining SAR, UV and IR analysis, using SENTINEL constellation.
The final goal of this work is to demonstrate the validity of the here-proposed detection algorithm for the SENTINEL-2 high-resolution data (and possibly exportable on LANDSAT-8 imagery), with the aim to provide to the volcanological community-from single researcher to volcanic observatories and monitoring centers-a useful and solid tool to evaluate the presence and size of active hot spot, to detect possible volcanic precursor thermal signals, and to track the thermal and spatial evolution of active volcanic portions, such as single vents, fumaroles fields or lava bodies.

SENTINEL-2 Products
The SENTINEL-2 mission consists of the two platforms, 2A and 2B, launched in June 2015 and March 2017, respectively, by the European Space Agency (Copernicus program). They are located at 786 km altitude and placed in the same sun-synchronous polar orbit at 180° from each other [49]. The Multispectral Instrument (MSI) onboard captures multispectral images at 13 different bands, spanning from visible (VIS) to short-wave infrared (SWIR) wavelengths, with spatial resolution of This is one of the first works, even if embryonic, to specifically design a volcano thermal-detection algorithm. using S2 high-resolution imagery, tested on several different volcanic cases through a multi-year-long analysis. Moreover, it is the first attempt to compare the two different thermal signals derived from SENTINEL-2 MSI and MODIS sensors on volcanic environments. The algorithm was developed in order to be applicable for volcanic-monitoring purposes at a global scale and to be integrated into thermal space-based systems, such as MIROVA, and in multiparametric applications. In this respect, the algorithm here presented works operationally nowadays on MOUNTS (Monitoring Unrest from Space, http://www.mounts-project.com/home; [30]), a multiparametric monitoring satellite system combining SAR, UV and IR analysis, using SENTINEL constellation.
The final goal of this work is to demonstrate the validity of the here-proposed detection algorithm for the SENTINEL-2 high-resolution data (and possibly exportable on LANDSAT-8 imagery), with the aim to provide to the volcanological community-from single researcher to volcanic observatories and monitoring centers-a useful and solid tool to evaluate the presence and size of active hot spot, to detect possible volcanic precursor thermal signals, and to track the thermal and spatial evolution of active volcanic portions, such as single vents, fumaroles fields or lava bodies.

SENTINEL-2 Products
The SENTINEL-2 mission consists of the two platforms, 2A and 2B, launched in June 2015 and March 2017, respectively, by the European Space Agency (Copernicus program). They are located Remote Sens. 2020, 12, 820 5 of 32 at 786 km altitude and placed in the same sun-synchronous polar orbit at 180 • from each other [49]. The Multispectral Instrument (MSI) onboard captures multispectral images at 13 different bands, spanning from visible (VIS) to short-wave infrared (SWIR) wavelengths, with spatial resolution of 10, 20 and 60 m ( Table 1). The revisit time, considering both orbiting platforms, spans from 5 days at equator latitudes to 2-3 days at midlatitudes [50]. The SENTINEL-2 constellation provides coherence with LANDSAT-type image data, contributing to ongoing multispectral observation [50,51]. Here, we used the SENTINEL-2 Level-1C products, composed of 100 × 100 km 2 tiles (or granules) orthorectified in UTM/WGS84 projection, where the radiometric data are furnished in Top of Atmosphere (TOA) reflectances (ρ), together with the conversion parameters for radiances calculation.

Data Access
SENTINEL-2 data are originally distributed by the ESA Copernicus Open Access Hub [52], and in the past years, other visualization and storage platforms came out, which can be of great help to whoever takes the first steps to manage satellite data. Amongst these, the cloud storage service of Amazon Web Service S3 (AWS-S3, [53]) hosts SENTINEL-2 data, which are added regularly, usually within few hours after they are available on Copernicus OpenHub. One of the advantages of this service (AWS-S3) consists in the possibility to download just single bands of interest, instead of the entire SENTINEL-2 L1C product, greatly reducing the amount of data to be downloaded for specific applications. By considering three bands with 20 m of resolution, as used by the here-proposed algorithm, this opportunity allows us to reduce the data volume from~600 MB (full SENTINEL-2 Level-1C size) to that of the three downloaded bands (i.e., less than~100 MB). Conversely, it is important to point out that, since summer 2018, this service has imposed a pay request for each download ($0.05-$0.09 per GB, [54]; roughly, less than $1 for a one-year thermal monitoring, using the three SWIR bands per volcano). After downloading the bands of interest in JPEG 2000 format and the supporting metadata files (Granule_Metadata_File and Product_Level_Metadata containing respectively the geometric/satellite and sensing data), granules products are cropped in a mask of 10 km size (501 × 501 pixel), centered on the coordinates of each volcano summit of interest. For consistency with the MIROVA system (see [25]), we used the Global Volcanism Program volcano database [55] as a reference for both volcano names and coordinates.

Hot-Spot Algorithm
According to several authors [43][44][45][46], for the hot-spot detection during daytime, SWIR reflectances are preferred to radiances because (i) they are normalized to the incoming solar radiation, reducing the effects of Earth-sun-satellite geometry; and (ii) they allow for better compatibility between different satellite sensors (SENTINEL-2 MSI and LANDSAT-8 OLI, mainly), making the algorithm easily exportable to other sensors [43]. The presented algorithm is based on the analysis on the TOA Remote Sens. 2020, 12, 820 6 of 32 reflectances of three bands in the SWIR regions, ρ8a (865 nm), ρ11 (1610 nm) and ρ12 (2190 nm), all having a 20 m pixel size resolution ( Figure 2). These bands are chosen because the nearest to higher wavelengths and falling in the SWIR region, extremely sensitive to radiations from very hot surfaces (e.g., fires and fresh lava; [21]). In addition, the selected S2 bands are the equivalent of bands 7-6-5 OLI LANDSAT-8 sensor, previously used for fires detection [43][44][45]. The algorithm is essentially divided into three main steps focused in detecting all the hot-spot-contaminated pixels present in the analyzed image ( Figure 3). In the first step, a series of logical test is used to identify pixels potentially contaminated by hot spots. The second step calculates a Thermal Index (hereby considered a proxy of the pixel integrated temperature) for each alerted pixel and group the adjacent pixels into distinct clusters. In the third step, each cluster is analyzed spatially and statistically in order to reduce effects of thermal halo, diffractions spikes or "blurring" (adjacency effects, i.e., the scattering of light reflected from the nearby land into the sensor's field of view, due to contrast between a target pixel and its neighborhood; [56]).
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 33 algorithm easily exportable to other sensors [43]. The presented algorithm is based on the analysis on the TOA reflectances of three bands in the SWIR regions, ρ8a (865 nm), ρ11 (1610 nm) and ρ12 (2190 nm), all having a 20 m pixel size resolution ( Figure 2). These bands are chosen because the nearest to higher wavelengths and falling in the SWIR region, extremely sensitive to radiations from very hot surfaces (e.g., fires and fresh lava; [21]). In addition, the selected S2 bands are the equivalent of bands 7-6-5 OLI LANDSAT-8 sensor, previously used for fires detection [43][44][45]. The algorithm is essentially divided into three main steps focused in detecting all the hot-spot-contaminated pixels present in the analyzed image ( Figure 3). In the first step, a series of logical test is used to identify pixels potentially contaminated by hot spots. The second step calculates a Thermal Index (hereby considered a proxy of the pixel integrated temperature) for each alerted pixel and group the adjacent pixels into distinct clusters. In the third step, each cluster is analyzed spatially and statistically in order to reduce effects of thermal halo, diffractions spikes or "blurring" (adjacency effects, i.e., the scattering of light reflected from the nearby land into the sensor's field of view, due to contrast between a target pixel and its neighborhood; [56]).    algorithm easily exportable to other sensors [43]. The presented algorithm is based on the analysis on the TOA reflectances of three bands in the SWIR regions, ρ8a (865 nm), ρ11 (1610 nm) and ρ12 (2190 nm), all having a 20 m pixel size resolution ( Figure 2). These bands are chosen because the nearest to higher wavelengths and falling in the SWIR region, extremely sensitive to radiations from very hot surfaces (e.g., fires and fresh lava; [21]). In addition, the selected S2 bands are the equivalent of bands 7-6-5 OLI LANDSAT-8 sensor, previously used for fires detection [43][44][45]. The algorithm is essentially divided into three main steps focused in detecting all the hot-spot-contaminated pixels present in the analyzed image ( Figure 3). In the first step, a series of logical test is used to identify pixels potentially contaminated by hot spots. The second step calculates a Thermal Index (hereby considered a proxy of the pixel integrated temperature) for each alerted pixel and group the adjacent pixels into distinct clusters. In the third step, each cluster is analyzed spatially and statistically in order to reduce effects of thermal halo, diffractions spikes or "blurring" (adjacency effects, i.e., the scattering of light reflected from the nearby land into the sensor's field of view, due to contrast between a target pixel and its neighborhood; [56]).

Step 1: Spectral Principles
This initial step is an implementation of the Hotmap algorithm [44], developed to detect hot-spot-contaminated pixels within OLI-LANDSAT 8 daytime images.
Accordingly, hot-spot-contaminated pixels are individuated by using reflectance ratios and single-band thresholds grouped to define four distinct logical tests, named α, β, S and γ: A pixel fulfilling at least one of these tests is flagged as hot-spot-contaminated and is taken into account during the further steps.
In this procedure, we keep the original Hotmap approach, where tests α and β are designed to (i) reveal very few false alarms and detect at least one hot "warm" pixel (α condition), and (ii) identify "hot" pixels (β condition). However, we modified the band ratios' thresholds in order to better avoid false alarms over highly reflective surfaces (i.e., clouds and dry soils): parameters α and β are structured in a first-ratio portion (ρ12 must be 40% and 20% greater than ρ11 and ρ8a for α, respectively; ρ11 twofold of ρ8a for β) and in minimal conditions in the three bands, to avoid too-low reflectances in bands 12, 11 and 8a.
In addition, we added two other conditions which allow us to consider pixels having α and/or β corrupted by saturation of bands 12 and 11 (S condition), as well as very reflective pixels (γ condition), surrounded by other alerted pixels (satisfying α and β). The thresholds for S test condition represent the saturation values, retrieved empirically by analyzing the limits of reflectance values for the SWIR bands in different volcanic thermal active bodies. The γ condition triggers just the very hot detected pixels, which are placed into an area surrounded by α and β parameters; this is necessary to detect pixels located in the inner part of very-high-radiance-emitting events (such as large lava flows or very thermally active lava lakes) otherwise discarded by the previous tests. It should be noted that some single-band thresholds in condition S and γ are higher than 1: this condition is necessary since specular effect on land surface and clouds can lead to reflectance values of S2 higher than 1 [51]. The spectral thresholds proposed for Step 1 are derived by using an empirical approach, starting with those band ratios proposed by previous hot-spot-detection systems [43] and testing on different volcanic cases (as shown later) the best achievable combination both to avoid non-volcanic false alarms and to detect low thermal anomalies. The output of Step 1 is thus a logical matrix named ALERT (composed by 0 and 1), having the same dimension of the S2 cropped image (501 × 501 pixel), and containing all the pixels classified as alerted (with value = 1). If the ALERT matrix is composed only by zeros, no alerted pixels are founded.

Step 2: Thermal Index
The ALERT mask is scanned in order to group all adjacent pixels potentially alerted into distinct clusters, identified as each connected component in the binary image. As explained above, it is not practical to use the spectral radiances to accurately calculate the area or integrated temperature the sub-pixel hot emitter because, over volcanic targets, the bands 12 and 11 are often saturated. For this reason, in this work, we introduced an empirical Thermal Index (TI), defined as the sum of the reflectance recorded in bands 8, 11 and 12: Thermal Index(TI)= ρ8a + ρ11 + ρ12 Remote Sens. 2020, 12, 820 8 of 32 Here, the TI is considered to be a proxy for the pixel integrated temperature, based on the assumption that increasing the size or temperature of the sub-pixels' hot target will produce an increase in the reflectance values and hence in the TI. This Thermal Index is calculated for all the pixels alerted during Step 1. Each cluster is then analyzed separately from the others during the Step 3.

Step 3: Spatial and Statistical Principles
Step 3 consists in the spatial and statistical analysis of TI values, applied to each cluster identified in Step 2 (see scheme in Figure 2). The aim of this step is to remove, or at least reduce, the number of alerted pixels caused by "blurring", diffraction spikes and/or halo effects, which in some cases distort the effective spatial pattern of the anomaly and cause an evident increase in the size of the alerted clusters (i.e., excess of hot area, which may occur in presence of clouds or large thermal anomalies), especially in case of intense thermal anomalies.
Due to their small size, all clusters composed by less or equal than 9 pixels (3×3 box) are considered unaffected by the above-described effects and are immediately classified as hot-spot-contaminated, independently on the TI values. This condition allows us to detect and keep unchanged small-size thermal anomalies associated to hot degassing cracks or eventually occurring at the beginning of a renovated volcanic activity phase (i.e., Figure 4a-d). For instances, fumaroles at Chaitén lava dome (Chile, [57]) are perfectly localized by our algorithm, which detect several small (< 9 pixels) clusters composed by α-triggered warm pixels only (Figure 4 a,b). Instead, glowing inside the main vent of Anak Krakatau Island (Indonesia, [58]) triggered a central hot pixel (β condition) and a surrounding perimeter of eight warm pixels (α condition), related to a new Strombolian activity producing incandescence and grey plumes at the beginning of the eruptive phase culminating with the December 2018 island flank collapse (Figure 4c,d; see [59]). These cases definitively reveal the meaning of the 9 pixels threshold, working as a minimal criterion to detect thermal signals produced by a very hot and localized source, with a hotter core extending heat in the nearby area of at least of eight pixels. As discussed later, this feature proves to be extremely useful to detect volcanic precursor thermal signals and/or weak but persistent hot spots.
If the cluster is composed of more than 9 pixels, the algorithm applies an investigation of the frequency distribution of TI values (see flow chart in Figures 3, 5a and 6a). In particular, it calculates the arithmetic mean of TI (TI mean ), the 30th percentile of TI (TI 30% ), and recognizes the value of max departure of the observed distribution from the normal one (TI flex ), for each cluster. The frequency distribution plot of the Thermal Index allows us to observe how the different pixels are thermally distributed and where the observed distribution diverges from the theoretical normal one (Figures 5a and 6a). The purpose of this step is to automatically recognize the lower tail of the observed thermal distribution, by defining a "flex" threshold where the distribution of hottest pixels in the cluster moves from the less hot pixels. This analysis allows us to automatically recognize the colder pixels in the lower tail of the TI distribution, which we ascribe to the thermal halo or "blurring effects" surrounding the real hot spots. These pixels are detected by defining a contextual TI threshold (TI thres ), which varies from cluster to cluster, depending on the analyzed TI distribution and according to the following conditions: TI thresh = TI flex (for TI flex < TI mean ) TI thresh = TI 30% (for TI flex > TI mean ) where TI flex is the TI value exhibiting the largest difference between the observed and the theoretical normal distribution (Figure 5a).
(see Figure 6a). In these cases, where the TIflex > TImean (Equation (7)) in the cluster (green dotted line and black line respectively in Figure 6a), we choose a conservative and fixed solution, setting TIthresh at the 30th percentile of thermal distribution (TI30%, yellow dotted line in Figure 6a). This threshold is set empirically because we observed in lots of thermal distribution cases, particularly those related to limited extended lava bodies, that the 30% simulates a good fitting cut in the divergency between the lower tail and the normal theoretical distribution (cross between blue pixels distribution and the normal theoretical distribution black dotted line Figure 6a).   In red and blue, hot selected and discarded pixels. Dotted black, solid black, dotted yellow, dotted green and red lines represent, respectively, the theoretical normal distribution model, TI mean , TI flex , TI 30% and TI thres hot threshold, fitting with the TI flex in this case. (b) SENTINEL-2 RGB (12-11-8a) image zoomed over cluster analyzed; in green, red and withe, alfa-triggered, beta-triggered and discarded pixels; respectively. Distribution plot of the Thermal Index (TI) of the analyzed cluster. In red and blue, hot selected and discarded pixels. Dotted black, solid black, dotted yellow, dotted green and red lines represent, respectively, the theoretical normal distribution model, TImean, TIflex, TI30% and TIthres hot threshold, fitting with the TIflex in this case. (b) SENTINEL-2 RGB (12-11-8a) image zoomed over cluster analyzed; in green, red and withe, alfa-triggered, beta-triggered and discarded pixels; respectively.  In red and blue, hot selected and discarded pixels. Dotted black, solid black, dotted yellow, dotted green and red lines represent, respectively, the theoretical normal distribution model, TI mean , TI flex , TI 30% and TI thres hot threshold, fitting with the TI flex in this case. (b) SENTINEL-2 RGB (12-11-8a) image zoomed over cluster analyzed; in green, red and withe, alfa-triggered, beta-triggered and discarded pixels; respectively.
We define the two conditions expressed in Equations (6) and (7), because the clear presence of a single sharp flex in the distribution of a thermal cluster able to cut the lower emitting pixels is not unequivocal in all cases. This could be mainly due to a limited extension of the hot body, which may bring to a complex thermal distribution with different potential "flexes" (Figure 6), or to a masking effects of thermal anomalous region by clouds/plume presence. Indeed, the process could lead to the definition of a too-high thermal flex, which does not differentiate the lower tail but the higher one (see Figure 6a). In these cases, where the TI flex > TI mean (Equation (7)) in the cluster (green dotted line and black line respectively in Figure 6a), we choose a conservative and fixed solution, setting TI thresh at the 30th percentile of thermal distribution (TI 30% , yellow dotted line in Figure 6a). This threshold is set empirically because we observed in lots of thermal distribution cases, particularly those related to limited extended lava bodies, that the 30% simulates a good fitting cut in the divergency between the lower tail and the normal theoretical distribution (cross between blue pixels distribution and the normal theoretical distribution black dotted line Figure 6a).
The benefit of the frequency distribution analysis is that it does not act as a single threshold over the entire image, but works with a contextual filter based on the thermal distribution of each cluster itself. This means that, for instance, the threshold value for a high thermal anomaly produced by a widespread lava flow will not be the same as the one related to a breakout over an active lava-tube, to hot-degassing fumaroles fields or to hot materials at the surface of a growing lava dome. The discarding process of lower anomalous pixels based on the flex definition cuts the colder boundaries surrounding the hotter core and defines a more-fitting shape of thermal anomalies, as shown in Figures 5b and 6b (where white represents the alerted pixels discarded; green and red are the α-triggered and β-triggered conditions pixels, respectively).
The differentiated-steps procedure described above was designed to work on a wide spectrum of volcanic activity, from hot fumaroles to lava domes, lava flows and lava lake. It has a high sensitivity to low, small thermal anomalies, which is extremely useful to detect both volcanic precursor signals and/or weak persistent hot spots. It is also able to effectively isolate the hottest pixels in wider clusters, discarding anomalous pixels alerted by thermal halo or clouds refraction effects.
In this regard, Figure 7 displays two explicit cases to show how the algorithm works overall and its potential for volcanological purposes.

Results
The reliability of the algorithm described above was tested on several different volcanoes, first by quantitively evaluating the algorithm performance to detect thermal anomalies and to avoid false alerts, and then by comparing the number of hot pixels detected by the algorithm with MODISderived radiant-heat-flux timeseries. Here, we present the results of these analysis.
We investigated eight volcanoes case studies, located in very different geographic contexts and characterized by four exemplary end-members volcanic heat sources: lava flows, lava lake, lava domes and open-vents. The eight case studies are examined from the 1st January 2016 to 1st October 2019 timespan: Kliuchevskoi (Kamchatka peninsula, Russia; [64]) and Etna (Sicily island, Italy; [61]) for lava-flow activity, Erta Ale (Afar region, Ethiopia; [60]) and Masaya (Nicaragua; [65]) for lava-lake activity, Stromboli (Sicily island, Italy; [62]) and Villarrica (Chile; [66]) for open-vent activity in lowviscosity systems, and Bezymianny (Kamchatka peninsula, Russia; [67]) and Láscar (Chile; [68]) as lava-dome type with high-viscosity lavas extruded (cf. with Figure 1). The algorithm is therefore tested on a wide variety of volcanic systems, in order to assess its performance and limits in detecting various thermal sources. Because the chosen volcanoes show persistent but variable activity, we tested the reliability of this latter on detecting transitions between different eruptive styles over time. Regardless, geographical locations of the volcanoes are variable (Figure 1), with diverse contributions of weather conditions and cloud-coverage effects. A small-sized light-brown plume, probably of volcanic origin, overlays the anomalies triggered by summit strombolian activity at Stromboli volcano [62] during January 2018. The presence of this semitransparent cloud affects the thermal detection, because it produces a halo effect, expanding the possible hot-spotted pixels (Figure 7a). The contextual algorithm here proposed recognizes a clear TI flex in the thermal distribution of the cluster and excludes all pixels interested by cloud reflections, thus precisely isolating the hottest portion of the cluster, mainly colored with yellow to orange and bright red tones (Figure 7c). Figure 7b shows the performance of the Hotmap algorithm [43], which overestimates by about 4 times the hot area detected by our algorithm (Figure 7c).
The strong thermal anomalies occurring over the Yasur's pyroclastic cone [63], composed by two distinct active vents, indicate a sustained strombolian activity or spatter producing incandescent bombs as far as the crater rim, underscoring the hazardous nature of get closer to the volcano summit ( Figure 7d). The extremely powerful heat-release affects the MSI sensor and induces it to produce a "cross-shaped" thermal effect, called diffraction spikes, with four arms radiating from the bright volcanic sources in the SENTINEL-2 image. We noticed that these artifacts occurred commonly in very energetic and saturating thermal anomalies. As we can see in Figure 7e,f, the here-proposed contextual algorithm, compared to Hotmap algorithm, reduces the number of hot pixels related to these diffraction wings, detecting hot materials just into the two crater areas. Therefore, the concurrence of both spatial and statistical filters (Step 3 in the algorithm) particularly allowed to exclude hot-spot pixels not directly related to the volcanic hot spot (i.e., clouds coverage) or triggered by instrument optics effects combined to intense thermal emissions (i.e., diffraction spikes). This enhancement is extremely relevant for monitoring purposes, particularly for automated system applications, since pixels affected by halo effects might have been wrongly interpreted as hot materials emplacing outside the summit crater terrace of Stromboli or the Yasur crater rim, with relevant implications for hazard assessment.

Results
The reliability of the algorithm described above was tested on several different volcanoes, first by quantitively evaluating the algorithm performance to detect thermal anomalies and to avoid false alerts, and then by comparing the number of hot pixels detected by the algorithm with MODIS-derived radiant-heat-flux timeseries. Here, we present the results of these analysis.
We investigated eight volcanoes case studies, located in very different geographic contexts and characterized by four exemplary end-members volcanic heat sources: lava flows, lava lake, lava domes and open-vents. The eight case studies are examined from the 1st January 2016 to 1st October 2019 timespan: Kliuchevskoi (Kamchatka peninsula, Russia; [64]) and Etna (Sicily island, Italy; [61]) for lava-flow activity, Erta Ale (Afar region, Ethiopia; [60]) and Masaya (Nicaragua; [65]) for lava-lake activity, Stromboli (Sicily island, Italy; [62]) and Villarrica (Chile; [66]) for open-vent activity in low-viscosity systems, and Bezymianny (Kamchatka peninsula, Russia; [67]) and Láscar (Chile; [68]) as lava-dome type with high-viscosity lavas extruded (cf. with Figure 1). The algorithm is therefore tested on a wide variety of volcanic systems, in order to assess its performance and limits in detecting various thermal sources. Because the chosen volcanoes show persistent but variable activity, we tested the reliability of this latter on detecting transitions between different eruptive styles over time. Regardless, geographical locations of the volcanoes are variable (Figure 1), with diverse contributions of weather conditions and cloud-coverage effects.

S2 Algorithm Evaluation
The hot-spot-detection algorithm has the main aim to detect the most of thermal anomalies produced by volcanic activity, by minimizing errors related to false alerts, which is crucial to performing operationally over a global scale. To define when a pixel is volcanologically alerted or not is not a univocal concept [43], and this kind of analysis brings a certain degree of subjectivity and relativity. In this section, we present the performance of the proposed S2 algorithm based on the visual inspection all the SENTINEL-2 images, for a total of 2.211 products, acquired during the January 2016 to October 2019 period, over the abovementioned eight volcanic case studies.
According to our visual inspection, we classify all the "alerting conditions" in the following ways: (i) True volcanic alert: a pixel anomaly detected by the algorithm which is expressly and visually related to volcanic activity (hot degassing, lava body exposed, hot eruptive materials exposed and possibly confirmed by literature or consistent with the background knowledge of volcano activity) and that shows a visible thermal glowing (from dark-red to bright-white in color); (ii) Fires or anthropogenic alert: a pixel anomaly detected by the algorithm expressly and visually related to wildfire occurrence and/or located near of human-settled areas (usually these latter anomalies could be far from the source of volcanic activity, show dark-red tonalities and have a very small areal extent of 1-2 pixels); (iii) False alert: a pixel anomaly detected by the algorithm, visually related to cloud coverage, secondary cloud fringes and reflections effects over hot sources (also volcanic), blurring, high-reflectivity effects of land/sea surfaces; (iv) S2-MSI-derived false alert: a pixel anomaly detected by the algorithm and visually triggered by abnormal colored pixels related to stripe artifacts of the MSI sensor.
Notably, one "alerting condition", as defined above, cannot exclude the occurrence in the same S2 image of one or more others "alerting condition(s)", so that different kinds of alerts could take place concurrently. For a successful detection, the algorithm must detect a true volcanic alert, and this capability underlines the reliability of the algorithm to recognize thermal volcanic activity. On the contrary, the algorithm inaccuracy is expressed by how many false alerts have been counted in the overall number of alerts detected. In Table 2, we summarize the evaluation results.
In the upper part of the table, results about detection capabilities of the algorithm are presented. The successful detection percentages (sixth column in upper part of Table 2) show high values for all the case studies, with percentages spanning between ca. 85% and 97% (Bezymianny and Kliuchevskoi as minimum, Masaya as maximum case studies), highlighting the excellent algorithm accuracy on sensing thermal volcanic activity. The main controlling factor appears to be the weather conditions, with the only two percentages below 90% related to the Russian Kamchatka's volcanoes, located at high-latitudes and with a robust cloud cover contribution. Particularly at Bezymianny and Etna, the algorithm missed diverse still visible volcanic hot spot due to a commonly thick cloud coverage, in first case mainly, and to persistent white steam and degassing over summit vents, partially hiding hot anomalous pixels. Notably, fire/human-induced alerts are excluded from this analysis, because they do not represent volcanogenic thermal anomalies, but nevertheless they are not considered as false, because they are likewise triggered by hot-spot sources, thermally comparable with volcano-related ones. Understandably, fires/human alerts are commonplace in Etna and Stromboli, where fires are likely and human activities are located on the volcano's slopes, and Masaya, where the nearby villages surrounding Masaya city fall into 5 km from volcano crater.
In the second part of the table, false-detection amounts per volcano are listed, distinguished between false alerts and MSI artifacts. Algorithm false detection percentages (second column in lower part of Table 2) exhibit generally extreme low values, spanning from no fake alerts to 3.45%, with the two higher percentages are to be found on Etna and Kliuchevskoi. In some case studies (e.g., Erta Ale, Masaya, Villarrica and Láscar), no false alerts are detected, and this is mostly due to absence of cloud coverage and its related disturbing effects. Higher percentages of false detection occur if we also consider fake effects induced by artifact on MSI detector sensing. These effects, known as "spectral response non-uniformity" [69], consist in oblique soft-edged darker or brighter stripes near the detector margins, that induce localized irregular spectral response and anomalous pixel colors, some of that of red tones. Particularly, fake anomalous pixels are produced when stripes interact with clouds coverage, and this issue is more noticeable still in Etna and Kliuchevskoi cases and absent in ones with cloud-free conditions. These artifacts are tiles and sensor dependent and for that reason does not show the same occurrence in the geographical areas of volcanic case studies analyzed.
The evaluation results summarized in Table 2 indicate an overall great reliability of the S2 algorithm here presented for working on a global-scale level, considering the eight different volcanic case studies and their variabilities in thermal activity and geographical location. The remarkable successful detection percentages underline the high sensitivity of the algorithm to detect different thermal hot-spot sizes and natures. Moreover, if we do not consider the MSI-related artifacts, which are not dependent by how algorithm works and neither avoidable nor finally resolvable by our application, false-detection occurrence has certainly low frequencies. Expectedly, the most affecting factor algorithm turns out to be the cloud-cover-related effects, both in true-and false-alerts detection. Table 2. Results of evaluation of the SENTINEL-2 hot-spot algorithm. 1 N. images analyzed; 2 N. images with at least one alert detected by the algorithm; 3 N. images with at least one true volcanic alert detected by the algorithm; 4 N. images with at least one true volcanic alert; 5 N. images with at least one true volcanic alert missed by the algorithm; 6 Percentage of successful detections = (Algorithm volcanic alerts/True volcanic alerts) * 100; 7 N. images with at least one fire/human-related alert detected by the algorithm; 8 N. images with at least one false alert detected by the algorithm. 9 Percentage of false detections = (Algorithm false alerts/Algorithm alerts) * 100; 10 N. images with at least one false alert related to S2-MSI sensor artifacts detected by the algorithm; 11 Percentage of false detections considering also the S2-MSI derived false alerts = ((Algorithm false alerts + S2-MSI derived false alerts)/Algorithm alerts) * 100.

SENTINEL-2 and MODIS-MIROVA Timeseries Comparison
Here, we show the almost four-year-long timeseries comparison between the number of hot pixels, S2Pix, detected by the algorithm application on SENTINEL-2 MSI images and the Volcanic Radiative Power elaborated by the MIROVA system on MODIS data, in order to explore the accuracy of the algorithm to detect and track volcanic thermal activity.
MIROVA is an automated global hot-spot-detection system based on the near-real-time processing of infrared MODIS-MIR data (http://www.mirovaweb.it/; [25]). Its thermal algorithm incorporates spectral and spatial filters, providing an accurate sensitivity to lower thermal anomalies. It gives a quantification of the Volcanic Radiative Power (heat flux, in Watt) through a hybrid algorithm based on MIR radiance data analysis (~3.9 µm) recorded at the moderate spatial resolution of 1 km per pixel of MODIS satellites images. MIROVA currently monitors over 220 volcanoes, providing real-time post-processed products, such as lava effusion rates, in support to eruptive crisis management [20]. Figures 8-11 show the heat flux in Watt calculated by MIROVA system (on the left y scale, with blue stem, from 10 5 W to 10 10 W), and the number of hot pixels S2Pix detected by the new algorithm (on the right y scale, with red dots, from 1 to 10,000) on the eight volcanic case studies. No atmospheric correction was performed on S2-MSI and MODIS images here presented, in order to evaluate the purest detection capabilities. The MODIS dataset is filtered to exclude images with poor viewing geometry (Zenith > 40 • ) and possible deformation effects of the projected thermal anomalies (see [25]). The MIR method [70] applied by MIROVA to MODIS images detects the thermal flow radiated from the surfaces with T > 500 K solely and return the VRP with a standard error of 30%. VRP values represent just the thermal output of high-temperature features, produced by the arrival of magma at the surface or at very shallow levels [20]). The comparison with MIROVA data is therefore relevant for the S2 algorithm proposed here, because both aim to detect the hottest portion, directly related to the ascent of magmatic fluids subaerially exposed. The first measures the volcanogenic heat flux, and the second measures the number of hot pixels, i.e., the radiant emitting hot area.
Lava flows represent the most prominent example of volcanic heat sources detectable by satellite, as depicted in Figure 8a-h, showing Etna and Kliuchevskoi volcanoes.
In the case of Etna, whose thermal activity varied by several orders of magnitude during the analyzed period (Figure 8a), it is possible to observe how the detected S2Pix follow the VRP thermal trend with high accuracy and mimic its general trend. In fact, there is a strong correspondence between the two analyzed thermal signals, both during the so-called "low thermal regime", related to strombolian and degassing activity at the summit craters (characterized by VRP < 10 8 W (cf. [37]) and S2Pix <~100, Figure 8c), and during the "high thermal regime" related to high-sustained strombolian activity or major lava effusion phases (Figure 8b,d and yellow bands in Figure 8a, characterized by VRP > 10 8 W and S2Pix > 100). Even if energetic events such as summit overflows or fountaining episodes at Etna are generally short-lived, each peak detected by MIROVA is also constrained by S2Pix hot pixels detections, also thanks to 2-3 days of revisit time of SENTINEL-2 above Etna.
Kliuchevskoi volcano is another good example of lava flow-forming eruption, with an excellent correlation between MODIS and S2 related thermal signals, despite the high latitude and strong cloud coverage in the region. In this case, the main thermal activity was recorded between March and November 2016, with a rising trend related to an increased lava effusion from the summit craters (max VRP > 10 9 W, S2Pix >1000; Figure 8e,f,g). A sharp drop in the thermal activity tracked by MODIS and S2 marked the end of the eruptive phase on November 2016, and it was followed by a month-long slow decline of thermal emissions detected by MIROVA (VRP < 10 6 W; Figure 8h). Despite the low revisit frequency at these latitudes, spanning from 5 to 10 days, the effusive phase and its increasing trend was well detected by the algorithm, with a growing number of detections from hundreds to a peak of 1200 S2Pix (Figure 8e). Afterward, the lack of a continuous satellite acquisition and the more sporadic thermal activity of Kliuchevskoi makes the S2 data sparse, although sporadic hot spots were detected during short-lived, low-level resuming phases (May-June 18 and July-August 19; Figure 8e,h).  Lava-lake cases are represented by Erta Ale and Masaya volcanoes, as shown in Figure 9a-h. Erta Ale volcano is well-known for its long-lived lava-lake activity hosted within one (or two) summit pit-crater(s) [71]. However, since January 2017, a large flank eruption (still ongoing at the time of this writing), has produced lava flows, with an estimated area of at least 26 km 2 , in June 2019 [72]. Considering the drastic change in the eruptive style, the consistency between MIROVA and S2 thermal signals is remarkable both in terms of recorded intensity and trends (Figure 9a). The correlation between the two signals is very high both during the lava-lake activity (before January 2017; Figure 9b), with S2Pix < 100 and VRP < 10 8 W, except some overflows, as well as during the lava flows' production (after January 2017; Figure 9c,d), with a peak phase in S2Pix (max > 1000 pix.) and VRP (> 10 9 W). Some S2Pix values are underestimated due to cloud coverage that partially masked the satellite sensor detection (S2Pix~10, visual inspection of those data). Moreover, a series of short-lived overflows occurred during the second half of 2016 and the beginning of the effusive eruption on 22 January 2017 (VRP > 7000 MW) were not detected by S2. This clearly highlights the limitations of this system in detecting short-lived yet very intense events. The high number of S2 hot-spot detections and the strong correlation observed with VRP values are even favored by the good atmospheric conditions in this African area, with a very low number of cloudy periods and, thus, optimal conditions satellite observations. Remote Sens. 2020, 12, x FOR PEER REVIEW 18 of 33   A similar behavior can be seen in the Masaya case study (Figure 9e), where a persistent thermal output has been detected by MIROVA since the resumption of the lava lake in December 2015 [8]. Here the S2Pix (generally between 10 to 200 pixels) maintains a good overall match with VRP (~10 7 -10 8 W), tracking the formation and rise of the lava lake until March 2017 (Figure 9e,f,g), followed by a gentle decline until now; notwithstanding, the S2 hot spots appear more scattered and, in some cases, seems to be higher to the corresponding VRP (Figure 9e). This discrepancy could be partially due to (i) a greater influence of cloud coverage in scattering, masking or, at times, increasing of hot spots due to reflection effects by low and thin clouds (blurring effects not completely removed), and (ii) an inner variability in thermal output emission of the Masaya lava-lake system, with different energetic phases related to a rapid fluid dynamic behavior rather than a more stable and slower surface observable at Erta Ale, as testified by some works on Masaya lava-lake dynamics (see [73] and references therein).

Discussion
The relationship between the number of S2 pixels and the VRP heat flux provides important information on the thermal structure of sources and related volcanic processes at the origin of the anomalies detected by SENTINEL-2 MSI and MODIS. For each of the eight cases described above, we selected only the cloud-free S2 images (through visual inspection), and we considered the Open-vent activity case studies are represented by Stromboli and Villarrica volcanoes, in Figure 10a-h. Both cases present interesting results, considering that open vents volcanoes should be characterized highly variable activities, with continuous gas emissions, mild explosive behavior and with the magma column top often visible and considered as the last portion of the shallow plumbing system [74].
At Stromboli, our algorithm detects almost persistent thermal anomalies associated to the continuous but variable strombolian activity (May-August 2017, December 2017, January 2018) as inferred by MIROVA-derived VRP (Figure 10a). The S2Pix and the VRP show excellent agreement during effusive activity (July-August 2019, Figure 10d), with the threshold of~50 MW (see [75]), or 50 S2Pix, representing the transition between strombolian and effusive regimes (summit overflows and/or flank effusions, in Figure 10a as yellow bands and Figure 10d). During periods of weak strombolian activity the high sensitivity of our algorithm to very low thermal emissions (January-March 2016) is testified by few hot pixels almost continuously detected (S2Pix < 10), with no or weak MIROVA thermal anomalies (VRP~10 6 W; Figure 10b,c). It is also clear how the thermal signal produced by strombolian open-vent activity results in a variability in the hot area detected by S2, spanning from a limited S2Pix < 10 to a dozen of hot-spots.
At Villarrica volcano, the MIROVA VRP measurements indicate persistent thermal emissions with values comprised between 10 7 and 10 8 W (Figure 10e-h), with variations in the volcanic heat flux probably related to the fluctuations of the lava level and intensity of strombolian activity [76]. Until December 2017, the S2 algorithm successfully recognizes this thermal trend, with a variable hot-spots number, from a few hot pixels to S2Pix > 100, that nicely fit the VRP trend (Figure 10e). Afterward, the VRP shows at least two major cycles of decreasing and increasing thermal emission (Figure 10e). In relation to these cycles, the S2Pix appears quite stable with a value around 10 pix., whilst they also track the higher thermal emitting phases of June-September 2018 (Figure 10g) and of September 2019. This suggests that, during the lowering phases, the area occupied by the thermal anomaly remains roughly constant (i.e., the bottom of the crater), while the thermal flux decreases greatly, in response to the lowering of the effective temperature of the hot target.
A similar behavior is thus found in these two open vent cases, with a strict correlation between VRP-S2Pix during major activity phases and the sensibility to recognize hot emitting area by S2 even when VRP clearly decreases.
Considering the lava-dome activity types, or more generally open-vents activity with high-viscosity magmas [74], we analyze Bezymianny and Láscar volcanoes, in Figure 11a-h. During the analyzed timespan, Bezymianny activities consist of sustained fumarolic emissions, periodically interrupted by phases of lava-dome growth, strong explosions and eventually by extrusion viscous lava flows (Figure 11a-d). This variability is well represented by VRP data, showing a "thermal baseline" (VRP < 10 7 W), and peaks in thermal flux occurring when lava domes/flows or hot pyroclastic deposits are exposed on the volcano summit and flanks (Figure 11a-d). This activity is well tracked also by the S2Pix, which follows the same overall trend, with particularly good accuracy during the higher activity phases (i.e., peak 16 March 2019, with~7 × 10 8 W of VRP and S2Pix = 515, where hot material is dispersed over the top of volcano and surrounding slopes Figure 11d). During lower thermal emission periods or decreasing phase, such as after a strong explosion occurred in June 2017 (Figure 11b), the S2Pix agrees with the overall VRP trend, but with a clear dispersion due to the high cloud coverage in the region.
Instead, over Láscar dome volcano, we find an extremely stable and persistent thermal signal, partially owed to excellent weather conditions in that area, detected by both MIROVA and S2 algorithm. Here, the S2Pix mimics very well the VRP trend both in the 2016-2017 phase, (with a VRP < 10 7 W and a S2Pix < 10 pixels; Figure 11e), and during phases of increased thermal activity, as occurred from November 2018 (Figure 11h). It is interesting to note that, during the first half of 2018, no thermal activity was detected by either of the systems (Figure 11g). Moreover, the S2 algorithm was able to detect the sharp rise in the thermal signal produced by Láscar lava dome during October to November 2018 (from S2Pix~3 to S2Pix > 20 pixels; Figure 11h).

Discussion
The relationship between the number of S2 pixels and the VRP heat flux provides important information on the thermal structure of sources and related volcanic processes at the origin of the anomalies detected by SENTINEL-2 MSI and MODIS. For each of the eight cases described above, we selected only the cloud-free S2 images (through visual inspection), and we considered the maximum VRP recorded by MODIS images (filtered data with Zenith angle < 40 • ) in a time window of ±24 h, from each S2 acquisition. This allowed us to associate the S2Pix to a VRP value measured almost simultaneously, limiting the discrepancies due to sudden changes in activity or cloud cover. The aim is to compare the two best satellite images, acquired when the hot spots and the heat flux should represent and measure the same thermal volcanogenic process. This kind of analysis was similarly applied by some authors to test the correlation between MODIS and other thermal satellite related signals (i.e., ASTER; see [35]).
In order to evaluate the correlation between SENTINEL-2 and MODIS-MIROVA data, we examined the eight precedent cases by comparing the intensity of the heat flux produced (VRP) and the hot area responsible for that thermal radiation (S2Pix).
In Figure 12, we plotted the S2Pix as a function of the corresponding VRP, with dashed lines representing "isotherm" curves (from 10 5 to 10 8 W/S2Pix). Considering that VRP represents the radiant flux in Watt and the S2Pix is a proxy of the hot area, the region of the plot where the observed data fall is an indicator of the temperature of the source, which is useful to assess the thermal behavior of the investigated volcanic phenomena. In this regard, for each case, we modeled the relationship (red lines in Figure 12) according to the following equation: where A and B are the slope (Watt per pixel) and intercept (minimum radiant power for one alerted S2 pixel) of the model. According to Equation (8), all the isotherm curves in the graphs have A = 1. Beginning from volcanoes dominated by lava flows, we observe how the relationship is represented by data that follow an isotherm trend and linearly grow in area and heat flux produced. Even if the numbers of the useful data are strikingly different, both Etna and Kliuchevskoi show this similar behavior (Figure 12a,b). In particular, for Etna, the distribution of the data is representative both of Strombolian (major clustering about 10-100 S2Pix and 10 6 -10 8 W) and effusive activity, with major hot pixels (100-1000 S2Pix) and VRP (>10 8 W). This fits with what is expected for the thermal signal of a lava flow: the higher the area occupied by the hot body is, the greater the radiative power is [25]. Kliuchevskoi cluster has a lower number of data points (few cloud free images) and is more dispersed, but it is strictly related to the effusion phase described in Figure 8e, where the VRP-S2Pix combination detects the rising in thermal signal and the consecutive decline following the cessation of lava-flow feeding.
At lava lakes, two different behaviors are observed at Erta Ale and Masaya (Figure 12c,d). This results from the fact that Erta Ale dataset is composed of two different eruptive styles and associated thermal regimes (cf. Figure 9a). A first regime of lava-lake thermal emission only is characterized by a lower VRP/S2Pix ratio and some higher VRP-lower S2Pix points associated to overflow events. Conversely, the higher thermal regime is related to the large effusive phase (reaching S2Pix > 10 2 and VRP > 10 8 W) and is represented by the dense cluster with isothermal behavior, not by chance, as in Etna and Kliuchevskoi cases (Figure 12a-d). Masaya shows a different thermal behavior, with data clustered in the same region, around S2Pix~10 2 and 10 7 W < VRP < 10 8 W, and not aligned along an isotherm. This feature suggests the presence of a thermal source (the lava lake), essentially confined within the deep crater, with very limited variations in both the hot area and the heat flux.  Lava dome examples draw clusters in the lower leftmost portion of the plots (Figure 12e,f), thus indicating small-size anomalies with an overall low heat radiation involved in dome-forming eruptions, such as degassing, surface cracks on the dome body, collapse events or extrusion phases. The Bezymianny volcano exhibits thermal behavior with most detections having a VRP < 10 7 W and S2Pix < 20, except for a few points related to explosive events or post-explosive exposure of hot ejecta (Figure 12e). Láscar displays two well-clustered point clouds, the first being sparser, around~10 6 W, and the second being denser, with a VRP up to 10 7 W (Figure 12f), that in turn represent the two thermal regimes already observed in the timeseries (Figure 11e). The overall distribution seems partly to cut the isotherms, indicating that periods of higher heat flux are not accompanied by a proportional increase in the hot area exposed.
At open-vent systems, we observe the most scattered distributions among the case studies. Both Stromboli and Villarrica seem to draw a trend with a sharp increase in heat flux, from about 10 5 to VRP > 10 8 W, still maintaining a relatively small amount of hot area exposed, generally with S2Pix < 10 2 (Figure 12g,h). Only a few detections at Stromboli seem to represent a "isothermal" linear relationship between VRP and S2Pix, that are in fact related to an effusive phase occurred during July-August 2019 (cf. with Figure 10a). This behavior could represent an inner character of open-vent volcanoes, where an intensification in the VRP could be related to the increase in the convective dynamics in the upper portion of the magma column, inducing a rising of the magma column level and an increase in vent(s) temperature [20,75], while retaining a hot area quite constant.
These considerations allow us to outline some key points about the algorithm here proposed and its reliability in the detection of volcanic hot spots.
First of all, we observe that the SENTINEL-2 hot-spots number has a high correlation with MIROVA analysis (R 2 > 0.8), particularly in low-viscous lava flows cases (i.e., Etna, Erta Ale and Kliuchevskoi; Figure 12a-c), where hot emitting bodies have the possibility to expand over large areas, to emit greater amount of heat flux and, thus, to trigger the sensitivity both of MODIS and MSI/S2 detectors, even though theirs spectral, spatial and revisit time differences.
Conversely, a diverse relationship takes place, when the heat source is partially morphologically constrained, and when the radiative processes are colder or more subtle to detect.
If the number of hot pixels cannot grow due to morphological constrain, the relationship necessarily diverges from an isotherm behavior and become more steep (increasing A parameter in Equation (8)). This occurs particularly at open-vent systems (i.e., Masaya, Villarrica and, in part, Stromboli examples; Figure 12d,g,h) where changes in the heat flux measured by MIROVA (A > 1) are not accompanied by the same evident variations in the hot area exposed to the atmosphere (the vent(s)). Thus, the increasing/decreasing in A value could, in turn, be useful to prove a variability in the temperature, and the consequent emitted heat flux, of the volcanic source.
When the emitter has a small size dimension, such as for lava domes, the VRP/S2Pix distribution substantially falls in the lower left region of the plots, indicating a clear colder thermal behavior, as testified by Bezymianny and Láscar cases. The points distribution can extend to higher VRP and higher S2Pix (following a more linear rule) only during major extrusive phases (Figure 12e,f) or in correspondence of pyroclastic flows produced generally after strong explosions.
Therefore, the proposed analysis turns out to be a useful approach to explore the thermal behavior and evolution of the volcanic processes. Accordingly, the different correlations observed and the relative values of R 2 and the occurrence of a linearity between VRP and S2Pix (A parameter) indicate variable relationships between temperature and hot area of volcanic hot emitting body. It is correct and meaningful, from a volcanological point of view, that the amount of VRP produced by a lava dome through hot fumaroles and that of an expanding lava flow by radiation, for examples, have a different related number of hot pixels and thus a different extension of the thermal source. This discrepancy represents the variation in temperature, and we suggest that this method of analysis could be potentially used to explore and classify the thermal signals of a variety of volcanic activities by space-based methods. The comparison between the number of hot pixels detected and the MIROVA signal is thus relevant to both validate the S2 algorithm detection capability and to infer useful insights about volcanic processes and their thermal features.

High-Spatial-Resolution Sensitivity
The most important advantage of SENTINEL-2 MSI imagery products is the high spatial resolution they offer. This feature strongly affects the detection capability and marks a sharp improvement compared to the moderate spatial resolution imagery, such as MODIS. In fact, sensors with higher spatial resolution are more skillful to sense small hot targets, because they cover a bigger proportion of the field of view of the detector element [43]. Detecting small thermal anomalies, minor volcanic events or changes in the locations of subtle hot targets could be extremely useful to notice thermal precursor signals and/or variations in weak but persistent hot spots. Here, we show two examples that clearly reflect this proficiency.
In Figure 13, timeseries of SENTINEL-2 S2Pix and MIROVA VRP are presented, for Merapi (Java Island and Indonesia; [77]) and Ol Doinyo Lengai (Tanzania, [78]) volcanoes. Merapi is the most active Indonesian volcano, located~25 km from the densely populated town of Yogyakarta, and hosts in a 200 m deep crater, a growing basaltic-andesitic dome regularly disrupted by explosions (up to VEI 4 level) causing pyroclastic flows [79]. Notably, an explosive eruption in 2010 killed about 300 people in the southern area of the volcanic edifice [80].

High-Spatial-Resolution Sensitivity
The most important advantage of SENTINEL-2 MSI imagery products is the high spatial resolution they offer. This feature strongly affects the detection capability and marks a sharp improvement compared to the moderate spatial resolution imagery, such as MODIS. In fact, sensors with higher spatial resolution are more skillful to sense small hot targets, because they cover a bigger proportion of the field of view of the detector element [43]. Detecting small thermal anomalies, minor volcanic events or changes in the locations of subtle hot targets could be extremely useful to notice thermal precursor signals and/or variations in weak but persistent hot spots. Here, we show two examples that clearly reflect this proficiency.
In Figure 13, timeseries of SENTINEL-2 S2Pix and MIROVA VRP are presented, for Merapi (Java Island and Indonesia; [77]) and Ol Doinyo Lengai (Tanzania, [78]) volcanoes. Merapi is the most active Indonesian volcano, located ~25 km from the densely populated town of Yogyakarta, and hosts in a 200 m deep crater, a growing basaltic-andesitic dome regularly disrupted by explosions (up to VEI 4 level) causing pyroclastic flows [79]. Notably, an explosive eruption in 2010 killed about 300 people in the southern area of the volcanic edifice [80].
Ol Doinyo Lengai (OL) volcano is famous for its unique active carbonatite volcanism, producing fluid and low temperature (495-590 °C) natrocarbonatite lavas [81]. Its recent summit thermal features included fumaroles, open vents, or small-scale cooling lava from pools or hornitos [82]. Some works already focused on the thermal analysis of OL activity using MODIS imagery, illustrating the limitations of original MODVOLC algorithm to routinely detect the low-intensity thermal anomalies [82,83] and underlining that these restrictions are mainly due to low temperatures, small area extents and low effusion rates of Ol Doinyo Lengai flows.  Ol Doinyo Lengai (OL) volcano is famous for its unique active carbonatite volcanism, producing fluid and low temperature (495-590 • C) natrocarbonatite lavas [81]. Its recent summit thermal features included fumaroles, open vents, or small-scale cooling lava from pools or hornitos [82]. Some works already focused on the thermal analysis of OL activity using MODIS imagery, illustrating the limitations of original MODVOLC algorithm to routinely detect the low-intensity thermal anomalies [82,83] and underlining that these restrictions are mainly due to low temperatures, small area extents and low effusion rates of Ol Doinyo Lengai flows.
Our analysis shows that the MSI-S2 20 m resolution images allow the detection of low thermal hot spots at these two volcanoes, otherwise undetected by MODIS 1 km resolution images (Figure 13). At Merapi (Figure 13a), the algorithm detected very low thermal anomalies (S2Pix = 1) since January 2016 and discontinuous but frequent hot spots during all 2017-2018 (max. S2Pix = 20, the 16/09/2018), when MIROVA system did not detect any (Figure 13b). The first MODIS anomaly is measured during December 2018 (very low < 10 6 W), followed during June-August 2019 by sporadic low anomalies (max. 2*10 6 W). In the analyzed period, the SENTINEL-2 SWIR thermal signal precedes by at least two years the first MIROVA detection, giving regarding the thermal activity persistence and locations during the dome evolution.
Similarly, at Ol Doinyo Lengai, the first SENTINEL-2 detection occurs in September 2016 and is followed by several others hot-spot detections, between 1 < S2Pix < 10, during June-September 2017, April-October 2018 (Figure 13c). Notably, the first thermal activity was detected by MIROVA in September 2018 and then intermittently during 2019 (Figure 13c,d). Here, too, the first SENTINEL-2 hot spot anticipated by almost two years the MIROVA detections, demonstrating the sensitivity of the algorithm to detect the small, low-temperature signals produced by the OL volcano.
Despite the lack of the thermally sensitive MIR bands of MODIS, the advantages using the high-spatial-resolution SWIR wavelengths to detect very hot surfaces are thus outstanding, because any hot spot will most likely cover a major portion of a SWIR pixel, resulting in a clear thermal signal detected.
These two examples demonstrate the improved ability of the proposed algorithm, to detect small and low thermal emissions which are undetected by systems like MIROVA. This ability is of great interest in order to identify possible thermal precursors at explosive and high-risk volcanoes, such as Merapi. Moreover, the low temperature lavas of Ol Doinyo Lengai shows the ability of the proposed algorithm to successfully detect volcanic thermal anomalies in a very wide variety of volcanic settings.

Reflection Effects and False Anomalies
Due to the high sensitivity of the algorithm to detect subtle and low-temperature thermal anomalies, some issues arise in specific cases, generally related to clouds. One of the most difficult cases is when dense and cirrus clouds partially overlap and/or surround volcanogenic-heat-emission sources. Cirrus clouds are thin and semitransparent, able to create thermal halo when covering the heat source. As demonstrated beforehand (see Section 2.3; Figure 7a-c), the proposed contextual algorithm isolates the hottest part of the triggered cluster, which is falsely enlarged by the thermal halo.
In the case of dense clouds covering intense thermal radiations, strong reflection, along with diffraction effects, may occur. In Figure 14a,b, a SENTINEL-2 RGB image of Nyiragongo lava-lake volcano (DR Congo, [84]) in the 8a-11-12 bands is shown, with the relative hot pixels triggered in the α and β conditions by the algorithm. The combined effect of the strong thermal emission produces a "cross-shaped" diffraction spike, and the presence of dense clouds is striking. Here, pixels obscured by clouds become deep red by approaching lava lake because of the sharp reflection and are triggered, mainly by the α condition, as hot pixels by the algorithm (Figure 14b). Moreover, the powerful thermal emission of Nyiragongo lava lake generates an enlargement of the apparent hot area, even when not covered directly by clouds, due to diffraction spike effect on MSI sensor (cf. Figure 7d-f). These pixels are also recognized as hot and alerted by the α condition. Though the thermal emission by Nyiragongo lava lake is not an artifact (the relative measured MIROVA heat flux is about 7 × 10 8 W), the result clearly overestimates the hot area dimension (Figure 14b). This issue results from the high sensitivity set for the α condition, that, however, as previously demonstrated, allows the detection of low thermal signal. Overall, the main aim being to offer a volcanic global applicable algorithm, able to detect hot spots, including hot-fumaroles and low-temperatures lavas, entails that partial overestimation can occur in specific settings, such as Nyiragongo. A possible solution is to use the SENTINEL-2 clouds mask [85] or to apply a more restrictive thermal algorithm just for a specific volcanic setting, for example limiting the detection to β condition triggered pixels in order to isolate the hottest portion of the detected hot area.
Remote Sens. 2020, 12, x FOR PEER REVIEW 27 of 33 restrictive thermal algorithm just for a specific volcanic setting, for example limiting the detection to β condition triggered pixels in order to isolate the hottest portion of the detected hot area.

Conclusions and Perspectives
We developed a new algorithm for detecting, locating and counting hot thermal anomalies in volcanic environments, with a global applicability, using SENTINEL-2 satellites images. The algorithm is based on a multispectral hybrid approach, whereby SWIR bands 12, 11 and 8a are spectrally, spatially and statistically elaborated, in order to enhance the presence of subpixel hot spots, with a spatial resolution of 20 meters that represents, nowadays, a remarkable profit for volcanic studies.
To explore the algorithm efficiency, the S2-derived number of hot pixel (S2Pix), detected at eight different volcanoes, was compared with the Volcanic Radiative Power (VRP) recorded by the MODIS-MIROVA thermal dataset. The results demonstrate an extremely coherent match in a variety of volcanic settings, involving hot bodies of various temperatures, spatial extent and typology. The match is particularly striking for wide intense thermal emissions such as lava flows, and more generally proves to give complementary information about thermal volcanic status, even if some limitations are still present mainly due to clouds coverage disturbs. Moreover, we have not confined our analysis to MSI SENTINEL-2 SWIR and MODIS-MIROVA MIR trend comparisons, but we explored what these correlations, between number of hot pixels detected (hot area) and the heat flux, could decrypt about thermal features related to different volcanic processes.
Significantly, the algorithm presented here is, as far we know, the first SENTINEL-2 multispectral based volcanic thermal detection process that runs operationally. In fact, the algorithm is part of the multiplatform MOUNTS volcanic monitoring system online since the beginning of 2018 (http://www.mounts-project.com/home; [30]). This system uses the SENTINEL constellation (-1, -2 and -5P) to retrieve and display key parameters volcano monitoring timeseries, such as deformation, heat anomalies and SO2 flux, in near real time.
In this work, we demonstrated, through an almost four-year-long investigation of eight case studies, how SENTINEL-2 thermal signal analysis can enhance the study and monitoring of several volcanic processes, including lava-flow morphometric evaluations, extrusion phases and growing dynamics of lava domes, lava-lake pulsations, fumarolic activity, multiple active craters thermal activity, thermal precursors and new hot spots' presence, or periods of magma column rising in openvent volcanoes.

Conclusions and Perspectives
We developed a new algorithm for detecting, locating and counting hot thermal anomalies in volcanic environments, with a global applicability, using SENTINEL-2 satellites images. The algorithm is based on a multispectral hybrid approach, whereby SWIR bands 12, 11 and 8a are spectrally, spatially and statistically elaborated, in order to enhance the presence of subpixel hot spots, with a spatial resolution of 20 meters that represents, nowadays, a remarkable profit for volcanic studies.
To explore the algorithm efficiency, the S2-derived number of hot pixel (S2Pix), detected at eight different volcanoes, was compared with the Volcanic Radiative Power (VRP) recorded by the MODIS-MIROVA thermal dataset. The results demonstrate an extremely coherent match in a variety of volcanic settings, involving hot bodies of various temperatures, spatial extent and typology. The match is particularly striking for wide intense thermal emissions such as lava flows, and more generally proves to give complementary information about thermal volcanic status, even if some limitations are still present mainly due to clouds coverage disturbs. Moreover, we have not confined our analysis to MSI SENTINEL-2 SWIR and MODIS-MIROVA MIR trend comparisons, but we explored what these correlations, between number of hot pixels detected (hot area) and the heat flux, could decrypt about thermal features related to different volcanic processes.
Significantly, the algorithm presented here is, as far we know, the first SENTINEL-2 multispectral based volcanic thermal detection process that runs operationally. In fact, the algorithm is part of the multiplatform MOUNTS volcanic monitoring system online since the beginning of 2018 (http://www.mounts-project.com/home; [30]). This system uses the SENTINEL constellation (-1, -2 and -5P) to retrieve and display key parameters volcano monitoring timeseries, such as deformation, heat anomalies and SO 2 flux, in near real time.
In this work, we demonstrated, through an almost four-year-long investigation of eight case studies, how SENTINEL-2 thermal signal analysis can enhance the study and monitoring of several volcanic processes, including lava-flow morphometric evaluations, extrusion phases and growing dynamics of lava domes, lava-lake pulsations, fumarolic activity, multiple active craters thermal activity, thermal precursors and new hot spots' presence, or periods of magma column rising in open-vent volcanoes.
A future perspective for monitoring purposes is to build an integration of the two thermal datasets, joining both the high-spatial-resolution potentialities of MSI SENTINEL-2 and high temporal resolution of MODIS data processed by the MIROVA system, in order to provide a specifically devoted product for the volcanic thermal activity characterization. In Figure 15, an example is provided with an application to Popocatépetl volcano (Mexico, [86]). The two SENTINEL-2 RGB images allow to visualize the thermal anomaly location and status into the crater, with a double zoom visualization 2 × 2 km and 10 × 10 km and immediate estimation of the number of hot pixels and maximum distance from the summit of volcano, while the plots on short (two months) and long (two years) timespan allow to evaluate the thermal activity and to track the evolution of hot area exposed and heat flux radiated. The two RGB S2 images with a different detail level, the number hot pixels and the heat flux produced by volcanic activity, summarized in a near real-time produced output, can provide an easy-to-understand product, rich of relevant qualitative (location and presence of hot spot) and quantitative (heat flux, hot area exposed) information both on the past and current thermal activity state of volcanoes.
Remote Sens. 2020, 12, x FOR PEER REVIEW 28 of 33 A future perspective for monitoring purposes is to build an integration of the two thermal datasets, joining both the high-spatial-resolution potentialities of MSI SENTINEL-2 and high temporal resolution of MODIS data processed by the MIROVA system, in order to provide a specifically devoted product for the volcanic thermal activity characterization. In Figure 15, an example is provided with an application to Popocatépetl volcano (Mexico, [86]). The two SENTINEL-2 RGB images allow to visualize the thermal anomaly location and status into the crater, with a double zoom visualization 2 × 2 km and 10 × 10 km and immediate estimation of the number of hot pixels and maximum distance from the summit of volcano, while the plots on short (two months) and long (two years) timespan allow to evaluate the thermal activity and to track the evolution of hot area exposed and heat flux radiated. The two RGB S2 images with a different detail level, the number hot pixels and the heat flux produced by volcanic activity, summarized in a near real-time produced output, can provide an easy-to-understand product, rich of relevant qualitative (location and presence of hot spot) and quantitative (heat flux, hot area exposed) information both on the past and current thermal activity state of volcanoes.  11 and 8a). On the top-left, a 2 × 2 km zoomed-in view over crater area; on the bottom-left, a 10 × 10 granule with an estimation of the number of hot pixels and the maximum distance of hot pixels from the volcano summit [86]. On the right, two thermal y-logarithmic timeseries, two months and two years on the top and bottom respectively, with MIROVA thermal flux in Watt, blue stems, and S2Pix in red dots.
Considering the SENTINEL-2 global coverage, and the possible improvement by integrating the LANDSAT 8 OLI data (with 30 m/pixel resolution in the NIR/SWIR bands), a similar output of that exposed in Figure 15, could be produced over the most active volcanoes with a global scale and with a revisit time of a few days. These outputs could be made available to observatories, monitoring centers and local authorities, in order to improve already existing monitoring systems and to contribute to the thermal surveillance of volcanic activity, as already partially suggested by some recent attempts [20,87].  11 and 8a). On the top-left, a 2 × 2 km zoomed-in view over crater area; on the bottom-left, a 10 × 10 granule with an estimation of the number of hot pixels and the maximum distance of hot pixels from the volcano summit [86]. On the right, two thermal y-logarithmic timeseries, two months and two years on the top and bottom respectively, with MIROVA thermal flux in Watt, blue stems, and S2Pix in red dots.
Considering the SENTINEL-2 global coverage, and the possible improvement by integrating the LANDSAT 8 OLI data (with 30 m/pixel resolution in the NIR/SWIR bands), a similar output of that exposed in Figure 15, could be produced over the most active volcanoes with a global scale and with a revisit time of a few days. These outputs could be made available to observatories, monitoring centers and local authorities, in order to improve already existing monitoring systems and to contribute to the thermal surveillance of volcanic activity, as already partially suggested by some recent attempts [20,87].
The current availability of satellite sensors with InfraRed detection capabilities in diverse wavelengths (MIR, TIR and SWIR) gives the possibility to characterize and monitor volcanic activity with an unprecedent level of detail. The MSI SENTINEL-2 multichannel thermal detection algorithm proposed here, even if represents a first-stage effort, and its correlation with MODIS-MIROVA analysis, fits into this view and aims to contribute to the monitoring and understanding of volcanic activity on a global scale.