The Transition from MODIS to VIIRS for Global Volcano Thermal Monitoring

The Moderate Resolution Imaging Spectroradiometer (MODIS) is one of the most-used sensors for monitoring volcanoes and has been providing time series of Volcanic Radiative Power (VRP) on a global scale for two decades now. In this work, we analyzed the data provided by the Visible Infrared Imaging Radiometer Suite (VIIRS) by using the Middle Infrared Observation of Volcanic Activity (MIROVA) algorithm, originally developed to analyze MODIS data. The resulting VRP is compared with both the MIROVAMODIS data as well as with the Fire Radiative Power (FRP), distributed by the Fire Information for Resource Management System (FIRMS). The analysis on 9 active volcanoes reveals that VIIRS data analyzed with the MIROVA algorithm allows detecting ~60% more alerts than MODIS, due to a greater number of overpasses (+30%) and improved quality of VIIRS radiance data. Furthermore, the comparison with the nighttime FIRMS database indicates greater effectiveness of the MIROVA algorithm in detecting low-intensity (<10 MW) thermal anomalies (up to 90% more alerts than FIRMS). These results confirm the great potential of VIIRS to complement, replace and improve MODIS capabilities for global volcano thermal monitoring, because of the future end of Terra and Aqua Earth-observing satellite mission of National Aeronautics and Space Administration’s (NASA).


Introduction
Volcanic activity causes the variation of numerous geophysical and geochemical parameters that characterize the state of a volcano. Among these, the satellite-derived thermal flux is increasingly used to detect signs of unrest [1][2][3][4] or to follow the evolution of an eruption once the activity has started [5][6][7].
The first applications of this discipline date back to the 60s and 70s through the analysis of the High-Resolution Infrared Radiometer (HRIR), a sensor mounted on NASA's meteorological satellites Nimbus-1 and Nimbus-2 [8]. Later in the mid-1980s, several studies demonstrated how the data provided by the Thematic Mapper (TM) of the Landsat missions [9,10] and by the Advanced Very High-Resolution Radiometer (AVHRR) of the Tiros missions [11] can be successfully applied for studying thermal emissions; these sensors have been used even recently, for example, to quantify thermal budgets of eruptions [12] and to build thermal satellite monitoring tools [13][14][15]. The development of algorithms capable of detecting and monitoring the different types of activities allowed to expand this field of study [16], also thanks to the availability of thermal data acquired from sensors such as the Along Track Scanning Radiometer (ATSR) and the Advanced Baseline Imager (ABI) mounted on the Geostationary Operational Environmental Satellite (GOES) [9,11,17,18].
representing one of the most valid candidates to continue the analyses of volcanic thermal flux on a global scale [46]. VIIRS applications for detecting and studying global wildfire activity and its impacts are well-known and validated, also by merging data from different satellites see [47]. However, in the last years, several authors use VIIRS-derived Fire Radiative Power (FRP) data from freely-downloadable Fire Information for Resource Management System (FIRMS; [48]) to study volcanic eruptions and activity [49][50][51]. However, there is still no work that has validated the FIRMS data for volcanological applications. Recently, few authors worked on the development of algorithms and systems dedicated to the thermal monitoring application of volcanic activity by using VIIRS data [52][53][54], but none of them is either operative on a global scale yet.
In this work, we test the exportability of the MIROVA algorithm, currently ingesting NRT MODIS images [55], to VIIRS images. The nighttime VIIRS data acquired during 2012-2020 at two persistent active volcanoes, Láscar (Chile) and Erta Ale (Ethiopia), were processed using the MIROVA algorithm to build time series of VRP (here called VRP VIIRS ). This dataset was then compared with that obtained by the same algorithm applied to the MODIS data (VRP MODIS ) as well as with that provided by the FIRMS algorithm [56] (this last resulting in time-series of Fire Radiative Power; FRP VIIRS ). The three datasets are then analysed in terms of the number of alerts, the statistical distribution of VRP/FRP, and the cumulative thermal energies.
Finally, we expand our analysis by testing the MIROVA VIIRS algorithm on a temporally limited dataset (1 year: 2021) which includes both data (day and night) acquired by the two VIIRS sensors currently in orbit (see Section 3.1). This test is performed over on 9 volcanoes characterized by very different types of activity (e.g., lava flows, lava lakes, lava domes, open vent activities) and located in very different environmental contexts ( Figure 1). These preliminary results, provide a solid basis to permit the continuity of MIROVA and others hot-spot detection systems through the analysis of VIIRS data toward a global thermal volcano monitoring for the next decade(s). to MODIS thus representing one of the most valid candidates to continue the analyses of volcanic thermal flux on a global scale [46]. VIIRS applications for detecting and studying global wildfire activity and its impacts are well-known and validated, also by merging data from different satellites see [47]. However, in the last years, several authors use VIIRS-derived Fire Radiative Power (FRP) data from freely-downloadable Fire Information for Resource Management System (FIRMS; [48]) to study volcanic eruptions and activity [49][50][51]. However, there is still no work that has validated the FIRMS data for volcanological applications. Recently, few authors worked on the development of algorithms and systems dedicated to the thermal monitoring application of volcanic activity by using VIIRS data [52][53][54], but none of them is either operative on a global scale yet.
In this work, we test the exportability of the MIROVA algorithm, currently ingesting NRT MODIS images [55], to VIIRS images. The nighttime VIIRS data acquired during 2012-2020 at two persistent active volcanoes, Láscar (Chile) and Erta Ale (Ethiopia), were processed using the MIROVA algorithm to build time series of VRP (here called VRPVIIRS). This dataset was then compared with that obtained by the same algorithm applied to the MODIS data (VRPMODIS) as well as with that provided by the FIRMS algorithm [56] (this last resulting in time-series of Fire Radiative Power; FRPVIIRS). The three datasets are then analysed in terms of the number of alerts, the statistical distribution of VRP/FRP, and the cumulative thermal energies.
Finally, we expand our analysis by testing the MIROVAVIIRS algorithm on a temporally limited dataset (1 year: 2021) which includes both data (day and night) acquired by the two VIIRS sensors currently in orbit (see Section 3.1). This test is performed over on 9 volcanoes characterized by very different types of activity (e.g., lava flows, lava lakes, lava domes, open vent activities) and located in very different environmental contexts (Figure 1). These preliminary results, provide a solid basis to permit the continuity of MIROVA and others hot-spot detection systems through the analysis of VIIRS data toward a global thermal volcano monitoring for the next decade(s).

Case Studies
The thermal activity investigated in this work is principally those retrieved on Láscar and Erta Ale volcanoes ( Figure 1). We choose these two targets for an intercomparison of the three datasets (MIROVA VIIRS , MIROVA MODIS , FIRMS VIIRS ) due to the following motivation: 1.
both volcanoes were in persistent activity for the entire duration of the analyzed period (2012-2020). This makes thermal datasets large enough to allow a robust correlation (>3000 images analyzed at each volcano); 2.
both volcanoes are located in arid areas with low cloud fraction, which favors the high alert detection frequency, and reduces the noise in radiative power time series; 3.
both volcanoes are characterized by evident fluctuations of heat flux during the considered period. This feature allows to evaluate how the VRP/FRP responds to sudden variations of the heat flux over several orders of magnitude; 4.
the two volcanoes are characterized by different kinds of activity (lava lake and basaltic lava flow for Erta Ale, and persistent fumarolic degassing and explosive activity for Láscar); thus we may evaluate the algorithm's performance on volcanic thermal sources having different intensity and temperature distributions.
In the next paragraphs, we give an overview of the geological background and recent activity of the investigated volcanoes.

Láscar
Láscar (23.37 • S-67.73 • W) is a 5592 m high andesitic-dacitic stratovolcano located in the northern Chilean Andes ( Figure 2) and it's the most active volcano of the Andean Central Volcanic Zone (ACVZ) [9]. The ACVZ is formed by the subduction of the Nazca plate under the South American plate in a back-arc context [57].

Erta Ale
Erta Ale (13.6° N-40.67° E) is the most active volcano in Ethiopia and one of the most active in Africa ( Figure 3). This basaltic shield volcano lies on the about 120 km-NNW elongated Erta 'Ale spreading segment within the Afar depression region, which represents the triple junctions between Nubian, Somalian and Arabian plane plate [74,75]. The 600 m-high edifice is cut by a 0.9 × 1.6 km elliptical summit caldera which hosts two pit craters [76,77]. For almost a century one or two lava lakes have been observed within these pit-crater(s), which made Erta Ale the volcano hosting the longest-lived lava lake(s) on Earth [76,78,79]. However, in January 2017, massive overflows from the lava lake inside the south pit-crater were followed by the first observed fissure eruption, who took place outside the caldera, along the southeast rift zone (Figure 3b,c). The effusive activity persisted until February-March 2020 and produced a large lava field covering ~28 km 2 The edifice of Láscar consists of two overlapping cones and currently hosts three summit craters that partially overlap. The westernmost of these craters is the one currently active, also known as "Active Crater" [58][59][60]. In the last century, this volcano has been characterized by cyclic lava dome extrusions and collapses, accompanied by Vulcanian to Plinian eruptions [58,61]. The last major eruption (VEI 4) occurred between January and August 1993 and was characterized by a series of phreatic to Plinian explosions that produced columns up to 25 km high and more than 0.1 km 3 of ash [62]. Since then, Láscar  [60,63,64]. The activity of the last decade shows a thermal pattern defining three yearly-long cycles having a similar trend. These cycles have been described by [65] and partially observed in dedicated periodic reports [66].
Because of its geographical conditions, the activity of Láscar volcano has been investigated mainly using satellite remote sensing data [9,10,63,[67][68][69][70]. According to these authors, the source of thermal emission is located mainly at the bottom of the Active Crater (800 m wide, 400 m deep) where extensive fumarolic areas (having average temperature estimated between 300 • C [59,71,72] and 600 • C [73]) have been observed, eventually associated to the presence of permeable andesitic lava dome(s).

Erta Ale
Erta Ale (13.6 • N-40.67 • E) is the most active volcano in Ethiopia and one of the most active in Africa ( Figure 3). This basaltic shield volcano lies on the about 120 km-NNW elongated Erta 'Ale spreading segment within the Afar depression region, which represents the triple junctions between Nubian, Somalian and Arabian plane plate [74,75]. The 600 m-high edifice is cut by a 0.9 × 1.6 km elliptical summit caldera which hosts two pit craters [76,77]. For almost a century one or two lava lakes have been observed within these pit-crater(s), which made Erta Ale the volcano hosting the longest-lived lava lake(s) on Earth [76,78,79]. However, in January 2017, massive overflows from the lava lake inside the south pit-crater were followed by the first observed fissure eruption, who took place outside the caldera, along the southeast rift zone (Figure 3b,c). The effusive activity persisted until February-March 2020 and produced a large lava field covering~28 km 2 (data updated June 2019) [77,[80][81][82]. This effusive phase has been investigated using various sensors (i.e., MSI Sentinel-2 and OLI Landsat-8), to evaluate new algorithms and new multi-sensor approaches for monitoring volcanic activity (cf with [83][84][85]). On the whole, all the authors agree in associating the persistent thermal emissions recorded until 2017 to the lava lake(s) within the pit crater(s) and/or to the occurrence of episodic overflows [80,[86][87][88]. On the other hand, the highest thermal anomalies detected between January 2017 and March 2020 are related to the emplacement of the lava flows that occurred during the exceptional flank eruption. Since March 2020, several satellite monitoring systems identified discontinuous thermal activity sourced a new from the two pit craters [82].

Sensor and Products
The Visible Infrared Imaging Radiometer Suite (VIIRS) is a whiskbroom scanning radiometer mounted on two platforms: the Suomi National Polar-Orbiting Partnership's (SNPP) and the Joint Polar Satellite System's (JPSS) JPSS-1 (NOAA-20 or N20) orbiting since October 2011 and November 2017, respectively [89,90]. The Suomi-NPP mission is a "bridge" mission born to continue NASA's Earth Observing System (EOS) missions

Sensor and Products
The Visible Infrared Imaging Radiometer Suite (VIIRS) is a whiskbroom scanning radiometer mounted on two platforms: the Suomi National Polar-Orbiting Partnership's (SNPP) and the Joint Polar Satellite System's (JPSS) JPSS-1 (NOAA-20 or N20) orbiting since October 2011 and November 2017, respectively [89,90]. The Suomi-NPP mission is a "bridge" mission born to continue NASA's Earth Observing System (EOS) missions (which included Terra and Aqua satellites). VIIRS installation's key objective is to guarantee a continuous acquisition of the parameters acquired by MODIS. Both satellites (Suomi-NPP and JPSS-1) allow full daily coverage of Earth thanks to their polar orbit, at a nominal altitude of 824 km [45,91]. VIIRS field-of-view (FOV) is 112.56 • (13 • higher than MODIS) in the cross-track direction and its swath width is 3060 km. Different from MODIS, the onboard acquisition scheme of VIIRS is characterized by an aggregation function that reduces pixel footprint growth across the scan by aggregating two or three raw samples. This function increases the pixel resolution at the edge of the images (Table 1) and extends the possibility of high-quality acquisitions over a higher satellite zenith angle [45]. VIIRS acquires in 22 spectral bands covering the EM spectrum between 0.412 µm and 12.01 µm, including 16 moderate-resolution bands (M-bands), with spatial resolution of 750 m, 5 imaging resolution bands (I-bands), with a spatial resolution of 375 m and one panchromatic Day-Night Band (DNB), with a 750 m spatial resolution throughout the scan. 11 of the 16 M-bands are Reflective Solar Bands (RSBs), 5 are Thermal Emission Bands (TEBs), while the I-bands include 3 RSBs and 2 TEBs [45].
In this work we used M13 and M15 M-bands, covering respectively the MIR (3.973-4.128 µm) and TIR (10.263-11.263 µm) portion of the spectrum, with a 750 m spatial resolution, which represent the corresponding MODIS bands 21/22 (MIR) and 31 (TIR), respectively. Table 1 shows the main characteristics of the VIIRS and MODIS sensors [92,93].
For this work, we used Suomi-NPP VIIRS Level 1B products which include the radiances data product (VNPMOD02-Moderate Resolution 6-Min L1B Swath 750 m) and the geolocation data product (VNPMOD03-Moderate Resolution Terrain Corrected Geolocation 6-Min L1 Swath 750 m Light V001). These two products are distributed by NASA's Level-1 and Atmosphere Archive & Distribution System-Distributed Active Archive Center (LAADS-DAAC) [94], in netCDF4/HDF5 format files with sizes of~200 MB and~55 MB, respectively. We elaborated only the nigh time data acquired by Suomi-NPP from 19th January 2012 to 31st December 2020, which consists of a bulk dataset of about 0.5 Tb per volcano. To make the comparison with MODIS as homogeneous as possible we only used MODIS-Aqua nighttime acquisitions elaborated by MIROVA over the same time window.

MIROVA Algorithm and Calculation of Volcanic Radiative Power (VRP)
The Middle InfraRed Observation of Volcanic Activity (MIROVA) is a near-real-time monitoring system based on the processing of MODIS InfraRed data [28]. The system operates in near-real-time (NRT) for a list of~220 active volcanoes and provides VRP time series used for monitoring ongoing volcanic activity [7]. The algorithm elaborates satellite images using spectral and spatial filters in search for high-temperature thermal anomalies, sourced by volcanic features. The NRT processing chain is made of 4 successive steps: (i) download; (ii) resampling; (iii) hot-spot detection and (iv) calculation of the VRP. As described below, steps 1 to 4 have been slightly modified from [28] and adapted to VIIRS due to its spatial, geometric, and spectral differences with MODIS (Table 1).
We downloaded VIIRS L1B granules from the LAADS-DAAC system [94] and we extracted MIR (band M13) and TIR (band M15) radiances and the geolocation data related to our target volcanoes. Resampling is performed in a UTM 51 × 51 km grid, centered on the volcano summit (consistent with MODIS MIROVA images) by keeping the nominal resolution of 750 m. This results in matrices of 67 × 67 pixels rather than 51 × 51 pixels obtained from MODIS (Figures 2 and 3).
The hot-spot detection algorithm is the same used for MODIS and based on spectral indices and statistical thresholds that automatically identify and localize pixels containing high-temperature thermal anomalies (see [28] for details). Once the hot-spot-contaminated pixels have been detected, we calculate the VRP by adapting the MIR method [31] to the VIIRS features: where A pix is the pixel surface in km 2 (equal to 0.5625 for VIIRS M-bands). The constant 1.97 × 10 7 represents the best-fit wavelength-specific coefficient that takes into account the relationship between the Stefan-Boltzmann law and Plank's Radiation law for a range of temperature between 600 and 1600 K [31,47]. In (1), the parameter ∆L MIR represents the excess of MIR radiance, calculated as: where L MIRhot is the radiance of the alerted pixel/s and L MIRbk is the background radiance, calculated as the arithmetic mean of pixels surrounding the active one/s [28,30,47]. The VRP is calculated in Watts (W) and represents a combined measurement of the area of the volcanic emitter having an effective radiating temperature higher than 600 K [28].

FIRMS Database and Fire Radiative Power (FRP FIRMS )
The two volcanoes were analyzed by using the thermal anomalies provided by the FIRMS system over the same regions and time windows. In particular, we download the Suomi-NPP alerts, detected by the new VIIRS 375 m active fire detection algorithm [95], which combine the two MIR bands (at 750 m and 375 m) available on the VIIRS sensors. These data are archived since January 2012 and downloadable at https://firms.modaps. eosdis.nasa.gov/ (accessed on 16 December 2021), through temporal and spatial queries.
The downloaded data consists of a single file (.csv) which includes the Fire Radiative Power (FRP VIIRS ) values for each hot-spot contaminated pixel alerted in the period and the region of interest. Therefore, to be consistent with the VRP measured by the MIROVA system, for each satellite overpasses we aggregated the FRP of individual pixels to obtain the total FRP FIRMS value for that image. It should be emphasized that the hot-spot detection algorithm underlying the FIRMS VIIRS system [95] is different from MIROVA, this last being expressly developed to work on volcanic areas. However, the method of calculating the FRP VIIRS is essentially the same as that used to calculate the VRP VIIRS , both based on the MIR method described by Equation (1).

Uncertainties and Limits
The estimation of VRP through the MIROVA algorithm, or FRP through the FIRMS algorithm, is subject to several limitations that condition their direct application for volcano monitoring.
First of all, there is a standard error of ±30% associated with the estimation of the radiant power through the MIR method [31]. This error affects the estimates of VRP/FRP, although the coefficients used in Equation (1) Table 1). Additional uncertainty in the estimation of the VRP/FRP may derive from the method used to detect the hot pixels (i.e., the hot-spot detection algorithm) as well as in the estimation of the method to estimate background radiance values (Equation (2)).
Several external factors may also contribute to increasing the uncertainty of the VRP. Indeed, the presence of meteorological and volcanic clouds, the satellite viewing angle, while false alerts, fires, or anthropogenic heat sources, may also introduce noise in the time series and a higher uncertainty in the single-point interpretation. For this reason, during an eruptive crisis the data provided by MIROVA, as well as that eventually obtained from the FIRMS, must always be evaluated and interpreted by an end-user that takes into account the real acquisition conditions of each hot spot [7]. However, for a comparison between the time series, in this work, we have left the datasets "as they are", that is without applying image inspections or filters that discard cloudy scenes or scenes acquired in unfavorable geometric conditions (e.g., high satellite zeniths). Actually, under these conditions, we test the potential efficiency of the algorithm in NRT applications where such supervision is not applied.

Results
In this chapter, we present the results obtained at Láscar and Erta Ale volcanoes by comparing the three datasets (VRP VIIRS , VRP MODIS, and FRP FIRMS , all in Watt) as obtained by the MIROVA and FIRMS algorithms.

Láscar
The Láscar time series derived by both sensors and algorithms ( Figure 4) are characterized by nearly continuous detection of thermal anomalies resulting in a number of alerts (N alert ) equal to 1221 for MIROVA VIIRS , 875 for MIROVA MODIS, and 627 for FIRMS VIIRS ( Table 2). These data suggest that over Láscar volcano the MIROVA VIIRS algorithm detects 40% more alerts than MIROVA MODIS and 95% more than FIRMS VIIRS . By taking into account the number of MODIS and VIIRS overpasses (N pass ) the above data translates into a frequency of alert detection (f% = N alert /N pass ) of 30%, 28%, and 15% for MIROVA VIIRS , MIROVA MODIS, and FIRMS VIIRS , respectively (see Table 2). These percentages suggest that the higher number of MIROVA VIIRS alerts is not simply attributable to the higher number of VIIRS overpasses, but also to a greater sensitivity of VIIRS data to the presence of small thermal anomalies.
The trends defined by all the time series (Figure 4a-c) are very coherent and show at least three main eruptive cycles characterized by a sudden peaking phase (with VRP ≥ 8 MW), followed by a waning phase lasting several years (with VRP decreasing below 1 MW). These three cycles are part of a major phase of Láscar activity (identified as "Phase 4" in [65]) which began in April 2013 in correspondence with our first cycle. Based on MIROVA VIIRS data analysis, we identify "Cycle 1" between 5 April 2013 and 29 October 2015, "Cycle 2" between 30 October 2015 and 22 November 2018, and "Cycle 3" between 23 November 2018 and 31 December 2020 (still ongoing). The highest VRP values were recorded during Cycle 2, after a phreatic eruption occurred on 30 October 2015 [63], and produced VRP values between 13 MW and 16 MW (see Figure 4 for details). Notably, the time gap between the thermal onset of MIROVA VIIRS and MIROVA MODIS (i.e., a sharp increase of VRP) characterizing each cycle was within 24 h, with a difference in the VRP value reached in the subsequent peak phases of less than 2 MW (Figure 4). This outlines that the three datasets provided consistent indications of abrupt changes in volcanic activity in both the timing (±24 h) and magnitude (±6%) of thermal emissions. However, it can be noted that the monthly-long phases preceding the climax of the three recognized cycles are characterized by few FIRMS alerts, while MIROVA VIIRS detections outline the persistence of a very-low volcano-related thermal activity (Figure 4). The black dashed lines represent the date on which a sudden change in the VRP, associated with the beginning of a cycle, was detected by the three datasets. The arrows highlight the peaking VRP recorded during the following days. Note that, even if the high thermal phases reach their maximum detections on a different day, the inferred onset is detected by all systems within a time window of 24 h, with few differences in magnitude, probably due to instrumental and algorithm-related features.
The statistical analysis of the VRP dataset (log-transformed) may be useful to distinguish different thermal regimes and highlight threshold values separating different types of activity [96][97][98]. For long-term monitoring applications, it is therefore important that the datasets are complementary also in this aspect and that they allow the same thresholds to be calculated based on their statistical distributions.
The frequency distributions of MIROVA-derived VRP for Láscar are very similar for both VIIRS and MODIS, with data showing a range of VRP spanning from 0.1 to about 20 MW (Figure 5a). In general, there is a greater number of MIROVAVIIRS alerts in each class in agreement with the greater number of VIIRS passages compared to MODIS (Table 2). In the probability plots ( Figure 5b) the trends defined by these two datasets are also very similar both defining the presence of two different regimes separated by a threshold value of ~5 MW (logVRP ~6.7 in Figure 5a). We ascribe this threshold as separating two types of activity characterizing Láscar: (1) the near-continuous degassing activity (VRP < 5 MW) characterizing inter-eruptive periods, and (2) the main eruptive phases (VRP > 5 MW) associated with a moderate explosive activity eventually accompanied by episodes of dome extrusion or collapse.
The presence of these two regimes, on the other hand, is not detectable from the FIRMSVIIRS dataset, neither from the data distribution nor from the associated probability plot (Figure 5c,d). Indeed, there is a clear lack of FIRMS detection below 5 MW (logVRP ~6.7 in Figure 5c) suggesting that this algorithm is not able to fully detect the weak thermal anomalies inside the crater of Láscar (associated with high-temperature degassing activity). This result emerges clearly in the probability plot (Figure 5d) where there is a clear underestimation of FIRMSVIIRS measurements for logVRP values <6.7, which does not allow to separate the two distinct thermal regimes outlined by MIROVAVIIRS datasets. . The black dashed lines represent the date on which a sudden change in the VRP, associated with the beginning of a cycle, was detected by the three datasets. The arrows highlight the peaking VRP recorded during the following days. Note that, even if the high thermal phases reach their maximum detections on a different day, the inferred onset is detected by all systems within a time window of 24 h, with few differences in magnitude, probably due to instrumental and algorithm-related features.
To exclude the possibility that these weak anomalies were false, we evaluated their accuracy through a visual inspection of all the alerts detected during 2012-2013 and we verified that they were always located inside the Lascar's crater. The presence of a nearlycontinuous low-level thermal activity during inter-eruptive periods was already retrieved and validated by previous analysis with other sensors (e.g., ASTER analysis; [3]).
The statistical analysis of the VRP dataset (log-transformed) may be useful to distinguish different thermal regimes and highlight threshold values separating different types of activity [96][97][98]. For long-term monitoring applications, it is therefore important that the datasets are complementary also in this aspect and that they allow the same thresholds to be calculated based on their statistical distributions.
The frequency distributions of MIROVA-derived VRP for Láscar are very similar for both VIIRS and MODIS, with data showing a range of VRP spanning from 0.1 to about 20 MW (Figure 5a). In general, there is a greater number of MIROVA VIIRS alerts in each class in agreement with the greater number of VIIRS passages compared to MODIS (Table 2). In the probability plots ( Figure 5b) the trends defined by these two datasets are also very similar both defining the presence of two different regimes separated by a threshold value of~5 MW (logVRP~6.7 in Figure 5a). We ascribe this threshold as separating two types of activity characterizing Láscar: (1) the near-continuous degassing activity (VRP < 5 MW) characterizing inter-eruptive periods, and (2) the main eruptive phases (VRP > 5 MW) associated with a moderate explosive activity eventually accompanied by episodes of dome extrusion or collapse. Table 2. Comparison between the datasets MIROVA VIIRS , MIROVA MODIS and FIRMS VIIRS for Láscar and Erta Ale volcanoes. The summary includes the number of overpasses (N pass ), the number of alerts (N alerts ) including the frequency (f%) of alerts (N alert /N pass ), the N alert ratio (ratio between the number of alerts detected by the different algorithms) as well as the mean and maximum VRP/FRP obtained for each dataset. The presence of these two regimes, on the other hand, is not detectable from the FIRMS VIIRS dataset, neither from the data distribution nor from the associated probability plot (Figure 5c,d). Indeed, there is a clear lack of FIRMS detection below 5 MW (logVRP 6.7 in Figure 5c) suggesting that this algorithm is not able to fully detect the weak thermal anomalies inside the crater of Láscar (associated with high-temperature degassing activity). This result emerges clearly in the probability plot ( Figure 5d) where there is a clear underestimation of FIRMS VIIRS measurements for logVRP values <6.7, which does not allow to separate the two distinct thermal regimes outlined by MIROVA VIIRS datasets. To better correlate the three time series, we calculated for each sensor the average values of VRP on a weekly scale (VRPw), as shown in Figure 6a, on a logarithmic scale. The three weekly datasets have been correlated by calculating the best-fit linear coefficient (m) and the Pearson correlation coefficient (r). From these analyses, it results that the MIROVA VIIRS vs. MIROVA MODIS and MIROVA VIIRS vs. FIRMS VIIRS are characterized by coefficient r equal to~0.82 and~0.87, respectively (Figure 6b,c). Data are slightly dispersed around the 1:1 ratio line, with a linear best-fit coefficient m of~0.90 for MIROVA VIIRS vs. FIRMS VIIRS correlation, suggesting that MIROVA-derived weekly mean are slightly lower than FIRMS-derived FRP (Figure 6b). On the contrary, the coefficient m = 1.05 found for MIROVA-derived VRP (VIIRS vs. MODIS) indicates that the weekly mean for VIIRS is slightly higher than those obtained with MODIS (Figure 6c).

Erta Ale
The long-lasting activity of the Erta Ale is well represented by the continuous alerts detected by both sensors (Figure 7), corresponding to 2711, 1868 and 2816 for MIROVA VIIRS , MIROVA MODIS, and FIRMS VIIRS , respectively (see Table 2). For this volcano, the MIROVA VIIRS has therefore detected~45% more alerts than the MIROVA MODIS and~3% fewer than the FIRMS VIIRS . Taking into consideration the number of satellite overpasses, the three systems had alert frequencies (f%) equal to ca. 70%, for MIROVA VIIRS , 63% MIROVA MODIS , and 72% for FIRMS VIIRS . To better correlate the three time series, we calculated for each sensor the average values of VRP on a weekly scale (VRPw), as shown in Figure 6a, on a logarithmic scale. The three weekly datasets have been correlated by calculating the best-fit linear coefficient (m) and the Pearson correlation coefficient (r). From these analyses, it results that the MIROVAVIIRS vs. MIROVAMODIS and MIROVAVIIRS vs. FIRMSVIIRS are characterized by coefficient r equal to ~0.82 and ~0.87, respectively (Figure 6b,c). Data are slightly dispersed around the 1:1 ratio line, with a linear best-fit coefficient m of ~0.90 for MIROVAVIIRS vs. FIRMSVIIRS correlation, suggesting that MIROVA-derived weekly mean are slightly lower than FIRMS-derived FRP (Figure 6b). On the contrary, the coefficient m=1.05 found for MIROVA-derived VRP (VIIRS vs. MODIS) indicates that the weekly mean for VIIRS is slightly higher than those obtained with MODIS (Figure 6c). According to the eruptive activity that occurred at Erta Ale during the analyzed period, the dataset may be subdivided into two main phases (  Figure 7) probably outlining the occurrence of minor overflows from the lava lake(s). Successively, the sudden increase of VRP (occurred on 17 January 2017, marks the beginning of the lateral fissure eruption [77,81] and was recorded by the sudden increase of VRP by all the three systems (dashed line in Figure 7). The propagation of the lava flow in the following days led the VRP to increase rapidly, reaching peak values of 6.7 GW (MIROVA VIIRS ), 6.8 GW (MIROVA MODIS ), and 4.9 GW (FIRMS VIIRS ) by 21-22 January (Table 2 and Figure 7).  Table 2. Comparison between the datasets MIROVAVIIRS, MIROVAMODIS and FIRMSVIIRS for Láscar and Erta Ale volcanoes. The summary includes the number of overpasses (Npass), the number of alerts (Nalerts) including the frequency (f%) of alerts (Nalert/Npass), the Nalert ratio (ratio between the number of alerts detected by the different algorithms) as well as the mean and maximum VRP/FRP obtained for each dataset.  As shown in Figure 8a,c, Erta Ale's data are characterized by very similar distributions with the highest number of detections occurring between 25 and 32 MW, for all systems (logVRP equal to 7.4-7.5 in Figure 8c). As for Láscar, our data suggest that VIIRS detected a higher number of anomalies than MODIS (+7%; Figure 8a) likely due to the higher number of overpasses (Table 2). On the other hand, the distributions of data from MIROVA VIIRS and FIRMS VIIRS are also very similar (Figure 8b). However, for logVRP < 6.7 MIROVA algorithm seems to be more sensitive, as confirmed by the very low number of FIRMSderived detections below the 1 MW (logVRP = 6 in Figure 8c).

Volcano
All the three datasets outline that about 70% of detections are lower than 50 MW (logVRP equal to~7.7; Figure 8b,d), a thermal regime that we ascribe to the typical lava lake activity characterizing Erta Ale volcano until early 2017 (Figure 7). Two other regimes emerge in all the analyzed datasets and are separated by a threshold of 1 GW (logVRP = 9 in Figure 8b,d). These two regimes are likely associated with the bulk effusive activity (50 MW < VRP < 1 GW; 22% of the data; Figure 8) and with the peaking phases of the effusive activity (VRP > 1 GW; less than 4% of the data). As previously cited, these different phases are well recorded by both sensors and algorithms, also showing a time homogeneity observable in the time series. Figure 7) probably outlining the occurrence of minor overflows from the lava lake(s). Successively, the sudden increase of VRP (occurred on 17 January 2017, marks the beginning of the lateral fissure eruption [77,81] and was recorded by the sudden increase of VRP by all the three systems (dashed line in Figure 7). The propagation of the lava flow in the following days led the VRP to increase rapidly, reaching peak values of 6.7 GW (MIROVAVIIRS), 6.8 GW (MIROVAMODIS), and 4.9 GW (FIRMSVIIRS) by 21-22 January (Table  2 and Figure 7).  Accordingly, also for Erta Ale, there is an excellent correspondence between the three datasets, with a tendency of MIROVA VIIRS to produce weekly values slightly lower than the other systems.

Cumulative Radiant Energy (VRE) Calculation via VRP Datasets
A further element of comparison between the three datasets is provided by the cumulative curves of the radiated volcanic energy (VRE in Figure 10). The VRE is a fundamental parameter since it allows to compare eruptions and during the effusive eruptions it is directly correlated to the volume of erupted lava [30]. It is therefore useful to evaluate if the three systems provide comparable VRE values over time.
VRE was calculated by using the trapezoidal integration of VRP/FRP time series, whose cumulative curves are shown in Figure 10. In Table 3 we reported the VRE values for the two analyzed volcanoes, also subdivided based on the eruptive phases identified above. At Láscar we observe that the VRE VIIRS is about 18% higher than VRE MODIS and 21% lower than VRE FIRMS . These differences are somehow coherent with the correlation coefficient (m) characterizing the weekly values (as shown in Figure 6) but are possibly biased also by the number and distribution of VRP values in each dataset. Indeed, in the case of FIRMS, the lack of many low-intensity detections (Figure 5c) produces an integrated VRE curve higher than in the other datasets. On the other hand, no remarkable variations are observed in the ratio of the VRE throughout the eruption which remains roughly constant during all three cycles (Table 3). As shown in Figure 8a,c, Erta Ale's data are characterized by very similar distributions with the highest number of detections occurring between 25 and 32 MW, for all systems (logVRP equal to 7.4-7.5 in Figure 8c). As for Láscar, our data suggest that VIIRS detected a higher number of anomalies than MODIS (+7%; Figure 8a) likely due to the higher number of overpasses (Table 2). On the other hand, the distributions of data from MIROVAVIIRS and FIRMSVIIRS are also very similar (Figure 8b). However, for logVRP < 6.7 MIROVA algorithm seems to be more sensitive, as confirmed by the very low number of FIRMS-derived detections below the 1 MW (logVRP = 6 in Figure 8c).
All the three datasets outline that about 70% of detections are lower than 50 MW (logVRP equal to ~7.7; Figure 8b,d), a thermal regime that we ascribe to the typical lava lake activity characterizing Erta Ale volcano until early 2017 (Figure 7). Two other regimes emerge in all the analyzed datasets and are separated by a threshold of 1 GW (logVRP = 9 in Figure 8b,d). These two regimes are likely associated with the bulk effusive activity (50 MW < VRP < 1 GW; 22% of the data; Figure 8) and with the peaking phases of the effusive activity (VRP > 1 GW; less than 4% of the data). As previously cited, these different phases are well recorded by both sensors and algorithms, also showing a time homogeneity observable in the time series.  Similarly, at Erta Ale we found that the dataset producing the higher VRE is FIRMS VIIRS, followed by MIROVA MODIS and MIROVA VIIRS (Figure 10b). The relationship between the energies is partially influenced by the coefficient m found for Erta Ale (Figure 9), which is indicative of VRP VIIRS values slightly lower than VRP MODIS and FRP VIIRS . However, we found that most of the discrepancy between the cumulative curves is essentially produced during the effusive phase (Table 3 and Figure 10b). In fact, during the effusive eruption, there was a tendency for FIRMS VIIRS to measure energies slightly higher than the  (Figure 9a). This variation suggests that the type of activity (i.e., lava lake vs. lava flow) influenced the correlation between the analyzed datasets, even if the bulk of three VREs remains within a difference of ±20%. (from less than 1 MW to more than 1 GW). The linear regression analysis (Figure 9b,c) confirms the strong correlation between the three datasets (r~0.94 for both MIROVAVIIRS vs. MIROVAMODIS and MIROVAVIIRS vs. FIRMSVIIRS) with best-fit coefficients equal to m=0.85 (for MIROVAVIIRS vs. FIRMSVIIRS) and m=0.91 (for MIROVAVIIRS vs. MIROVAMODIS). Accordingly, also for Erta Ale, there is an excellent correspondence between the three datasets, with a tendency of MIROVAVIIRS to produce weekly values slightly lower than the other systems.

Cumulative Radiant Energy (VRE) Calculation via VRP Datasets
A further element of comparison between the three datasets is provided by the cumulative curves of the radiated volcanic energy (VRE in Figure 10). The VRE is a fundamental parameter since it allows to compare eruptions and during the effusive eruptions it is directly correlated to the volume of erupted lava [30]. It is therefore useful to evaluate if the three systems provide comparable VRE values over time.
VRE was calculated by using the trapezoidal integration of VRP/FRP time series, whose cumulative curves are shown in Figure 10. In Table 3 we reported the VRE values for the two analyzed volcanoes, also subdivided based on the eruptive phases identified above. At Láscar we observe that the VREVIIRS is about 18% higher than VREMODIS and 21% lower than VREFIRMS. These differences are somehow coherent with the correlation coefficient (m) characterizing the weekly values (as shown in Figure 6) but are possibly biased also by the number and distribution of VRP values in each dataset. Indeed, in the case of FIRMS, the lack of many low-intensity detections (Figure 5c) produces an  (Table 3). Similarly, at Erta Ale we found that the dataset producing the higher VR FIRMSVIIRS, followed by MIROVAMODIS and MIROVAVIIRS (Figure 10b). The relation between the energies is partially influenced by the coefficient m found for Erta Ale (Fi 9), which is indicative of VRPVIIRS values slightly lower than VRPMODIS and FRP However, we found that most of the discrepancy between the cumulative curve essentially produced during the effusive phase (Table 3 and Figure 10b). In fact, du the effusive eruption, there was a tendency for FIRMSVIIRS to measure energies slig higher than the MIROVAVIIRS (Figure 9a). This variation suggests that the type of act (i.e., lava lake vs. lava flow) influenced the correlation between the analyzed data even if the bulk of three VREs remains within a difference of ±20%.

Toward Global Volcano Thermal Monitoring Using VIIRS
The results presented above clearly highlight as the VIIRS and MODIS imagery, processed by the MIROVA system, produce consistent results in terms of VRP values, statistical distribution, the timing of major changes in volcanic activity, and cumulative thermal energy (VRE). These analyses were carried out on two volcanoes (Láscar and Erta Ale) specifically chosen to compare the nighttime long datasets also with the data provided by the FIRMS system, which is based on a different algorithm. However, to test the exportability of the MIROVA VIIRS algorithm for a near-real-time and near-global application (as currently made for MODIS), we have to consider the multitude of environmental contexts in which volcanoes are located, as well as the full potential of the VIIRS suite, now constituted by two sensors each one providing nighttime and daytime images.
For this purpose, we analyzed a one-year-long dataset (all daytime and nighttime images between January-December 2021) acquired over 9 persistently active volcanoes ( Figure 1) by the pair of VIIRS sensors currently onboard Suomi-NPP and NOAA-20 satellites, respectively. The resulting VRP time series (Supplementary material) are then compared with those provided by the operational MIROVA MODIS system that includes both night and daytime imagery from Terra and Aqua satellites [7].
The chosen volcanic targets (Table 4 and Figure 11) show thermal emissions ranging from about 1 MW to more than 1 GW and embed a variety of volcanic activity that includes: silicic lava domes (Shiveluch, Sabancaya, Nevados de Chillan), andesitic lava domes (Bagana), lava lakes (Nyiragongo), open vent (Etna, Manam), and mafic lava flow (Kilauea, Pacaya). Some of these volcanoes have experienced important eruptive phases in the considered period (see Supplementary Material), such as Etna [99,100], Nyiragongo volcano [101], and Kilauea [102]. However, in addition to their recent activity, these volcanoes have been chosen because they are representative of a wide range of environmental and meteorological conditions spanning from ice-covered high-latitude volcanoes (e.g., Sheveluch) to tropical islands (e.g., Bagana, Kilauea). We, therefore, focus on evaluating the performance of the MIROVA VIIRS algorithm at each volcano by analyzing: (i) the number of alerts, (ii) the correlations coefficient between weekly VRP values, (iii) the nighttime vs. daytime performance.
The number of alerts: Overall, our analysis reveals that the VIIRS sensors provide about 30% more overpasses than the MODIS (17487 vs. 13416; Table 4). This is due to the larger swath of VIIRS (Table 1) that allows, in some cases, imaging a target volcano even on multiple adjacent orbits. Although this high number of passages is characterized by frequent acquisitions with very high satellite zeniths (many VIIRS passages are characterized by angles greater than 55 • , the limit of MODIS), the total number of VIIRS alerts is on average 62% more than MODIS (6293 vs. 3883; Table 4), indicating an increased efficiency in detecting hot-spots. It should be emphasized that the increase of VIIRS detections is inversely proportional to the average VRP (Table 4), and growth from 30-50% at highly radiating volcanoes (such as Kilauea, Etna, Nyiragongo, and Pacaya), to more than 100% on volcanoes characterized by low thermal flux (e.g., Manam, Bagana, Sabancaya). This makes the VIIRS more effective than MODIS in monitoring background thermal emissions at persistently active volcanoes (Supplementary Material), but also more suitable for detecting thermal precursors at restless volcanoes [7].
We suggest that the higher efficiency of VIIRS, with respect to MODIS, is related to a combination of improved instrument features such as the higher spatial resolution, the lower noise equivalent temperature (NedT), and the already-mentioned aggregation function that allows keeping the spatial resolution below 1.5 km also in extreme viewing conditions (see Table 1). All these features increase the possibility of detecting small sub-pixel hot-spots and ensure a better location of the alerted pixels, too (Figures 2 and 3).
Correlation Coefficient: The correlation found for every single volcano as well as the total dataset (all volcanoes together) is shown in Figure 11. The Pearson coefficient (r) confirms the good correlation between VIIRS-and MODIS-derived VRP, with values ranging from 0.62 (Sabancaya; Figure 11b) to 0.99 (Shiveluch; Figure 11a). Similarly, the best-fitting coefficient (m) is found to vary between 0.78 (Kilauea; Figure 11h) and 1.14 (Shiveluch; Figure 11a), possibly reflecting the influence of variable environmental, meteorological, and volcanological conditions characterizing such a relationship at each distinct volcano. Taken as a whole the bulk correlation ( Figure 11j) give a Pearson coefficient r = 0.95 and a regression coefficient m = 0.87, delineating an overall discrepancy between VRP VIIRS and VRP MODIS of~13%. These results underline the excellent correlation between the two datasets although it is also evident that local factors can play a role in determining a case-by-case variance in the estimate of weekly VRP (typically less than 20%). Table 4. Comparison between the 2021 datasets of VIIRS and MODIS for selected volcanoes. The summary includes the correlation coefficient shown in Figure 11, the number of overpasses (N pass ), the number (N alerts ), the frequency (f% = N alert /N pass ) of alerts (subdivided into nighttime and daytime data), the ratio between VIIRS and MODIS detections, and the mean VRP obtained by averaging the weekly dataset showed in Figure 11 and Supplementary Material. As we have discussed previously, the local correlation coefficients are likely conditioned by the number and distribution of the data, by the type of volcanic activity (stable or fluctuating) as well as by the environmental conditions that may condition the performance of the underlying algorithms.
Nighttime vs daytime data: The comparative analysis confirms that the performance of MIROVA algorithm decreases when applied to the daytime data. In particular, we observed a reduction in the frequency of alerts (f% = N alert / N pass ) from 51% (nighttime) to 20% (daytime) for the VIIRS, compared to a reduction from 43% (nighttime) to 16% (daytime) measured with the MODIS ( Table 4). As described in [26], this strong reduction is due to the higher thresholds of the daytime algorithm specifically settled to avoid false alerts due to the contribution of solar heating and reflection. The reduction is particularly evident in volcanoes characterized by a general low heat emission (e.g., Manam, Bagana, Sabancaya) for which thermal anomalies are detected mostly at night.

Conclusions
In this paper, we presented VRP time series, and cumulative VRE data obtained on two selected target volcanoes, by applying an adapted version of the MIROVA algorithm to VIIRS imagery. We thus analyzed the efficiency of the MIROVA VIIRS algorithm by comparing the obtained time series with those provided by the operative MIROVA MODIS algorithm [55] as well as by the FRP time series provided by the FIRMS system [48]. The comparative analysis at these two specific volcanoes allowed us to conclude that: • by using the MIROVA algorithm, the VIIRS sensor detects~40% more alerts than MODIS (on average). This difference is likely due to the greater number of VIIRS overpasses, compared to MODIS (+29% on average), but also to a better quality of the VIIRS images which make the hot-spot detection more efficient (e.g., better spatial resolution, better pixel aggregation at high satellite zenith angle, less NEdT); • the two MIROVA-derived datasets are highly correlated each other (r = 0.81-0.94 for Láscar and Erta Ale, respectively) and show a best-fit linear coefficient (m = VRP VIIRS /VRP MODIS ) around the 1:1 ratio, ranging from 0.91 (Erta Ale) to 1.05 (Láscar); the two datasets are also comparable in terms of VRP distributions (modal value ±5%), the timing of major events (within 24 h), and cumulative radiant energy (±18%); • the comparison of the VIIRS data processed with FIRMS algorithms instead reveals a better ability of the MIROVA algorithm to detect small thermal anomalies (<10 MW). At volcanoes characterized by low-amplitude thermal anomalies (such as Láscar), this difference translates into an increase of up to 95% of the number of alerts detected by MIROVA VIIRS compared to the FIRMS VIIRS . On the other hand, the two algorithms, appear to be equally efficient on volcanoes characterized by a more intense thermal activity (e.g., Erta Ale); Finally, to test the global application of the MIROVA VIIRS algorithm, we analyzed one year of nighttime and daytime data acquired by the two pairs of VIIRS and MODIS sensors. This analysis has been performed at 9 volcanoes located in different environments and characterized by very different volcanic activity. The bulk linear correlation coefficient suggests a best-fit value (m = VRP VIIRS /VRP MODIS ) equal to 0.87 (r = 0.95), with single volcanoes coefficients spanning from 0.78 to 1.14 ( Table 4). Based on this analysis we may conclude that the VIIRS detects from 32 to 108% more alerts than MODIS with an average of + 62% on the entire dataset. This percentage varies from volcano to volcano and, as previously mentioned, it is likely to be attributed to the improved features of the VIIRS instrument that allow detecting smaller thermal anomalies.
The above points confirm that the VIIRS instrument is fully consistent with MODIS and suitable for replacing and improving the latter in monitoring high-temperature volcanic heat emissions for the next decades. Moreover, our results suggest that the datasets processed by the MIROVA and FIRMS systems are comparable and provide consistent estimates on the radiant power sourced by volcanic eruptions. This cross-validation strengthens, even more, the value of the VRP which therefore constitutes a robust, multiplatform, multi-algorithm, indicator of volcanic thermal activity. On the other hand, the MIROVA algorithm, and in particular the MIROVA VIIRS version, turns out to be more efficient in detecting low-intensity thermal anomalies (<10 MW). This greater efficiency translates into a better ability to discriminate between background and anomaly values (cf. Figure 6) as well as to detect small signs of activity, otherwise undetected by the FIRMS. In terms of volcanic surveillance, this difference can be very significant at restless volcanoes as the first appearance of thermal activity can be detected days, weeks, or even months in advance compared to MODIS [7]. In this view, the application of the MIROVA algorithm to VIIRS data at 375 m resolution has all the potential to improve the detection capability of hot-spots during volcanic unrest phases as recently proved by [103].
Given the imminent end of the EOS missions, the future JPSS missions (the next scheduled for 2022, [104]) assume a key role. They will both fill the lack of MODIS data and potentially improve the monitoring capabilities with a higher number of daily acquisitions. On the other hand, our results show that there is margin for improvement in thermal remote sensing of volcanic activity, especially if future space missions, expressly dedicated to volcanic monitoring [46], will provide mid-infrared data with high spatial and high temporal resolutions in continuity with the VRP time series inherited from MODIS and VIIRS.