A Multi-Channel Algorithm for Mapping Volcanic Thermal Anomalies by Means of Sentinel-2 MSI and Landsat-8 OLI Data

The Multispectral Instrument (MSI) and the Operational Land Imager (OLI), respectively onboard Sentinel-2A/2B and Landsat 8 satellites, thanks to their features especially in terms of spatial/spectral resolution, represents two important instruments for investigating thermal volcanic activity from space. In this study, we used data from those sensors to test an original multichannel algorithm, which aims at mapping volcanic thermal anomalies at a global scale. The algorithm, named Normalized Hotspot Indices (NHI), combines two normalized indices, analyzing near infrared (NIR) and short wave infrared (SWIR) radiances, to identify hotspot pixels in daylight conditions. Results, achieved studying a number of active volcanoes located in different geographic areas and characterized by a different eruptive behavior, demonstrated the NHI capacity in mapping both subtle and more intense volcanic thermal anomalies despite some limitations (e.g., missed detections because of clouds/volcanic plumes). In addition, the study shows that the performance of NHI might be further increased using some additional spectral/spatial tests, in view of a possible usage of this algorithm within a known multi-temporal scheme of satellite data analysis. The low processing times and the straight forth exportability to data from other sensors make NHI, which is sensitive even to other high temperature sources, suited for mapping hot volcanic targets integrating information provided by current and well-established satellite-based volcanoes monitoring systems.


Introduction
Several studies have shown the important role of satellite observations in monitoring thermal volcanic activity, as a complement to field measurements or as a unique source of information, particularly in remote and inaccessible areas, where traditional surveillance systems often lack (e.g., [1][2][3][4][5][6][7][8][9]).
In this work, we present a single image multichannel algorithm for mapping volcanic thermal anomalies at a global scale, by means of infrared Sentinel 2 MSI and Landsat 8 OLI data. To test the algorithm, christened normalized hotspot indices (NHI), we investigated eight active volcanoes located in different geographic areas ( Figure 1) and characterized by a different eruptive behavior. The proposed method, whose detections are assessed here by means of independent ground and satellite-based observations, uses two normalized indices to map volcanic thermal anomalies in daylight conditions analyzing SWIR and NIR radiances. Advantages and limitations in using these indices, which are different from other ones as the normalized thermal Index [16] (NTI) and thermal anomaly index [37] (TAI), respectively applied to MODIS and ASTER data, are analyzed and discussed in the next sections.

Mt. Etna Eruptive Events of July 2018-February 2019
Mt. Etna (37.748 • N; 14.999 • E) is the largest volcano in Europe and one of most active volcanoes in the world. Its eruptive activity frequently occurs from summit craters, although some intense flank eruptions took place in the last decades (e.g., [38]).
Remote Sens. 2019, 11, 2876 3 of 25 In July 2018, after several months of low-level activity, an increase of Strombolian activity was recorded at Bocca Nuova (BN) (from BN-1 and BN-2 vents) and at North-East Crater (NEC). In late August, the E vent of the New South-East Crater (NSEC) complex emitted new lava flows. In the following three months, low-intensity Strombolian activity and intermittent ash emissions occurred from multiple vents at summit craters. In late November, a lava flow emerged from a small scoria cone inside the NSEC, continuing until the end of the month [39]. Between 4 and 24 December, Strombolian activity produced a fan-shaped lava field at NSEC [40].  December, a flank eruption occurred. In more detail, in the early morning of 24 December, a magmatic dyke intruded on the southeast flank of the volcano, leading to an intense seismic swarm culminating, two days after, with the strongest event (ML 4.8). This phenomenon caused extensive ground deformations along the southeastern flank of Mt. Etna, and superficial faulting having seriously damaged numerous roads and buildings [40]. At around 11:00 UTC, a main NNW-SSE eruptive fissure, 2 km in length, opened on the eastern flank of the volcano, at an altitude of 3100-2400 m above sea level (a.s.l.). Contemporary, another small fissure (~70 m) opened, between NEC and BN, at around 3000 m a.s.l., showing a weak Strombolian activity lasting a few tens of minutes only. The main fissure feeds some lava flows, extending up to about 1650 m a.s.l., inside the deserted Valle del Bove. A few hours before, and during the eruption, both BN and NEC produced a Strombolian activity, which was accompanied by the emission of dense ash columns. Since the beginning of January 2019, both NEC and BN vents emitted ash, while degassing and discontinuous explosive activity took place in February [41].

Eruptive Events of Stromboli (Italy) Volcano of March-June 2019
Stromboli (38.789 • N; 15.213 • E) is a stratovolcano located in the Aeolian Island (Italy). Strombolian eruptions and minor effusive activities generally characterize this volcano, although some important lava effusions occurred in recent years (e.g., [42][43][44]), During March-April 2018, a number of explosions occurred from active vents located at the summit crater area. Eruptive activity continued through the following two months. In January 2019, magma remained high inside the conduits, and intermittent Strombolian explosions occurred from multiple vents. Afterwards, volcanic activity gradually decreased in intensity returning to normal levels. Small to moderate Strombolian explosions then occurred, at irregular intervals, from several vents through March 2019 [45]. Since 12 June, volcanic activity increased remaining at a moderate level until one day before the strong paroxysm of 3 July, when a gas bubble arrived at the surface, generating two lateral blasts causing the emissions of lava fragments at high temperature, and injecting a number of bush fires along the volcano slopes [46].

Erta Ale (Africa) Eruptions of 2017
Erta Ale (13.6 • N; 40.67 • E) is a shield volcano located in the Afar region of northeastern Ethiopia (Africa). It is among the few volcanoes in the world showing a persistent lava lake, which has been active during most of the past decades [47]. Two active craters (i.e., Northern and Southern ones), within the larger oval-shaped summit caldera, exhibit periodic lava fountaining causing lava-lake overflows [48].
Since late December 2016, an effusive eruption took place from active craters. On the morning of 16 January, 2017, the lake overflowed the W rim of the crater. Five days later, a new fissure opened on the SE volcano flank (about 4 km from the caldera) emitting lava [48,49]. Activity at the fissure vent strongly increased during mid-April/early June 2017, when pahoehoe lava flows moved in both NE and SW direction [48].

Kilauea (Hawaii) Eruptive Activity of March 2018
Kilauea (19.421 • N; 155.287 • W) is a shield volcano located in the Hawaiian Islands (United States), which has been continuously erupting since 1983 when open lava lakes, and flows from the summit caldera and East Rift Zone, occurred (e.g., [50,51]). Eruptive activity generally occurs from two vents, the Halema'uma'u crater (within the caldera), and the Pu'u 'Ō'ō cone, which is located in the East Rift Zone at about 20 km from the summit (e.g., [52]).
On 26 March, 2018, after a brief swarm of small earthquakes recorded in the upper East Rift Zone, and a long-period earthquake swarm, a small lava flow began erupting onto the Pu'u 'O'o crater floor. Eruptive activity followed a small increase in seismic events at Pu'u 'O'o, and the opening of a first vent, forming a small pit about 35 m wide, five days before [52].

Sakurajima (Japan) Eruptive Events of September-October 2018
Sakurajima (31.593 • N; 130.657 • E) is an active volcano located in the Aira caldera, in southern Kyushu (Japan), which usually emits ash plumes and scatters blocks during eruptions from craters.
In July-December 2018, several ash emissions took place at the Minamidake crater. In September, when the maximum plume height was estimated around 2.3 km above the crater, incandescent blocks were emitted on volcano flanks. At the end of October, an observation flight revealed the occurrence of a fumarolic activity at the Showa crater. A dilute brown plume was observed at the Minamidake, which also showed an intermittent incandescence visible at the nighttime [53].

Popocatepetl (Mexico) Eruptive Activity of March 2019
Popocatépetl (19.023 • N; 98.622 • W) is a stratovolcano located in Central Mexico, which has been persistently active since 2005. Eruptive activity generally consists in the emission of gas/steam and ash plumes. Moreover, incandescent blocks are often scattered across the volcano flanks, and frequent growth of lava domes occur in the summit crater [54].
During March-August 2018, large emissions of volcanic SO 2 occurred; steam/gas emissions were then observed in the following months. On 3 March, 2019, the National Center for Prevention of Disasters (CENAPRED) reported the emission of ash and water vapor. Volcanic explosions prompted local authorities to issue a yellow alert (middle level on a three-color scale), implementing a security perimeter in a 12-km radius of the crater [55]. On 14 March, a dense ash plume drifting NNE was emitted [56]. On 26 March, incandescent fragments ignited fires on volcano flanks. Three days later, an eruption occurred, generating an ash plume moving in the SE direction [56,57].

Shiveluch (Kamchatka) Activity of March 2016-April 2019
Shiveluch is one of the most active andesitic volcanoes of Kamchatka (Russia), which generally shows Vulcanian explosive eruptions and periods of dome growth [58].
During March 2016-April 2017, Kamchatka Volcanic Eruption Response Team (KVERT), performing daily video-visual and satellite monitoring of volcanoes in the Russian Federation [59], reported a lava-dome extrusion onto the N flank, strong fumarolic activity, dome incandescence, ash explosions, and hot avalanches. Thermal anomalies were then observed by satellite from August 2017 through January 2018 [60]. Between December 2018 and April 2019, both KVERT and Middle InfraRed Observation of Volcanic Activity (MIROVA), which monitors active volcanoes in near-real time by means of infrared MODIS data [61], revealed an increase in thermal volcanic activity. In early January and through April 2019, lava dome continued to grow extruding viscous lava blocks. Additionally, strong fumarolic emissions, gas-and-steam plumes containing ash were reported [62].

Piton de la Fournaise (Reunion Islands) Eruptions of April 2018-June 2019
Piton de la Fournaise is a basaltic shield volcano located in the Réunion Island (Indian Ocean), whose eruptive activity generally consist of lava fountains and effusive eruptions producing large lava flows (e.g., [63]). In the evening of 27 April, 2018 a seismic swarm accompanied by rapid deformation occurred, indicating the magma migration towards the surface. A few hours later, a new flank eruption started from fissures located at Rivals crater and SW flank of the Dolomieu crater. Eruptive activity continued through May, when lava flows were generally confined to tubes, and a small area of incandescence affected the main crater [64]. In 2019, two important effusive events occurred. The first one, lasting from 18 February-10 March, 2019, was associated with the opening of multiple fissures on the E flank of the Dolomieu crater. The second one began on 11 June at the SSE outer slope of Dolomieu crater, and lasted about 48 h [65].

Sentinel-2 MSI Data
The Sentinel-2 mission was established within the European Space Agency (ESA) Copernicus program to perform multispectral high spatial resolution optical observations of land surfaces, inland and coastal waters (e.g., [66]). The mission is a constellation of two polar satellites (i.e., Sentinel-2A/2B), placed in the same sun-synchronous orbit and phased at 180 • to each other, orbiting the Earth at 786 km altitude [67]. The image swath width is about 290 km, which is larger than Landsat [68]. The Sentinel-2 constellation guarantees a revisit time of 5 days at the equator in cloud-free conditions, increasing up to 3-5 days at mid-latitudes [66]. Sentinel-2 satellites are equipped with the MSI instrument, which is a state of the art push-broom imager, measuring the solar radiance reflected from the Earth in 13 spectral bands from the visible (VIS) to SWIR region (i.e., 443-2190 nm), with a spatial resolution of 10-60 m (see Table 1) [67,68]. Sentinel-2 data, made available to users, include the Level-1C and Level-2A products distributed through the ESA Open Access Hub [69]. The Sentinel-2 L1C product is a calibrated top-of-atmosphere (TOA) reflectance product composed of 100 × 100 km 2 tiles in the Universal Transverse Mercator/World Geodetic System 1984 (UTM/WGS84) projection. The product, containing parameters to transform reflectance into radiances, is resampled with a constant ground sampling distance, based on native resolution of the MSI spectral bands [70,71]. The Level 2A product is a bottom-of-atmosphere (BOA) reflectance product, which is in the same projection of the associated Level 1C.

Landsat 8 OLI Data
Since 1972, the joint National Aeronautics and Space Administration (NASA)/U.S. Geological Survey (USGS) Landsat series have continuously acquired imagery of the Earth's land surface [72]. The Landsat 8 satellite, which was launched in 2013, orbits the Earth in a sun-synchronous, near-polar orbit at an altitude of 705 km. The revisit time is 16-days and the image swath is about 185 km. Thermal Infrared Sensor (TIRS) and OLI are the two instruments onboard this satellite.
TIRS is a two-band thermal sensor operating from 10.6-12.5 µm. OLI is a pushbroom sensor, having a 12-bit radiometric resolution, providing data in nine spectral channels, with a spatial resolution from 15 m (panchromatic) to 30 m in the visible/infrared bands (see Table 2) (e.g., [73]).
In this work, we used the Landsat 8 OLI Level 1 data products, made freely available online by USGS, for running the NHI algorithm. We converted digital numbers (DNs) to TOA spectral radiance by using the rescaling coefficients found in the metadata (MTL) file.

Background
Several algorithms have been, up to now, proposed to detect and map thermal anomalies (e.g., fires; lava flows) using high-spatial resolution satellite data. Among the most recent ones, HOTMAP [30] uses two binary parameters, analyzing the TOA reflectance measured in bands 5 (0.87 µm), 6 (1.61 µm), and 7 (2.2 µm) of the OLI sensor, to detect hotspots. The algorithm was compared to Autonomous Sciencecraft Experiment (ASE) developed for Hyperion [28] and to a contextual method originally designed for ASTER [74], showing better performance in detecting thermal anomalies at a global scale [30]. Other authors proposed an enhancement of HOTMAP, coupling a contextual analysis to the spectral ratios of SWIR signal, for better investigating volcanic thermal anomalies by means of Sentinel 2 MSI data [75]. The method runs operationally within the Monitoring Unrest from Space (MOUNT) platform, which monitors several active volcanoes in near-real time through a multi-sensor observing system [34]. These techniques, such as other novels algorithms of fire detection (e.g., those integrating a multi-temporal scheme of satellite data analysis [76] or using a three-band index analyzing both NIR and SWIR radiances [77]), confirmed the important role of high-spatial resolution satellite observations in supporting volcanic hazard and fire management activities. In this context, the algorithm presented in this work, using a simple detection scheme, which may be easily applied to data from different sensors, aims at performing a quick and efficient mapping of volcanic thermal anomalies in daylight conditions, regardless of geographic areas and periods of the year.

NHI Algorithm
Hot targets (e.g., lava flows) are particularly radiant in the SWIR spectral region (e.g., [5]). Moreover, higher temperature surfaces emit strongly at shorter rather than longer SWIR wavelengths (e.g., [78]). Hence, sensors having channels in the SWIR bands can be used to identify and map volcanic thermal anomalies (e.g., [78,79]). On the other hand, during daytime the reflected solar radiation represents a major issue (e.g., [80,81]). The NHI algorithm takes into account that mentioned above using two normalized indices to detect and map volcanic hotspots: where, L 2.2 , L 1.6 , and L 0.8 are the TOA radiances [W m −2 sr −1 µm −1 ] measured, for each pixel of the analyzed scene, at around 2.2 µm, 1.6 µm (SWIR), and 0.8 µm (NIR) wavelengths, in the relative MSI/OLI spectral bands. The index in Equation (1), analyzing SWIR radiances, takes into account that previously shown by other authors using airborne data [82]. The index in Equation (2), which uses the same (1) where, 2.2 , 1.6 , and 0.8 are the TOA radiances [W m −2 sr −1 µm −1 ] measured, for each pixel of the analyzed scene, at around 2.2 µ m, 1.6 µ m (SWIR), and 0.8 µ m (NIR) wavelengths, in the relative MSI/OLI spectral bands. The index in Equation (1), analyzing SWIR radiances, takes into account that previously shown by other authors using airborne data [82]. The index in Equation (2), which uses the same spectral bands of the Normalized Burnt Ratio (NBR) highlighting burn surfaces [83], should enable a more efficient identification of high-temperature surfaces. We assessed the NHI indices behavior by analyzing several Sentinel 2 MSI and Landsat 8 OLI scenes. An example is shown in Figure  The false color composite imagery, which was generated by combining bands 12(Red), 11(Green), and 8A(Blue) of MSI after converting TOA reflectance to radiance by means of the Sentinel Application Platform (SNAP) tool from ESA, emphasized hot features, clouds, background surfaces, and a volcanic plume, dispersing in the SW direction, in different colors. In particular, hotspot pixels, appearing yellow/red colored, could be well discriminated from other features as clouds (pixels from gray to white). Indeed, radiance changes recorded along the transect regions of Figure 2a showed The false color composite imagery, which was generated by combining bands 12(Red), 11(Green), and 8A(Blue) of MSI after converting TOA reflectance to radiance by means of the Sentinel Application Platform (SNAP) tool from ESA, emphasized hot features, clouds, background surfaces, and a volcanic plume, dispersing in the SW direction, in different colors. In particular, hotspot pixels, appearing yellow/red colored, could be well discriminated from other features as clouds (pixels from gray to white). Indeed, radiance changes recorded along the transect regions of Figure 2a showed that clouds and snow-covered areas leaded to the peak of L 0.8 radiance, highest L 1.6 values occurred over clouds (because of reflected solar radiance) and active summit vents (due to emitted radiance), while L 2.2 radiance increased over both crater area and the SE volcano flank (see Figure 2b).
By analyzing the behavior of NHI indices along the same transect regions of Figure 2a, we found that volcanic thermal anomalies could be well distinguished from other features as evident in Figure 3, showing the NHI imagery (top panels) and the spatial changes of NHI indices (middle-bottom panels).  Looking at the plots, it can be noted as while the index in Equation (1) increased up to positive values at crater borders and on SE flank of the volcano (Figure 3a), where lava effusion was significantly less intense than previous days, the index in Equation (2) became positive mainly over summit active vents (Figure 3b). This analysis, such as similar investigations, revealed a sort of complementary behavior of the two indices, suggesting to us to analyze pixels with values of > 0 or > 0 to identify and map volcanic thermal anomalies in the areas of Figure  1. The index in Equation (1) seems to be, in fact, more sensitive to low-medium intensity hotspots but less capable of mapping strong thermal anomalies than the index in Equation (2).  Looking at the plots, it can be noted as while the index in Equation (1) increased up to positive values at crater borders and on SE flank of the volcano (Figure 3a), where lava effusion was significantly less intense than previous days, the index in Equation (2) became positive mainly over summit active vents ( Figure 3b). This analysis, such as similar investigations, revealed a sort of complementary behavior of the two indices, suggesting to us to analyze pixels with values of NHI SWIR >0 or NHI SWNIR > 0 to identify and map volcanic thermal anomalies in the areas of Figure 1. The index in Equation (1) seems to be, in fact, more sensitive to low-medium intensity hotspots but less capable of mapping strong thermal anomalies than the index in Equation (2).

Results
In this section, we show results achieved analyzing hundreds of Sentinel 2 MSI and Landsat 8 OLI scenes by means of NHI algorithm; investigations were extended also to some neighbor non-volcanic areas of Figure 1, for better analyzing factors affecting NHI performance.
NHI maps shown here display hotspot pixels in two different colors. In detail, we depicted hotspot pixels in red (NHI SWNIR > 0) and yellow (NHI SWIR > 0), retrieving qualitative information about areas where thermal emissions were strong or less intense (see previous section).

Mt. Etna Eruptive Events of July 2018-February 2019
To assess the NHI potential in mapping thermal anomalies at Mt. Etna (Sicily, Italy) area, we processed all the Sentinel-2 MSI scenes acquired during July 2018-February 2019 at around 10:00 UTC.   As can be seen from the figure, NHI detected a thermal anomaly at the crater area since 9 March, when the cloud coverage was less significant. The thermal anomaly, including two clusters of hotspot During July-August 2018, NHI flagged several hotspot pixels at the crater area, and more in particular at BN-1 and BN-2, which were located on the bottom of Bocca Nuova (BN) crater, and at NEC, because of ongoing Strombolian activity. In addition, NHI correctly identified the lava effusion occurring at the end of August from the Bocca della Sella vent, which is located between the old and the New South-East Crater (NSEC) (see panels of 26-28 August). NHI maps of September-October indicated that both BN-1 and BN-2 were active, such as the NEC. Additionally, they showed the presence of a thermal anomaly at NSEC (see panel of 30 October), which was possibly associated to the fumarolic activity. Since the second week of November, the number of hotspot pixels (mostly those depicted in red) increased, revealing the presence of a more intense and extended thermal anomaly. In mid-December, NHI maps showed the occurrence of a thermal activity at NEC, BN, and NSEC, while those of 24-26 December marked the lava effusion occurring from a new fissure which opened in the upper part of Valle del Bove. Since the end of December, the number of hotspot pixels significantly decreased, revealing the intensity reduction of thermal volcanic activity. In particular, NHI flagged some hot spots pixels mainly at BN and NSEC (e.g., see panels of 2 and 27 February 2019), where a degassing activity was independently reported (see Section 2.1).
NHI detections described above were in good agreement with independent ground-based observations. On the other hand, although NHI provided reliable information about active vents, the manual inspection of satellite imagery (both visible and infrared ones) and the comparison with the RGB product (generated as Figure 2) revealed that in a few maps some background pixels were erroneously flagged as hotspots (e.g., as for panel of 26 December 2018). Pros and cons of the here proposed algorithm, emphasized also by results achieved using NHI for mapping thermal anomalies in the other volcanic areas of Figure 1, are analyzed and discussed in detail in Section 6.  As can be seen from the figure, NHI detected a thermal anomaly at the crater area since 9 March, when the cloud coverage was less significant. The thermal anomaly, including two clusters of hotspot pixels (see region marked in blue magnified at the bottom right side of the figure), appeared more extended and intense on the SW side of the summit crater area (e.g., see panel of 24 March), where presumably higher was the number of active vents. By analyzing the Sentinel-2 MSI scenes of June, we observed a general increase of thermal volcanic activity starting from the second week of the month (e.g., see panel of 17 June). Strombolian explosions, as well as a degassing activity occurring from multiple vents within the crater, corroborated information retrieved by satellite.

Activity of Stromboli Vlcano of March-June 2019
It is interesting to note that NHI identified a volcanic thermal anomaly, at the crater area, even one day before the strong paroxysm of 3 July (see Figure 5). This powerful event was the largest occurring since the 1930s, and caused the death of a hiker [46]. However, the detected thermal anomaly, which was partially masked by clouds (see panel of 2 July), did not represent a thermal precursor of that eruption. Indeed, the moderate Strombolian activity occurring since mid-June and continuing through 2 July [84] generated the thermal anomaly identified and mapped by NHI. In February, other clusters of hotspot pixels were recognizable, revealing since the end of same month a general increase of effusive activity. Afterwards, meteorological clouds mostly obscured the ROI, leading to the underestimation of volcanic thermal anomaly (e.g., see panel of 20 March).

Erta-Ale Euptions of 2017
We assessed the NHI potential in mapping thermal anomalies at Erta Ale also in other periods of the year. As an example, Figure 7 displays two NHI maps of December 2017 highlighting the presence of a thermal anomaly at the caldera, where the lava lake was persistent (Figure 7a), as well as East of analyzed satellite scenes (e.g., see region marked in green at the top-right side of Figure  7b).
Therefore, based on information provided by NHI, the lava flow emitted on SE flank of the volcano since January 2017 moved even in the NE direction fitting with that reported by volcanological bulletins (see Section 2.3.). In February, other clusters of hotspot pixels were recognizable, revealing since the end of same month a general increase of effusive activity. Afterwards, meteorological clouds mostly obscured the ROI, leading to the underestimation of volcanic thermal anomaly (e.g., see panel of 20 March).
We assessed the NHI potential in mapping thermal anomalies at Erta Ale also in other periods of the year. As an example, Figure 7 displays two NHI maps of December 2017 highlighting the presence of a thermal anomaly at the caldera, where the lava lake was persistent (Figure 7a), as well as East of analyzed satellite scenes (e.g., see region marked in green at the top-right side of Figure 7b).
Therefore, based on information provided by NHI, the lava flow emitted on SE flank of the volcano since January 2017 moved even in the NE direction fitting with that reported by volcanological bulletins (see Section 2.

Kilauea Eruptive Activity of March 2018
In Figure 8, we show two NHI maps generated from Sentinel-2 MSI data of March 2018, including the area of Kilauea (Hawaii Islands) volcano. Specifically, Figure 8a displays the NHI map of 4 March showing the presence of a thermal anomaly at Halema 'uma'u crater (see yellow pixels within the region marked in green), which was partially identified by satellite due to a thick cloud coverage and a volcanic plume dispersing in the SW direction. Map of Figure 8b, which was generated from satellite data of 29 March, shows that NHI identified the volcanic thermal anomaly in a more efficient way when the aforementioned features did not affect the crater area. In addition, the figure shows that NHI flagged several hotspot pixels at the East Rift Zone (see red/yellow pixels within the blue box magnified at the bottom right side of the figure), which were ascribable to the eruptive activity began in March 2018 at the Pu'u O'o cone (see Section 2.4). Figure 9 shows a sequence of NHI maps, generated from Sentinel 2 MSI data of October 2018, covering the area of Sakurajima volcano. Maps show that since the third week of October (i.e., when the cloud coverage was less significant), NHI flagged the occurrence of a moderate thermal activity at the Minamidake crater (see yellow pixels within the region marked in blue, magnified at the top

Kilauea Eruptive Activity of March 2018
In Figure 8, we show two NHI maps generated from Sentinel-2 MSI data of March 2018, including the area of Kilauea (Hawaii Islands) volcano. Specifically, Figure 8a displays the NHI map of 4 March showing the presence of a thermal anomaly at Halema 'uma'u crater (see yellow pixels within the region marked in green), which was partially identified by satellite due to a thick cloud coverage and a volcanic plume dispersing in the SW direction.

Kilauea Eruptive Activity of March 2018
In Figure 8, we show two NHI maps generated from Sentinel-2 MSI data of March 2018, including the area of Kilauea (Hawaii Islands) volcano. Specifically, Figure 8a displays the NHI map of 4 March showing the presence of a thermal anomaly at Halema 'uma'u crater (see yellow pixels within the region marked in green), which was partially identified by satellite due to a thick cloud coverage and a volcanic plume dispersing in the SW direction. Map of Figure 8b, which was generated from satellite data of 29 March, shows that NHI identified the volcanic thermal anomaly in a more efficient way when the aforementioned features did not affect the crater area. In addition, the figure shows that NHI flagged several hotspot pixels at the East Rift Zone (see red/yellow pixels within the blue box magnified at the bottom right side of the figure), which were ascribable to the eruptive activity began in March 2018 at the Pu'u O'o cone (see Section 2.4). Figure 9 shows a sequence of NHI maps, generated from Sentinel 2 MSI data of October 2018, covering the area of Sakurajima volcano. Maps show that since the third week of October (i.e., when the cloud coverage was less significant), NHI flagged the occurrence of a moderate thermal activity at the Minamidake crater (see yellow pixels within the region marked in blue, magnified at the top  Figure 8b, which was generated from satellite data of 29 March, shows that NHI identified the volcanic thermal anomaly in a more efficient way when the aforementioned features did not affect the crater area. In addition, the figure shows that NHI flagged several hotspot pixels at the East Rift Zone (see red/yellow pixels within the blue box magnified at the bottom right side of the figure), which were ascribable to the eruptive activity began in March 2018 at the Pu'u O'o cone (see Section 2.4). Figure 9 shows a sequence of NHI maps, generated from Sentinel 2 MSI data of October 2018, covering the area of Sakurajima volcano. Maps show that since the third week of October (i.e., when the cloud coverage was less significant), NHI flagged the occurrence of a moderate thermal activity at the Minamidake crater (see yellow pixels within the region marked in blue, magnified at the top right side of the figure). NHI did not detect, however, any hotspot pixel at the Showa crater where independent observations reported the occurrence of a fumarolic activity (see Section 2.5). The manual inspection of satellite imagery, and the comparison with the qualitative RGB product, performed as for Mt. Etna (see Section 5.1), corroborated that indicated by NHI. Thus, it is realistic to suppose that the fumarolic activity was not in progress at the time of analyzed satellite observations or it was particularly weak and then difficult to be identified by satellite.

Sakurajima Thermal Activity of October 2018
right side of the figure). NHI did not detect, however, any hotspot pixel at the Showa crater where independent observations reported the occurrence of a fumarolic activity (see Section 2.5). The manual inspection of satellite imagery, and the comparison with the qualitative RGB product, performed as for Mt. Etna (see Section 5.1), corroborated that indicated by NHI. Thus, it is realistic to suppose that the fumarolic activity was not in progress at the time of analyzed satellite observations or it was particularly weak and then difficult to be identified by satellite.
It is worth noting that by analyzing the entire satellite sub-scenes, we found that NHI flagged some hotspot pixels in areas located far from the eruptive center (e.g., see region marked in green on panel of 29 October). These pixels represented artefacts generated by possible cloud-edges effects, which should be taken into account when the algorithm is used. Specifically, maps show that a volcanic plume affected most of analyzed satellite scenes. As a result of this issue, NHI partially identified a thermal anomaly within the summit crater in early March (e.g., see panel of 3 March), while it did not flag any hotspot pixel ascribable to thermal volcanic activity at end of same month (see panel of 28 March). Nevertheless, NHI maps showed the presence of some thermal anomalies in the neighbor non-volcanic areas. We investigated these features using the same approach of Section 5.1. Based on those analyses, we found that thermal anomalies flagged outside the volcano edifice did not represent false positives ascribable to clouds (e.g., Section 5.5) or background pixels (Section 5.1), but were generated by other high temperature sources. In particular, possible fires, occurring both at the slopes (e.g., see pixels within the green box) and in the surrounding of the investigated volcanic area (green circles), generated the thermal anomalies identified and mapped by satellite. It is worth noting that by analyzing the entire satellite sub-scenes, we found that NHI flagged some hotspot pixels in areas located far from the eruptive center (e.g., see region marked in green on panel of 29 October). These pixels represented artefacts generated by possible cloud-edges effects, which should be taken into account when the algorithm is used. Specifically, maps show that a volcanic plume affected most of analyzed satellite scenes. As a result of this issue, NHI partially identified a thermal anomaly within the summit crater in early March (e.g., see panel of 3 March), while it did not flag any hotspot pixel ascribable to thermal volcanic activity at end of same month (see panel of 28 March). Nevertheless, NHI maps showed the presence of some thermal anomalies in the neighbor non-volcanic areas. We investigated these features using the same approach of Section 5.1. Based on those analyses, we found that thermal anomalies flagged outside the volcano edifice did not represent false positives ascribable to clouds (e.g., Section 5.5) or background pixels (Section 5.1), but were generated by other high temperature sources. In particular, possible fires, occurring both at the slopes (e.g., see pixels within the green box) and in the surrounding of the investigated volcanic area (green circles), generated the thermal anomalies identified and mapped by satellite. Figure 10. NHI maps from Sentinel-2 MSI data of March 2019 covering the Popocatepetl. Yellow/red pixels, within the region marked in blue, indicated an intra-crater thermal activity. Hotspot pixels within the areas marked in green were possibly associated to fires occurring at the volcano's slopes (green box) and far from the volcano edifice (green circles).

Shiveluch Thermal Activity of March 2016-April 2019
To assess the NHI capacity in mapping volcanic thermal anomalies also at the high latitude regions (i.e., in presence of a cold background), we investigated the Shiveluch thermal activity by means of Landsat 8 OLI data of March 2016April 2019.
Results of this analysis are shown in Figure 11a, revealing the occurrence of three major events (i.e., 2 March, 2016, 13 September, 2017 and 7 February, 2019) at the target area, and of a number of minor thermal volcanic activities, during the period of interest. Both major and minor events were consistent with information provided by the RGB product, which was generated using bands 7 (Red), 6 (Green), and 5 (Blue) of OLI. Figure 11b displays, for instance, the NHI map from Landsat 8 OLI data of 2 March 2016, showing the good agreement of detected thermal anomaly (see yellow pixels on bottom-left panel) with that emphasized by the false color composite imagery (see yellow/red colored pixels on bottom-right panel). Moreover, although NHI did not detect all the hotspot pixels associated to the eruptive activity in progress on that day probably due to the presence of a volcanic plume, the thermal anomaly in Figure 11a was consistent also with information provided by some well-established satellite-based monitoring systems (see Section 2.7). This agreement characterized also the thermal anomaly shown in Figure 11b, which was identified by NHI analyzing the Landsat 8 OLI data of 19 April, 2019. In particular, the detected thermal anomaly fits with the frequent and independent observation of volcanic hotspots, mostly located within the summit crater area, during the period January-April 2019 (e.g., [85]). volcanic plume Popocaté petl Figure 10. NHI maps from Sentinel-2 MSI data of March 2019 covering the Popocatepetl. Yellow/red pixels, within the region marked in blue, indicated an intra-crater thermal activity. Hotspot pixels within the areas marked in green were possibly associated to fires occurring at the volcano's slopes (green box) and far from the volcano edifice (green circles).

Shiveluch Thermal Activity of March 2016-April 2019
To assess the NHI capacity in mapping volcanic thermal anomalies also at the high latitude regions (i.e., in presence of a cold background), we investigated the Shiveluch thermal activity by means of Landsat 8 OLI data of March 2016-April 2019.
Results of this analysis are shown in Figure 11a, revealing the occurrence of three major events (i.e., 2 March, 2016, 13 September, 2017 and 7 February, 2019) at the target area, and of a number of minor thermal volcanic activities, during the period of interest. Both major and minor events were consistent with information provided by the RGB product, which was generated using bands 7 (Red), 6 (Green), and 5 (Blue) of OLI. Figure 11b displays, for instance, the NHI map from Landsat 8 OLI data of 2 March 2016, showing the good agreement of detected thermal anomaly (see yellow pixels on bottom-left panel) with that emphasized by the false color composite imagery (see yellow/red colored pixels on bottom-right panel). Moreover, although NHI did not detect all the hotspot pixels associated to the eruptive activity in progress on that day probably due to the presence of a volcanic plume, the thermal anomaly in Figure 11a was consistent also with information provided by some well-established satellite-based monitoring systems (see Section 2.7). This agreement characterized also the thermal anomaly shown in Figure 11b Figure 12 shows two NHI maps including the area of Piton de La Fournaise (Reunion Island), generated from Landsat 8 OLI data.

Piton de la Fournaise Lava Effusions of April 2018-June 2019
In detail, Figure 12a displays the NHI map of 12 May, 2018, where the detected thermal anomaly (see red pixels magnified at the bottom-right side of the figure) was consistent with a documented lava effusion began the previous month (see Section 2.8). Figure 12b shows the map of 24 February, 2019, providing information about the lava flow emitted on the E flank of the volcano. The comparison of those maps with the false color composite imagery corroborated information provided by NHI. Regarding the lava effusion of 11 June, 2019 (see Section 2.8), it was investigated using Sentinel 2 MSI data of that day, since the Landsat 8 OLI scene covering the target area was not available. Figure 12c displays the NHI map revealing the presence of a significant thermal anomaly (for intensity and spatial extent) affecting the SE flank of the volcano, which was mostly consistent with information provided by the RGB product (top-right panel). Main differences were ascribable to a few artefacts generated within the region marked in green (see yellow pixels), and to hotspot pixels undetected by NHI in the inner part of the lava flow (due to possible saturation effects of SWIR In detail, Figure 12a displays the NHI map of 12 May, 2018, where the detected thermal anomaly (see red pixels magnified at the bottom-right side of the figure) was consistent with a documented lava effusion began the previous month (see Section 2.8). Figure 12b shows the map of 24 February, 2019, providing information about the lava flow emitted on the E flank of the volcano. The comparison of those maps with the false color composite imagery corroborated information provided by NHI. Regarding the lava effusion of 11 June, 2019 (see Section 2.8), it was investigated using Sentinel 2 MSI data of that day, since the Landsat 8 OLI scene covering the target area was not available. Figure 12c displays the NHI map revealing the presence of a significant thermal anomaly (for intensity and spatial extent) affecting the SE flank of the volcano, which was mostly consistent with information provided by the RGB product (top-right panel). Main differences were ascribable to a few artefacts generated within the region marked in green (see yellow pixels), and to hotspot pixels undetected by NHI in the inner part of the lava flow (due to possible saturation effects of SWIR channels under investigation).
This example shows that NHI is a capable of well exploiting advantages of Landsat 8 OLI and Sentinel 2 MSI data integration despite some limitations.

Discussion
By testing the NHI algorithm in the volcanic areas of Figure 1, we found the good agreement of thermal anomalies identified using infrared Sentinel 2 MSI and Landsat 8 OLI data with volcanic activities reported by field reports. Additionally, independent satellite-based observations further corroborated information provided by NHI. In the case of Erta Ale eruptions, volcanic thermal anomalies were for instance consistent with those catalogued and collected by the ASTER Volcano Archive (AVA) system using ASTER data (e.g., see hotspot products in [86] temporally close to NHI maps of Figure 7). Concerning the Mt. Etna activity, NHI detections fitted with results of a recent independent study, performed using Sentinel 2 MSI and EOS-MODIS data [75] and with outputs of the highly sensitive RSTVOLC algorithm monitoring Italian volcanoes in near-real time (e.g., [87]). Figure 13 displays the temporal trend of hotspot pixels flagged by NHI ( Figure 13a) and RSTVOLC (Figure 13b) during the period July-February 2019. By looking at the figure, it can be noted as NHI detections were corroborated by the RSTVOLC ones (e.g., see thermal anomalies flagged during 2527 August and 2427 December, 2018 lava effusion periods). Additionally, although the low temporal sampling of Sentinel 2 MSI data did not enable the same continuity of NOAA/Metop (National Oceanic and Atmospheric Administration/Meteorological Operational Satellites)-AVHRR observations, the 20 m spatial resolution favored the identification of some minor events (e.g., early July 2018; February 2019). They included modest Strombolian activities generated by vents at the bottom of the summit craters, incandescent gas emissions, and fumaroles from fractures. This analysis confirmed advantages of integrating information from different satellite-based systems for

Discussion
By testing the NHI algorithm in the volcanic areas of Figure 1, we found the good agreement of thermal anomalies identified using infrared Sentinel 2 MSI and Landsat 8 OLI data with volcanic activities reported by field reports. Additionally, independent satellite-based observations further corroborated information provided by NHI. In the case of Erta Ale eruptions, volcanic thermal anomalies were for instance consistent with those catalogued and collected by the ASTER Volcano Archive (AVA) system using ASTER data (e.g., see hotspot products in [86] temporally close to NHI maps of Figure 7). Concerning the Mt. Etna activity, NHI detections fitted with results of a recent independent study, performed using Sentinel 2 MSI and EOS-MODIS data [75] and with outputs of the highly sensitive RST VOLC algorithm monitoring Italian volcanoes in near-real time (e.g., [87]). Figure 13 displays the temporal trend of hotspot pixels flagged by NHI (Figure 13a) and RST VOLC (Figure 13b) during the period July-February 2019. By looking at the figure, it can be noted as NHI detections were corroborated by the RST VOLC ones (e.g., see thermal anomalies flagged during 25-27 August and 24-27 December, 2018 lava effusion periods). Additionally, although the low temporal sampling of Sentinel 2 MSI data did not enable the same continuity of NOAA/Metop (National Oceanic and Atmospheric Administration/Meteorological Operational Satellites)-AVHRR observations, the 20 m spatial resolution favored the identification of some minor events (e.g., early July 2018; February 2019). They included modest Strombolian activities generated by vents at the bottom of the summit craters, incandescent gas emissions, and fumaroles from fractures. This analysis confirmed advantages of integrating information from different satellite-based systems for better investigating changes of thermal volcanic activity also in well-monitored volcanic areas (e.g., [88] better investigating changes of thermal volcanic activity also in well-monitored volcanic areas (e.g., [88]). Regarding the factors affecting the NHI performance, clouds had the major impact on results of this study. Indeed, about 38% of analyzed Sentinel 2 MSI scenes (about 90 in total) were almost completely overcast did not enabling the identification of volcanic thermal anomalies. Additionally, by the manual inspection of satellite imagery, and the comparison with RGB product, we found that NHI underestimated some volcanic thermal anomalies in presence of clouds/volcanic plumes (e.g., see Section 5.6) partially affecting the target areas and leading to negative values of used indices. On the other hand, background pixels showing radiance values slightly higher at 2.2 µ m than 1.6. µ m wavelength generated artefacts associated to positive values of the index in Equation (1) (e.g., see Section 5.1). To minimize this issue, an initial test on SWIR radiance might be used before computing the NHI indices as shown in Figure 14, displaying the thermal anomaly map generated, under Google Earth Engine [89] using the preliminary NHI tool, from Landsat OLI data of 24 February, 2019. It can be noted as artefacts in Figure 12c, i.e., false hotspot pixels within the region marked in green, could be efficiently removed optimizing the algorithm (e.g., filtering pixels with values of L2.2 < 3.0 W m −2 sr −1 µm −1 ). Hence, the confidence level of NHI detections may be further increased; in this direction, we could use also cloud-masking procedure to minimize false positives related to possible cloud-edge effects (see Section 5.5).
Since NHI is also sensitive to fires, a spatial distance filter might be integrated within the process to avoid the misinterpretation of detected thermal anomalies at the investigated volcanic areas. Bush fires generally develop along the volcano flanks due to the fallout of incandescent blocks on vegetated areas (e.g., see Section 2.6). Nevertheless, they may also occur because of anthropogenic causes, as Regarding the factors affecting the NHI performance, clouds had the major impact on results of this study. Indeed, about 38% of analyzed Sentinel 2 MSI scenes (about 90 in total) were almost completely overcast did not enabling the identification of volcanic thermal anomalies. Additionally, by the manual inspection of satellite imagery, and the comparison with RGB product, we found that NHI underestimated some volcanic thermal anomalies in presence of clouds/volcanic plumes (e.g., see Section 5.6) partially affecting the target areas and leading to negative values of used indices. On the other hand, background pixels showing radiance values slightly higher at 2.2 µm than 1.6. µm wavelength generated artefacts associated to positive values of the index in Equation (1) (e.g., see Section 5.1). To minimize this issue, an initial test on SWIR radiance might be used before computing the NHI indices as shown in Figure 14, displaying the thermal anomaly map generated, under Google Earth Engine [89] using the preliminary NHI tool, from Landsat OLI data of 24 February, 2019. It can be noted as artefacts in Figure 12c, i.e., false hotspot pixels within the region marked in green, could be efficiently removed optimizing the algorithm (e.g., filtering pixels with values of L 2.2 < 3.0 W m −2 sr −1 µm −1 ). Hence, the confidence level of NHI detections may be further increased; in this direction, we could use also cloud-masking procedure to minimize false positives related to possible cloud-edge effects (see Section 5.5). Since NHI is also sensitive to fires, a spatial distance filter might be integrated within the process to avoid the misinterpretation of detected thermal anomalies at the investigated volcanic areas. Bush fires generally develop along the volcano flanks due to the fallout of incandescent blocks on vegetated areas (e.g., see Section 2.6). Nevertheless, they may also occur because of anthropogenic causes, as for the Vesuvius (Italy) volcanic area on 12 July, 2017. Figure 15 shows the NHI map from Sentinel 2 MSI data of that day, showing the identification of several thermal anomalies on both flanks and slopes of this high-risk volcano, which were generated by a number of well-documented anthropogenic fires burning several hectares of vegetation (e.g., [77]). The use of a spatial filter analyzing distance of detected hotspot pixels from summit craters could then enable a better discrimination of volcanic thermal anomalies from bush fires, especially in remote areas where validation sources usually lack (e.g., [90]). On the other hand, Figure 15 confirms that NHI might potentially be used for mapping thermal anomalies from different sources (e.g., fires; gas flares), although further and more accurate investigations have to be performed to assess this potential, which is out of the scope of this study. for the Vesuvius (Italy) volcanic area on 12 July, 2017. Figure 15 shows the NHI map from Sentinel 2 MSI data of that day, showing the identification of several thermal anomalies on both flanks and slopes of this high-risk volcano, which were generated by a number of well-documented anthropogenic fires burning several hectares of vegetation (e.g., [77]). The use of a spatial filter analyzing distance of detected hotspot pixels from summit craters could then enable a better discrimination of volcanic thermal anomalies from bush fires, especially in remote areas where validation sources usually lack (e.g., [90]). On the other hand, Figure 15 confirms that NHI might potentially be used for mapping thermal anomalies from different sources (e.g., fires; gas flares), although further and more accurate investigations have to be performed to assess this potential, which is out of the scope of this study. The use of spectral/spatial filters described above should make NHI even more suited for mapping volcanic thermal anomalies at a global scale. In this context, the low processing times and the direct exportability to satellite data from other sensors as the Visible Infrared Imaging Radiometer Suite (VIIRS) represent other important advantages offered by the algorithm proposed and tested in this work.  The use of spectral/spatial filters described above should make NHI even more suited for mapping volcanic thermal anomalies at a global scale. In this context, the low processing times and the direct exportability to satellite data from other sensors as the Visible Infrared Imaging Radiometer Suite (VIIRS) represent other important advantages offered by the algorithm proposed and tested in this work.

Conclusions
In this paper, we have presented a single image multi-channel algorithm for mapping volcanic thermal anomalies at a global scale using Sentinel 2 MSI and Landsat 8 OLI data.
Results revealed the NHI potential in mapping both weak (e.g., fumarole fields) and more intense (e.g., lava bodies) volcanic thermal anomalies in spite of some limitations. Moreover, by combining two normalized indices, NHI provided qualitative information about intensity of detected thermal anomalies, in both remote (e.g., Erta Ale) and well-monitored (e.g., Mt. Etna; Stromboli) volcanic areas.
The use of additional spectral/spatial tests should further increase the NHI performance, enabling a more efficient identification of hot volcanic targets in daylight conditions. In this direction, once a sufficiently populated dataset of Sentinel 2 MSI/Landsat 8 OLI scenes will be available, the implementation of NHI within a well-established multi-temporal scheme of satellite data analysis (e.g., [91]) will be evaluated to increase the accuracy of detection, avoiding the use fixed threshold tests which are intrinsically affected by known limitations.
A NHI-based tool, which benefits from the computation capabilities of Google Earth Engine, a cloud-based platform for planetary-scale environmental data analysis whose catalog is continuously updated with a latency of about 24 h from data acquisition time (e.g., [92], is under development. The tool will enable a quick visualization of thermal anomalies flagged by NHI at volcanoes and periods of interest (e.g., see Figure 14), supporting end-users (e.g., volcano observatories) in retrieving information about location, shape, and spatial extent of volcanic thermal anomalies.
Finally, the low processing times make NHI (which could be potentially used for investigating hot targets from different sources) suited for operational contexts, by exploiting the increased revisit time of the combined Landsat 8 and Sentinel-2A/2B satellites. This way, NHI could complement observations from current satellite-based monitoring systems (e.g., those analyzing high temporal resolution satellite data), contributing to the surveillance of active volcanoes from space.

Funding:
The study did not receive any funding.