A Google Earth Engine Tool to Investigate, Map and Monitor Volcanic Thermal Anomalies at Global Scale by Means of Mid-High Spatial Resolution Satellite Data

: Several satellite-based systems have been developed over the years to study and monitor thermal volcanic activity. Most of them use high temporal resolution satellite data, provided by sensors like the Moderate Resolution Imaging Spectroradiometer (MODIS) that if on the one hand guarantee a continuous monitoring of active volcanic areas on the other hand are less suited to map thermal anomalies, and to provide accurate information about their features. The Multispectral Instrument (MSI) and the Operational Land Imager (OLI), respectively, onboard the Sentinel-2 and Landsat-8 satellites, providing Short-Wave Infrared (SWIR) data at 20 m (MSI) and 30 m (OLI) spatial resolution, may make an important contribution in this area. In this work, we present the ﬁrst Google Earth Engine (GEE) App to investigate, map and monitor volcanic thermal anomalies at global scale, integrating Landsat-8 OLI and Sentinel-2 MSI observations. This open tool, which implements the Normalized Hot spot Indices (NHI) algorithm, enables the analysis of more than 1400 active volcanoes, with very low processing times, thanks to the high GEE computational resources. Performance and limitations of the tool, such as its next upgrades, aiming at increasing the user-friendly experience and extending the temporal range of data analyses, are analyzed and discussed.


Introduction
Every year on average about 50-70 active volcanoes erupt, posing a serious threat at both local and global scale (e.g., [1,2]). Pyroclastic flows, lava flows, ashfall, gas emissions, tsunami and lahars induced by volcanic eruptions, may kill people and destroy infrastructure and vegetated areas (e.g., [2][3][4][5][6]). During strong explosive eruptions, huge quantities of gas and ash particles may be injected into the atmosphere, with a global effect on climate change (e.g., [4,7]). Volcanic ash may also cause significant air traffic disruptions, because of potential damage to aircraft engines, and to the exposed surfaces of the aircrafts [8,9]. Systems capable of guaranteeing a continuous and efficient monitoring of these phenomena are fundamental to mitigate volcanic risk.
In this direction, several methods have been developed over the years to monitor thermal volcanic activity in near-real time (e.g., [10][11][12][13][14][15][16]), exploiting advantages of satellite remote sensing, in terms of global coverage, and of multi-spectral data provided at different spatial/temporal resolutions.
In a previous paper [35], we tested the NHI algorithm over a number of active volcanoes located in different geographic areas, retrieving information about thermal anomalies in good agreement with independent field and satellite-based observations. These results encouraged us to implement the algorithm within the GEE environment, as discussed in detail in Section 2.3.

The Google Earth Engine (GEE) Platform
Google Earth Engine is a powerful cloud-based geo-spatial processing platform. It is addressed to a wide audience (e.g., scientists, policy makers, non-profit users), to facilitate the access to high-performance computing resources, for the analysis of huge georeferenced datasets. The GEE platform helps researchers in disseminating and sharing their results. Once an algorithm is developed, users can deliver systematic data products or deploy interactive applications, backed by the GEE's resources [40]. The GEE platform allows for the processing of long and densely populated time-series of satellite imagery, thanks to the access to different Earth-Observing (EO) data catalogues (e.g., U.S. Geological Survey-Landsat), and to the policies of EO imagery access (e.g., the Full, Open and Free data access policy of European Commission-European Space Agency Sentinel). In detail, GEE consists of large data catalogue archives (open and free accessible) of geospatial datasets, and in a set of library functions, which are optimized for the parallel processing of geospatial data. Users can handle both the pre-processed ready-to-use public data and their own data, through the internet-accessible application programming interfaces (API), and/or the web-based interactive development environment (IDE). The GEE data catalogue includes: (i) historical and current EO satellite images, which are updated on a daily basis; (ii) scientific datasets, such as meteorological and climatological data (e.g., climate forecasts); (iii) geophysical (e.g., land cover data, digital elevation models) and socio-economic data [41]. A wide range of functions compose the GEE's library. They include simple mathematical functions, common image processing operations, and algorithms developed for specific datasets (e.g., calibration algorithm for Landsat data). Although the GEE platform allows users to perform analyses at planetary scale, it imposes some limitations to avoid the monopolization of the system. They include the maximum duration of requests, the number of simultaneous requests per user, and the number of simultaneous executions for particular operations. A wide description of the GEE cloud computing environment, and a list of studies exploiting the GEE capabilities in different application fields, can be found in [41].

NHI Tool
The NHI tool, by exploiting the high computational resources and the extended public data archive of the GEE platform, enables the worldwide identification and mapping of volcanic thermal anomalies. The tool is accessible online, without any authentication, at the URL https://nicogenzano. users.earthengine.app/view/nhi-tool and through the open data portal of the Institute of Methodology for Environmental Analysis (IMAA) webpage [42].
Since the beta version (v.1.0) [35], the NHI tool has been developed as a GEE-App, including a group of features (i.e., GEE Feature Collection), corresponding to about 1400 Holocene's volcanoes listed by the Global Volcanism Program (GVP) [43] (see Figure 1). The NHI tool v.1.1 was made available online in January 2020, and included an initial test on SWIR2 radiance to avoid the generation of false positives due to differences in signal to noise ratio (SNR) and minimum radiance of used spectral bands (e.g., [44][45][46]). By testing the tool over different volcanic areas, we observed that false detections mostly affected S2-MSI scenes, in correspondence to cloud-edges mainly over sea areas. By the manual inspection of satellite imagery, we found that these features were related to multispectral misregistration of S2-MSI imagery, affecting high-altitude moving objects (e.g., planes, clouds) [47]. To minimize the impact of S2-MSI data issues on thermal anomaly identification, we implemented a spectral filter within the tool, which combined a normalized index analyzing radiance measured in bands 12 and 8A of the MSI sensor to a spectral test on visible radiance (see v.1.2 in Table  1). This spectral filter reduced the number of false positives at the expense of an increase in missed detections (see also Section 3.5). Hence, we better tuned the filter with the next tool release (v.1.3); the latter included an additional spectral test to better map very intense thermal anomalies yielding the saturation of both SWIR channels (see test for "extreme pixels" in Table 1). The NHI tool v.1.4 used here is similar to the previous one, apart from a function directly integrating L8-OLI and S2-MSI observations (see Section 3.2).
As for previous tool versions, users may access the client-side of the GEE App through the graphical user interface (GUI). In detail, users may run the process through the "Start NHI Analysis" button after selecting the volcano of interest, from a drop-down menu, and setting the temporal range and the distance buffer (DB). The system generates two charts of hot spot pixel number, one for each data collection (see plot 1 in Figure 2), except when, due to the imposed computation limits of the GEE platform, a high distance buffer (i.e., > 50 km for Landsat-8; > 30 km for Sentinel-2 data) is used [41]. Users may then visualize the additional plots of total SWIR radiance, which is calculated by summing radiances of detected hot spot pixels at the relative SWIR wavelength (plot 2 in Figure 2), and of total hot spot area (see Section 3.2). These products may be download in three different formats, i.e., comma-separated values (CSV), Scalable Vector Graphics (SVG) and Portable Network Graphic (PNG), for the offline analysis of results. Once the system generates these charts, users may overlap thermal anomalies to a static very-high resolution image or a topographic map of the volcano. Finally, the "Layers menu" area allows users to visualize five different layers including: i) the investigated area (white circle); ii) three groups of hot spot pixels, classified as mid-low intensity pixels (in yellow), high-intensity pixels (in red) and extreme pixels (in violet) (see Table 1); iii) the RGB product of the selected satellite scene (see Figure 2).
In the next section, we assess the capacity of the NHI tool (v.1.4) in investigating changes of thermal volcanic activity, in monitoring the space-time evolution of active lava flows, and in precisely localizing active vents, as a complement to information from other systems. To minimize the impact of S2-MSI data issues on thermal anomaly identification, we implemented a spectral filter within the tool, which combined a normalized index analyzing radiance measured in bands 12 and 8A of the MSI sensor to a spectral test on visible radiance (see v.1.2 in Table 1). This spectral filter reduced the number of false positives at the expense of an increase in missed detections (see also Section 3.5). Hence, we better tuned the filter with the next tool release (v.1.3); the latter included an additional spectral test to better map very intense thermal anomalies yielding the saturation of both SWIR channels (see test for "extreme pixels" in Table 1). The NHI tool v.1.4 used here is similar to the previous one, apart from a function directly integrating L8-OLI and S2-MSI observations (see Section 3.2). Table 1. Description of NHI tool versions. The test for extreme pixels takes into account the nominal saturation radiances of Short-Wave Infrared (SWIR) channel at 1.6 µm wavelength (e.g., 71.3 W/m 2 sr µm; for the Operational Land Imager (OLI)) (see [48,49]); (* tests applied only to pixels flagged by the index in Equation (2)). ND12-8A is a normalized index, using bands 12 and 8A of the Multispectral Instrument (MSI) sensor. As for previous tool versions, users may access the client-side of the GEE App through the graphical user interface (GUI). In detail, users may run the process through the "Start NHI Analysis" button after selecting the volcano of interest, from a drop-down menu, and setting the temporal range and the distance buffer (DB). The system generates two charts of hot spot pixel number, one for each data collection (see plot 1 in Figure 2), except when, due to the imposed computation limits of the GEE platform, a high distance buffer (i.e., >50 km for Landsat-8; >30 km for Sentinel-2 data) is used [41].

of 22
Users may then visualize the additional plots of total SWIR radiance, which is calculated by summing radiances of detected hot spot pixels at the relative SWIR wavelength (plot 2 in Figure 2), and of total hot spot area (see Section 3.2). These products may be download in three different formats, i.e., comma-separated values (CSV), Scalable Vector Graphics (SVG) and Portable Network Graphic (PNG), for the offline analysis of results. Once the system generates these charts, users may overlap thermal anomalies to a static very-high resolution image or a topographic map of the volcano. Finally, the "Layers menu" area allows users to visualize five different layers including: (i) the investigated area (white circle); (ii) three groups of hot spot pixels, classified as mid-low intensity pixels (in yellow), high-intensity pixels (in red) and extreme pixels (in violet) (see Table 1); (iii) the RGB product of the selected satellite scene (see Figure 2).   In the next section, we assess the capacity of the NHI tool (v.1.4) in investigating changes of thermal volcanic activity, in monitoring the space-time evolution of active lava flows, and in precisely localizing active vents, as a complement to information from other systems.

Identification of Unrest/Quiescience Periods of Active Volcanoes
To assess the capacity of the NHI tool in providing information about unrest and quiescence periods of active volcanoes, we show in Figure 3 the results of the Tolbachik (Kamchatka, Russia) volcano investigation, performed by analyzing L8-OLI and S2-MSI data of April 2013-May 2020.
target area (see thermal anomaly map in Figure 3). Thermal anomaly in Figure 3 is in agreement with results of a previous independent study, where authors tested and compared three different hot spot detection methods, including the more recent HOTMAP algorithm [31]. Moreover, information from the tool is consistent with the Advanced Land Imager ALI observation of 6 June 2013, showing the presence of a lava pond in the active vent, a fluid lava in the center of lava flow, and some multiple flow-front breakouts [52].  Table 1).
It is worth mentioning that the tool identified also the inner part of the main flow, while it slightly overestimated its borders, due to data artifacts (see red pixels surrounding NHI detections), possibly associated to the same effects discussed in [39]. In addition, it missed a few hot spot pixels (4-5%; see areas marked by the white circles) with SWIR2 radiance values slightly higher than the background. After June 2013, the hot spot pixel number drastically decreased also because of clouds, which strongly affected the target area, marking the reduction of thermal volcanic activity. Finally, results from L8-OLI data show that some thermal anomalies occurred also after the end of eruption (e.g., [53]). These features, which did not represent false positives (see also Table 2), may be ascribed to possible lava cooling effects recorded before the new quiescence phase of the volcano.  Table 1).
The most recent Tolbachick eruption began on 27 November 2012, after about 37 years of quiescence, and ended in mid-September 2013 [50]. Effusive activity was particularly intense especially during the first two days of the eruption, when the lava effusion rate was around 465 m 3 /s [51]. In late December 2012-June 2013, almost constant values of this parameter (~19 m 3 /s on average) were recorded. Afterwards, the lava effusion rate gradually decreased, until the end of the eruption [51,52].
Results from L8-OLI data (see chart on top-left side of the Figure 3) show the identification of a thermal anomaly at the Tolbachick volcano, associated to the above-mentioned lava effusion, from April 2013. In detail, the peak of hot spot pixel number (1482 pixels in total) refers to the satellite scene of 5 June 2013, when a major lava flow, moving in the SE direction, and some smaller flows affected the target area (see thermal anomaly map in Figure 3). Thermal anomaly in Figure 3 is in agreement with results of a previous independent study, where authors tested and compared three different hot spot detection methods, including the more recent HOTMAP algorithm [31]. Moreover, information from the tool is consistent with the Advanced Land Imager ALI observation of 6 June 2013, showing the presence of a lava pond in the active vent, a fluid lava in the center of lava flow, and some multiple flow-front breakouts [52].
It is worth mentioning that the tool identified also the inner part of the main flow, while it slightly overestimated its borders, due to data artifacts (see red pixels surrounding NHI detections), possibly associated to the same effects discussed in [39]. In addition, it missed a few hot spot pixels (~4-5%; see areas marked by the white circles) with SWIR2 radiance values slightly higher than the background. After June 2013, the hot spot pixel number drastically decreased also because of clouds, which strongly affected the target area, marking the reduction of thermal volcanic activity. Finally, results from L8-OLI data show that some thermal anomalies occurred also after the end of eruption (e.g., [53]). These features, which did not represent false positives (see also Table 2), may be ascribed to possible lava cooling effects recorded before the new quiescence phase of the volcano.
Regarding the outputs from S2-MSI data, they represented false positives associated to the striping of SWIR bands (e.g., [54]) and to residual multi-spectral misregistration issues, which were not Remote Sens. 2020, 12, 3232 7 of 22 completely removed by the spectral filter detailed in Table 1 (see Figure 4). The impact of these data issues on NHI tool performance is quantified in Section 3.5.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 22 Regarding the outputs from S2-MSI data, they represented false positives associated to the striping of SWIR bands (e.g., [54]) and to residual multi-spectral misregistration issues, which were not completely removed by the spectral filter detailed in Table 1 (see Figure 4). The impact of these data issues on NHI tool performance is quantified in Section 3.5.

Assessing Advantages of Data Integration in Monitoring Active Volcanoes
One of the main advantages in using S2-MSI and L8-OLI data is their integration that, if properly exploited, may increase the frequency of observation of active volcanic areas. Nevertheless, due to a number of factors (e.g., differences in spectral response functions), an effective integration of those data is generally not easy to perform, and may require specific analyses (e.g., corrections for atmospheric effects) to provide reliable results (e.g., [55][56][57]).
The NHI tool performs a rough integration of L8-OLI and S2-MSI observations through the computation of the total hot spot area, calculated from the hot spot pixel number and the pixel size. Figure 5 displays the chart of this parameter (see green curve) retrieved at Karangetang (Indonesia), where satellite observations represent an important source of information (e.g., [58]). Moreover, the figure shows the thermal anomaly identified on the S2-MSI scene of 30 December 2018, corresponding to the peak of total hot spot area, overlapped to the topographic map of the volcano. To emphasize advantages of data integration, Figure 6 displays curves from L8-OLI and S2-MSI data in two different colors. It is worth noting that the NHI tool well outlined the different phases of thermal volcanic activity, by combining different satellite observations. This is evident looking for instance at the period January 2017−November 2018, when a small thermal anomaly (900-2700 m 2 in size) affected, with a certain continuity, the target area. On 30 December 2018, a detected thermal anomaly covered a larger area (see map in Figure 5), revealing the presence of a lava flow, SW directed and extending up to about 1100 m elevation, for which we did not find independent ground-based observations.

Assessing Advantages of Data Integration in Monitoring Active Volcanoes
One of the main advantages in using S2-MSI and L8-OLI data is their integration that, if properly exploited, may increase the frequency of observation of active volcanic areas. Nevertheless, due to a number of factors (e.g., differences in spectral response functions), an effective integration of those data is generally not easy to perform, and may require specific analyses (e.g., corrections for atmospheric effects) to provide reliable results (e.g., [55][56][57]).
The NHI tool performs a rough integration of L8-OLI and S2-MSI observations through the computation of the total hot spot area, calculated from the hot spot pixel number and the pixel size. Figure 5 displays the chart of this parameter (see green curve) retrieved at Karangetang (Indonesia), where satellite observations represent an important source of information (e.g., [58]). Moreover, the figure shows the thermal anomaly identified on the S2-MSI scene of 30 December 2018, corresponding to the peak of total hot spot area, overlapped to the topographic map of the volcano. To emphasize advantages of data integration, Figure 6 displays curves from L8-OLI and S2-MSI data in two different colors. It is worth noting that the NHI tool well outlined the different phases of thermal volcanic activity, by combining different satellite observations. This is evident looking for instance at the period January 2017-November 2018, when a small thermal anomaly (~900-2700 m 2 in size) affected, with a certain continuity, the target area. On 30 December 2018, a detected thermal anomaly covered a larger area (see map in Figure 5), revealing the presence of a lava flow, SW directed and extending up to about 1100 m elevation, for which we did not find independent ground-based observations. Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 22 The tool flagged a similar thermal anomaly on some satellite scenes of January−February and August−November 2019; the smaller hot spot identified from January 2020 marked the reduction of thermal volcanic activity (see Figure 6). It should be stressed that the value of total hot spot area was generally comparable when retrieved from L8-OLI and S2-MSI observations close in time (e.g., as for 30 November 2018, when this parameter ranged in between 28,800 and 29,600 m 2 ). These results indicate that users may analyze the chart of total hotspot area to infer information about changes of thermal volcanic activity, despite some limitations discussed in Section 4.

Lava Flows Monitoring and Mapping
The NHI tool allows users to map and monitor areas inundated by lava integrating L8-OLI and S2-MSI observations. Figure 7 shows a sequence of thermal anomaly maps referring to the recent Stromboli (Italy) eruption of July 2019. This eruptive event, occurring after the strong paroxysm of 3  The tool flagged a similar thermal anomaly on some satellite scenes of January−February and August−November 2019; the smaller hot spot identified from January 2020 marked the reduction of thermal volcanic activity (see Figure 6). It should be stressed that the value of total hot spot area was generally comparable when retrieved from L8-OLI and S2-MSI observations close in time (e.g., as for 30 November 2018, when this parameter ranged in between 28,800 and 29,600 m 2 ). These results indicate that users may analyze the chart of total hotspot area to infer information about changes of thermal volcanic activity, despite some limitations discussed in Section 4.

Lava Flows Monitoring and Mapping
The NHI tool allows users to map and monitor areas inundated by lava integrating L8-OLI and S2-MSI observations. Figure 7 shows a sequence of thermal anomaly maps referring to the recent Stromboli (Italy) eruption of July 2019. This eruptive event, occurring after the strong paroxysm of 3 The tool flagged a similar thermal anomaly on some satellite scenes of January-February and August-November 2019; the smaller hot spot identified from January 2020 marked the reduction of thermal volcanic activity (see Figure 6). It should be stressed that the value of total hot spot area was generally comparable when retrieved from L8-OLI and S2-MSI observations close in time (e.g., as for 30 November 2018, when this parameter ranged in between 28,800 and 29,600 m 2 ). These results indicate that users may analyze the chart of total hotspot area to infer information about changes of thermal volcanic activity, despite some limitations discussed in Section 4.

Lava Flows Monitoring and Mapping
The NHI tool allows users to map and monitor areas inundated by lava integrating L8-OLI and S2-MSI observations. Figure 7 shows a sequence of thermal anomaly maps referring to the recent Stromboli (Italy) eruption of July 2019. This eruptive event, occurring after the strong paroxysm of Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 22 July (e.g., [35]), was previously analyzed through a satellite multi-sensor approach, which enabled the identification of eruption onset, the mapping of lava flows and the estimate of lava effusion rate [59]. In more detail, the figure shows that in the initial phase of the effusive eruption the lava flow, moving towards Sciara del Fuoco, was between 700 and 800 m elevations. On 12 July, it extended up to 500 m elevation, approaching five days later the Stromboli's coasts. On 23 July, clouds masked the distal part of the lava flow, which appeared spatially less extended than previous days. On 27 July, the lava flow field was mapped in between 600 and 800 m elevations; on the same day, some small fires probably occurred, as indicated by thermal anomalies flagged on the opposite flank of the volcano (see sparsely hot spot pixels in yellow). Three days later although clouds partially masked the lava flow, it was possible to map its front, which extended up to about 300 m elevation (see map of 30 July).
This example shows that the NHI tool could be used to assess outputs of numerical models simulating lava flow paths, and to support civil protection activities devoted to mitigating the effects of effusive eruptions on neighboring populated areas. The effects of clouds on thermal anomaly identification, the frequency of the combined L8-OLI and S2-MSI observations (see Section 1.1), and the lag of data availability within the GEE environment, are important factors to take into consideration when the NHI tool is used for these purposes.

Localization of Active Vents and Thermal Anomaly Assessment
The accurate localization of active vents, and the assessment of thermal anomalies flagged by other systems, are other relevant advantages offered by the GEE App presented and tested in this work. Figure 8 emphasizes this potential, displaying the thermal anomaly map from Landsat-8 OLI data of 28 June 2016, covering the Mt. Etna (Sicily; Italy) area. The figure shows the identification of a thermal anomaly of a low-mid intensity level (see yellow pixels) at the Voragine (VOR) crater, where a small vent of ~20-30 m in diameter, opened up on 7 August 2016 emitting high-temperature gases (i.e., up to over 600 C) [33]. In detail, the chart of hot spot pixel number shows that since the end of June 2016, i.e., after the powerful paroxysms of May, a less intense thermal anomaly affected the Mt. Etna crater area [33]. Thermal volcanic activity showed some intensity fluctuations (see chart of total SWIR radiance on bottom-left side of the figure), increasing a few days before the opening of new degassing vent at VOR. In more detail, the figure shows that in the initial phase of the effusive eruption the lava flow, moving towards Sciara del Fuoco, was between 700 and 800 m elevations. On 12 July, it extended up to 500 m elevation, approaching five days later the Stromboli's coasts. On 23 July, clouds masked the distal part of the lava flow, which appeared spatially less extended than previous days. On 27 July, the lava flow field was mapped in between 600 and 800 m elevations; on the same day, some small fires probably occurred, as indicated by thermal anomalies flagged on the opposite flank of the volcano (see sparsely hot spot pixels in yellow). Three days later although clouds partially masked the lava flow, it was possible to map its front, which extended up to about 300 m elevation (see map of 30 July).
This example shows that the NHI tool could be used to assess outputs of numerical models simulating lava flow paths, and to support civil protection activities devoted to mitigating the effects of effusive eruptions on neighboring populated areas. The effects of clouds on thermal anomaly identification, the frequency of the combined L8-OLI and S2-MSI observations (see Section 1), and the lag of data availability within the GEE environment, are important factors to take into consideration when the NHI tool is used for these purposes.

Localization of Active Vents and Thermal Anomaly Assessment
The accurate localization of active vents, and the assessment of thermal anomalies flagged by other systems, are other relevant advantages offered by the GEE App presented and tested in this work. Figure 8 emphasizes this potential, displaying the thermal anomaly map from Landsat-8 OLI data of 28 June 2016, covering the Mt. Etna (Sicily; Italy) area. The figure shows the identification of a thermal anomaly of a low-mid intensity level (see yellow pixels) at the Voragine (VOR) crater, where a small vent of~20-30 m in diameter, opened up on 7 August 2016 emitting high-temperature gases (i.e., up to over 600 • C) [33]. In detail, the chart of hot spot pixel number shows that since the end of June 2016, i.e., after the powerful paroxysms of May, a less intense thermal anomaly affected the Mt. Etna crater area [33]. Thermal volcanic activity showed some intensity fluctuations (see chart of total SWIR radiance on bottom-left side of the figure), increasing a few days before the opening of new degassing vent at VOR. These results are in agreement with the outcomes of a previous study, where RST VOLC detections, from infrared AVHRR and MODIS records, were assessed through the offline analysis of S2-MSI and L8-OLI data [33]. They show that the NHI tool may be used for an accurate localization of active vents also in periods of a less intense thermal volcanic activity (e.g., in the presence of thermal anomalies of a few MW), starting from alerts of systems monitoring active volcanoes in near-real time, by means of high-temporal resolution satellite data (see Section 1).
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 22 These results are in agreement with the outcomes of a previous study, where RSTVOLC detections, from infrared AVHRR and MODIS records, were assessed through the offline analysis of S2-MSI and L8-OLI data [33]. They show that the NHI tool may be used for an accurate localization of active vents also in periods of a less intense thermal volcanic activity (e.g., in the presence of thermal anomalies of a few MW), starting from alerts of systems monitoring active volcanoes in near-real time, by means of hightemporal resolution satellite data (see Section 1).

Quantification of the NHI Tool Performance
To quantify the performance of the NHI tool, we analyzed several active volcanoes with recent/no recent eruptions, over the same period April 2013 (Nov 2015 for S2-MSI data)−May 2020 (i.e., the same as Tolbachick in Figure 4). In particular, we estimated the false positive rate (FPR) according to equation (3), i.e., by the ratio of images with evidences of false detections f , revealed by the manual inspection of satellite imagery and information from volcanological reports, and the total number of available images for the volcanic area of interest: Table 2 shows that the FPR was particularly low when retrieved from L8-OLI data, regardless of used tool version. In more detail, the FPR ranged in between 0.0% at Acatenango (Guatemala) and 1.7% at Campi Flegrei (Italy), where false positives were ascribable to possible fires occurring in the surroundings of Pozzuoli town. The FPR was generally higher from S2-MSI data, reaching about 19.2% at Campi Flegrei when the NHI tool v.1.1 was used. In this high-risk volcanic area, which is characterized by a hydrothermal activity and whose last eruption occurred in 1538 AD (e.g., [60,61]), the tool v.1.2 was the most effective in reducing false positives occurring over cloud-edges and associated to S2-MSI data issues. The NHI tool v.1.3/1.4 was slightly less efficient in minimizing these issues (e.g., see Tolbachick). In general, the capacity of these tool versions in reducing false positives related to S2-MSI data issues was more evident at a higher distance buffer (DB), when both land and sea pixels were analyzed.

Quantification of the NHI Tool Performance
To quantify the performance of the NHI tool, we analyzed several active volcanoes with recent/no recent eruptions, over the same period April 2013 (Nov 2015 for S2-MSI data)-May 2020 (i.e., the same as Tolbachick in Figure 4). In particular, we estimated the false positive rate (FPR) according to Equation (3), i.e., by the ratio of images with evidences of false detections Ni f , revealed by the manual inspection of satellite imagery and information from volcanological reports, and the total number of available images Ni for the volcanic area of interest: Table 2 shows that the FPR was particularly low when retrieved from L8-OLI data, regardless of used tool version. In more detail, the FPR ranged in between 0.0% at Acatenango (Guatemala) and 1.7% at Campi Flegrei (Italy), where false positives were ascribable to possible fires occurring in the surroundings of Pozzuoli town. The FPR was generally higher from S2-MSI data, reaching about 19.2% at Campi Flegrei when the NHI tool v.1.1 was used. In this high-risk volcanic area, which is characterized by a hydrothermal activity and whose last eruption occurred in 1538 AD (e.g., [60,61]), the tool v.1.2 was the most effective in reducing false positives occurring over cloud-edges and associated to S2-MSI data issues. The NHI tool v.1.3/1.4 was slightly less efficient in minimizing these issues (e.g., see Tolbachick). In general, the capacity of these tool versions in reducing false positives related to S2-MSI data issues was more evident at a higher distance buffer (DB), when both land and sea pixels were analyzed. Regarding the false negative rate (FNR), we estimated this parameter according to Equation (4), i.e., by the ratio of satellite imagery with missed detections Ni m and the number of satellite scenes with evidences of thermal anomalies Ni T , revealed by the manual inspection of satellite imagery and the comparison with the RGB product: The first tool version, including only the initial test on SWIR radiance (see Table 1), showed the lowest FNR. The tool v.1.2 was more affected by missed detections. As an example, in areas such as Shishaldin (Alaska), the FNR increased about 25% in comparison with the previous tool version, while at Oldoynio Lengai (OL), which is the only volcano on Earth emitting natrocarbonatite lava [62], the increase in FNR was around 3.5% when compared to v.1.3/1.4 (FNR ≈ 52.5%). Therefore, the current NHI tool release (i.e., v.1.4) seems to offer the best trade-off between reliability and sensitivity, as indicated by the analysis of other active volcanoes (e.g., Mt Etna), which are not reported in Table 2. A different filter is however required to increase the confidence level of detection from S2-MSI data, without losing in sensitivity (e.g., FNR ≈ 44.4% at OL from L8-OLI data).

Discussion
By coupling an efficient single image algorithm to the high computational resources of the GEE platform, the NHI tool is capable of providing accurate information about volcanic thermal anomalies, with low processing times. Only a few seconds are, in fact, required to generate the tool outputs; the computational time may increase up to a few minutes when long-time series of satellite observations, at a high distance buffer, are analyzed (e.g., charts of Figure 3 were generated in less than one minute).
In this work, we investigated a number of active volcanoes characterized by a different thermal activity using a different DB as input, based on their relative size. At Campi Flegrei, which is a caldera of about 12 km in diameter [63], we used for instance a DB value of 6 km (see Table 2). The NHI tool helps users in setting the DB through a preview of the selected volcanic area. Nevertheless, in some volcanic regions this parameter may be more difficult to be set properly. As an example, by setting a DB ≥ 3.0 km at Acatenango (Guatemala), whose last eruption occurred in 1972 [64], we retrieved the same results for Fuego volcano, showing frequent eruptions and belonging to the same volcanic complex (Acatenango-Vulcán de Fuego massif). On the other hand, at Kilauea (Hawaii) a high DB value (50 km) was required to investigate the Lower East Rift Zone (LERZ) [65] (see Figure 9). In this case, the developed platform generated only outputs from L8-OLI data because of GEE restrictions (see Section 2).
This study has also shown that the NHI tool performs well in detecting and mapping lava flows, as observed at Nishinoshima (Japan) which erupted in November 2013 after 40 years of dormancy emitting lava [66], even thanks to the spectral test for "extreme pixels" described in Table 1. This study has also shown that the NHI tool performs well in detecting and mapping lava flows, as observed at Nishinoshima (Japan) which erupted in November 2013 after 40 years of dormancy emitting lava [66], even thanks to the spectral test for "extreme pixels" described in Table 1.
In more detail, we observed the underestimation of lava flows/lakes under certain conditions (e.g., pixel saturation effects, presence of clouds), while their slight overestimation was recorded mainly because of data artefacts as shown in Figure 3, and confirmed by the analysis of S2-MSI data. Figure 10 shows the lava flow from Sierra Negra (Ecuador) eruption of 16 July 2018 [67], confirming the low impact of data artefacts on thermal anomaly identification and the improvements achieved using the spectral test for extreme pixels (see pixels in violet in the inner part of lava flow).  In more detail, we observed the underestimation of lava flows/lakes under certain conditions (e.g., pixel saturation effects, presence of clouds), while their slight overestimation was recorded mainly because of data artefacts as shown in Figure 3, and confirmed by the analysis of S2-MSI data. Figure 10 shows the lava flow from Sierra Negra (Ecuador) eruption of 16 July 2018 [67], confirming the low impact of data artefacts on thermal anomaly identification and the improvements achieved using the spectral test for extreme pixels (see pixels in violet in the inner part of lava flow).
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 22 This study has also shown that the NHI tool performs well in detecting and mapping lava flows, as observed at Nishinoshima (Japan) which erupted in November 2013 after 40 years of dormancy emitting lava [66], even thanks to the spectral test for "extreme pixels" described in Table 1.
In more detail, we observed the underestimation of lava flows/lakes under certain conditions (e.g., pixel saturation effects, presence of clouds), while their slight overestimation was recorded mainly because of data artefacts as shown in Figure 3, and confirmed by the analysis of S2-MSI data. Figure 10 shows the lava flow from Sierra Negra (Ecuador) eruption of 16 July 2018 [67], confirming the low impact of data artefacts on thermal anomaly identification and the improvements achieved using the spectral test for extreme pixels (see pixels in violet in the inner part of lava flow).  Volcanic thermal features of a lower intensity level (e.g., hot fumaroles) may be more difficult to identify with the tool (sensitivity limits of the NHI algorithm will be assessed in a next work). Table 3 details L8-OLI with evidences of thermal anomalies at OL where mid-low intensity hot spots were prevalent, while Figure 11 shows an example of intra-crater thermal anomaly undetected by the tool. Table 3. L8-OLI data of April 2013-May 2020, with evidences of thermal anomalies at Oldoynio Lengai (Tanzania; Africa) revealed by the manual inspection of satellite imagery. Volcanic thermal features of a lower intensity level (e.g., hot fumaroles) may be more difficult to identify with the tool (sensitivity limits of the NHI algorithm will be assessed in a next work). Table  3 details L8-OLI with evidences of thermal anomalies at OL where mid-low intensity hot spots were prevalent, while Figure 11 shows an example of intra-crater thermal anomaly undetected by the tool. Table 3. L8-OLI data of April 2013−May 2020, with evidences of thermal anomalies at Oldoynio Lengai (Tanzania; Africa) revealed by the manual inspection of satellite imagery.  Despite the FNR higher than other areas, such as Suwanoseijima (Japan) where the FNR from L8-OLI data was estimated around 13%, the occurrence of some lava effusions inside the OL crater (e.g., [68]) indicates that events of natrocarbonatite lava, having temperatures in between 500 °C and 600 °C Despite the FNR higher than other areas, such as Suwanoseijima (Japan) where the FNR from L8-OLI data was estimated around 13%, the occurrence of some lava effusions inside the OL crater (e.g., [68]) indicates that events of natrocarbonatite lava, having temperatures in between 500 • C and 600 • C (e.g., [69]), may be detected/assessed by the tool. On the other hand, different satellite-based systems (e.g., those using TIR data) need to be used to identify low-temperature thermal anomalies, as for Campi Flegrei, where fumaroles reach temperatures up to 160 • C-165 • C (e.g., [61,70]).

Landsat-8 OLI data (YYYYMMDD) Detected (Yes/No
Among the outputs of the NHI tool, the chart of total SWIR radiance (which is generated without performing any correction for atmospheric effects), such as that of total hot spot area (see Section 3.2), require a careful analysis and interpretation. Figure 12 shows the chart of total SWIR radiance retrieved, from L8-OLI data of January 2016-April 2020, over the Nyiragongo (Africa) volcano, which hosts the largest active lava lake in the world [71]. It can be noted as despite the lower number of pixels flagged Remote Sens. 2020, 12, 3232 14 of 22 by the index in Equation (2) than Equation (1) (Figure 12a), the total SWIR1 radiance was systematically higher than SWIR2 (Figure 12b). This signal behavior, which is different from that recorded in other areas (e.g., see Figure 8), may be ascribed to the most intense hot spot pixels detected only by the index in Equation (2). These pixels yielded the OLI channel 7 saturation and showed radiance values on OLI channel 6 above the nominal saturation and close to the maximum observable radiance (96 W/m 2 sr µm) (e.g., [37]). Concerning the chart of total hot spot area, factors such as cloud coverage may lead to a misinterpretation of results even analyzing near contemporaneous satellite observations. In addition, the total hot spot area represents an overestimation of the area actually covered by hot targets (the tool considers each hot spot pixel as covered only by the hot component). Users should be aware of these limitations when the above-mentioned products are used.
Among the outputs of the NHI tool, the chart of total SWIR radiance (which is generated without performing any correction for atmospheric effects), such as that of total hot spot area (see Section 3.2), require a careful analysis and interpretation. Figure 12 shows the chart of total SWIR radiance retrieved, from L8-OLI data of January 2016−April 2020, over the Nyiragongo (Africa) volcano, which hosts the largest active lava lake in the world [71]. It can be noted as despite the lower number of pixels flagged by the index in equation (2) than equation (1) (Figure 12a), the total SWIR1 radiance was systematically higher than SWIR2 (Figure 12b). This signal behavior, which is different from that recorded in other areas (e.g., see Figure 8), may be ascribed to the most intense hot spot pixels detected only by the index in equation (2). These pixels yielded the OLI channel 7 saturation and showed radiance values on OLI channel 6 above the nominal saturation and close to the maximum observable radiance (96 W/m 2 sr m) (e.g., [37]). Concerning the chart of total hot spot area, factors such as cloud coverage may lead to a misinterpretation of results even analyzing near contemporaneous satellite observations. In addition, the total hot spot area represents an overestimation of the area actually covered by hot targets (the tool considers each hot spot pixel as covered only by the hot component). Users should be aware of these limitations when the above-mentioned products are used.

Future Perspectives
Thanks to the features of the GEE platforms, further improvements of the NHI tool are possible. These improvements will allows us to take into account the suggestions from users, provided through a feedback module directly linked to the tool; e.g., we could ingest more regions of interest within the tool to overcome limitations in S2-MSI data processing occurring at a high distance buffer.

Future Perspectives
Thanks to the features of the GEE platforms, further improvements of the NHI tool are possible. These improvements will allows us to take into account the suggestions from users, provided through a feedback module directly linked to the tool; e.g., we could ingest more regions of interest within the tool to overcome limitations in S2-MSI data processing occurring at a high distance buffer.
Some tool improvements under evaluation include: (i) the implementation of new functions to increase the user-friendly experience; (ii) the use of a procedure to correct SWIR radiances for atmospheric effects (e.g., [72]); (iii) the implementation of methods for thermal anomaly characterization; (iv) the ingestion of additional datasets to extend the temporal range of satellite data analyses.
Among the next tool upgrades, the use of Thematic Mapper (TM)/ Enhanced Thematic Mapper plus (ETM+) and ASTER data, could make more than 30 years of satellite observations available to users, enabling the analysis of worldwide volcanoes before the first availability of L8-OLI data. Figure 13 emphasizes this potential, showing the preliminary results of the Mount St. Helens (USA) and Eyjafjallajökull (Iceland) investigations, achieved using Landsat-5 TM and Landsat-7 ETM+ data, respectively. In detail, Figure 13a shows that an intra-crater thermal anomaly affected the Mount St. Helens on 14 October 2014. This feature was generated by a documented lava dome growth inside the crater [73,74], affecting other satellite scenes of 2004-2006 (see chart of hot spot pixel number referring to period 1985-2010). Figure 13b displays the lava flow emitted by the Eyjafjallajökull volcano on 3 April 2010 and generated by an intense effusive eruption which started the previous month after about two centuries of quiescence, previously investigated from space by means of infrared MODIS data [75]. It should be noted that, despite the striping in Landsat-7 ETM+ imagery (e.g., [76]), the analyzed time series (from 1999 to 2011) did not reveal any false detection (see quiescence periods correctly marked by the tool before and after the end of eruption). Thus, the NHI tool seems to perform well in using also TM/ETM+ data, although the impact of other issues (e.g., oversaturation effects [77]) needs to be better assessed. Additionally, even though the malfunctioning of ASTER SWIR detectors does not enable the analysis of satellite imagery acquired after 2008 (see Section 1), the ingestion of ASTER data within the tool could contribute to better analyzing past eruptions, thanks to the better dynamic range in the SWIR region (e.g., [78]).
Although the NHI tool was designed for volcanological applications, a similar GEE App could be developed to investigate and map other hot targets. Figure 14 shows an example of this potential, in reference to the gas flares, with temperatures in between 1600 and 2200 K [79]. In particular, the figure displays the chart of hot spot pixel number retrieved, at the Val d'Agri Oil Center (COVA) in Basilicata region (Southern Italy), from L8-OLI data of 2015-2020. The figure also shows the thermal anomaly identified on a satellite scene of 14 December 2019, corresponding to the peak of hot spot pixels, revealing the presence of two hot areas. In previous studies, gas-flaring activity at COVA was investigated and quantified using nighttime MODIS and Visible Infrared Imaging Radiometer Suite (VIIRS) data [80,81]. Information from L8-OLI/S2-MSI imagery could then integrate observations from other systems, allowing for estimates of the fire radiative power (FRP), and of volume of emitted gas, in daylight conditions (e.g., [82,83] Some tool improvements under evaluation include: i) the implementation of new functions to increase the user-friendly experience; ii) the use of a procedure to correct SWIR radiances for atmospheric effects (e.g., [72]); iii) the implementation of methods for thermal anomaly characterization; iv) the ingestion of additional datasets to extend the temporal range of satellite data analyses. Among the next tool upgrades, the use of Thematic Mapper (TM)/ Enhanced Thematic Mapper plus (ETM+) and ASTER data, could make more than 30 years of satellite observations available to users, enabling the analysis of worldwide volcanoes before the first availability of L8-OLI data. Figure 13 emphasizes this potential, showing the preliminary results of the Mount St. Helens (USA) and Eyjafjallajökull (Iceland) investigations, achieved using Landsat-5 TM and Landsat-7 ETM+ data, respectively. In detail, Figure 13a shows that an intra-crater thermal anomaly affected the Mount St. Helens on 14 October 2014. This feature was generated by a documented lava dome growth inside the crater [73,74], affecting other satellite scenes of 2004−2006 (see chart of hot spot pixel number referring to period 1985−2010). Figure 13b displays the lava flow emitted by the Eyjafjallajökull volcano on 3 April 2010 and generated by an intense effusive eruption which started the previous month after about two centuries of quiescence, previously investigated from space by means of infrared MODIS data [75]. It should be noted that, despite the striping in Landsat-7 ETM+ imagery (e.g., [76]), the analyzed time series (from 1999 to 2011) did not reveal any false detection (see quiescence periods correctly marked by the tool before and after the end of eruption). Thus, the NHI tool seems to perform well in using also TM/ETM+ data, although the impact of other issues (e.g., oversaturation effects [77]) needs to be better assessed. Additionally, even though the malfunctioning of ASTER SWIR detectors does not enable the analysis of satellite imagery acquired after 2008 (see Section 1.1), the ingestion of ASTER data within the tool could contribute to better analyzing past eruptions, thanks to the better dynamic range in the SWIR region (e.g., [78]). Although the NHI tool was designed for volcanological applications, a similar GEE App could be developed to investigate and map other hot targets. Figure 14 shows an example of this potential, in reference to the gas flares, with temperatures in between 1600 and 2200 K [79]. In particular, the figure displays the chart of hot spot pixel number retrieved, at the Val d'Agri Oil Center (COVA) in Basilicata region (Southern Italy), from L8-OLI data of 2015-2020. The figure also shows the thermal anomaly identified on a satellite scene of 14 December 2019, corresponding to the peak of hot spot pixels, revealing the presence of two hot areas. In previous studies, gas-flaring activity at COVA was investigated and quantified using nighttime MODIS and Visible Infrared Imaging Radiometer Suite (VIIRS) data [80,81]. Information from L8-OLI/S2-MSI imagery could then integrate observations from other systems, allowing for estimates of the fire radiative power (FRP), and of volume of emitted gas, in daylight conditions (e.g., [82,83]). Although the NHI tool was designed for volcanological applications, a similar GEE App could be developed to investigate and map other hot targets. Figure 14 shows an example of this potential, in reference to the gas flares, with temperatures in between 1600 and 2200 K [79]. In particular, the figure displays the chart of hot spot pixel number retrieved, at the Val d'Agri Oil Center (COVA) in Basilicata region (Southern Italy), from L8-OLI data of 2015-2020. The figure also shows the thermal anomaly identified on a satellite scene of 14 December 2019, corresponding to the peak of hot spot pixels, revealing the presence of two hot areas. In previous studies, gas-flaring activity at COVA was investigated and quantified using nighttime MODIS and Visible Infrared Imaging Radiometer Suite (VIIRS) data [80,81]. Information from L8-OLI/S2-MSI imagery could then integrate observations from other systems, allowing for estimates of the fire radiative power (FRP), and of volume of emitted gas, in daylight conditions (e.g., [82,83]). Since the NHI algorithm is sensitive to wildfires (e.g., [35]), we could develop a similar GEE App also to map forest fires at local/regional scale. Figure 15 shows forest fires (in red) identified on S2-MSI scene of 15 January 2020, affecting the New South Wales region (Australia). These features were among the numerous, large and devastating fire events occurring during the 2019-2020 Australian bushfire season (e.g., [84,85]). Although these features were slightly underestimated, as indicated by the manual inspection of satellite imagery, and even if other investigations are required to assess NHI performance under different fire regimes, the GEE App of Figure 15 could contribute in assessing risk of wildland-urban interface fires, complementing information from operational fire monitoring systems (e.g., [86]).

Eyjafjallajökull (Iceland)
Since the NHI algorithm is sensitive to wildfires (e.g., [35]), we could develop a similar GEE App also to map forest fires at local/regional scale. Figure 15 shows forest fires (in red) identified on S2-MSI scene of 15 January 2020, affecting the New South Wales region (Australia). These features were among the numerous, large and devastating fire events occurring during the 2019−2020 Australian bushfire season (e.g., [84,85]). Although these features were slightly underestimated, as indicated by the manual inspection of satellite imagery, and even if other investigations are required to assess NHI performance under different fire regimes, the GEE App of Figure 15 could contribute in assessing risk of wildlandurban interface fires, complementing information from operational fire monitoring systems (e.g., [86]).

Conclusions
In this paper, we have presented a powerful tool to analyze, map and monitor volcanic thermal anomalies at global scale, by means of medium-high spatial resolution satellite data.
The NHI tool, even if does not provide automated alerts about ongoing eruptions unlike other systems (e.g., [34]), is currently the only platform which allows for the investigation of more than 1400 active volcanic areas, through the integration of L8-OLI and S2-MSI observations. The tool, whose performance is strictly related to the used algorithm and the GEE computational resources, shows a low false positive rate especially from L8-OLI data, in spite of some limitations (e.g., slight dependence on data issues/artifacts; sensitivity reduction in presence of clouds/plumes). The next tool upgrades will include a new visualization mode to increase the user-friendly experience, the possible ingestion of some independent satellite products for volcano monitoring available in GEE (e.g., SO2 maps from Sentinel-5P), and the implementation of the prior Landsat-5/7 data collections to extend the temporal range of satellite data analyses. Moreover, the possible use of VIIRS data, available at 500 m spatial resolution in the Imagery resolution (I)-bands within GEE, is also under evaluation. A suite of multi-platform satellite observations will be then developed to analyze both past and ongoing eruptions, with an increased temporal coverage. A similar GEE App could be adapted to other hot targets (e.g., forest fires), having a strong impact on both environment and population.

Conclusions
In this paper, we have presented a powerful tool to analyze, map and monitor volcanic thermal anomalies at global scale, by means of medium-high spatial resolution satellite data.
The NHI tool, even if does not provide automated alerts about ongoing eruptions unlike other systems (e.g., [34]), is currently the only platform which allows for the investigation of more than 1400 active volcanic areas, through the integration of L8-OLI and S2-MSI observations. The tool, whose performance is strictly related to the used algorithm and the GEE computational resources, shows a low false positive rate especially from L8-OLI data, in spite of some limitations (e.g., slight dependence on data issues/artifacts; sensitivity reduction in presence of clouds/plumes). The next tool upgrades will include a new visualization mode to increase the user-friendly experience, the possible ingestion of some independent satellite products for volcano monitoring available in GEE (e.g., SO2 maps from Sentinel-5P), and the implementation of the prior Landsat-5/7 data collections to extend the temporal range of satellite data analyses. Moreover, the possible use of VIIRS data, available at 500 m spatial resolution in the Imagery resolution (I)-bands within GEE, is also under evaluation. A suite of multi-platform satellite observations will be then developed to analyze both past and ongoing eruptions, with an increased temporal coverage. A similar GEE App could be adapted to other hot targets (e.g., forest fires), having a strong impact on both environment and population.
In conclusion, this paper shows that the NHI tool may highly support scientists, volcano observatories, particularly those in developing countries where the expertise in data processing is often lacking [87], and even simple users (e.g., students), in investigating active volcanoes, complementing information from other systems.