A Statistical Approach to Satellite Time Series Analysis to Detect Changes in Thermal Activities: The Vulcano Island 2021 Crisis

Abstract

Vulcano belongs to the seven volcanic islands forming the Aeolian archipelago (Italy) and has the privilege to define an eruptive style as “Vulcanian”. It has to be considered as an active volcano as its most recent activity demonstrated. Starting by late spring 2021, the thermal state of the Vulcano summit area changed and the gas emission increased. During the summer and, in particular, starting from September, geophysical and geochemical signals, precisely those linked to the activity of the hydrothermal system that feeds the fumaroles of the Fossa crater, varied. The temperature of the gases emitted by the fumaroles on the crater rim has increased and the composition of the gases has showed an increase in CO₂ and SO₂ (carbon dioxide and sulfur dioxide) concentration. For such reasons, the authors decided to follow this event by analyzing the remotely sensed available data suitable for detecting changes in thermal state. By processing the TIRS (Landsat 8) and ASTER time series, two long-term surface temperature logs were obtained and, therefore, by adopting a statistical approach, an analysis in both space and time domains has emphasized a thermal signature since mid-September 2021.

1. Introduction

The island of Vulcano is the southernmost emerged volcanic edifice in the Aeolian archipelago. Vulcano has erupted in historical times and the most recent activity occurred between August 1888 and March 1890 [1]. Given the small size of the island (about 20 km²) and of the volcano-related active phenomena (fumaroles, hot muds, etc.), the observability by space represents an interesting challenge. Because of its characteristics, a huge volume of scientific literature has been written focusing on the different aspects: concerning the geology and the eruptive dynamics of Vulcano, only in recent years [2,3,4,5]; the gas dispersal, lava flows, pyroclastic density currents, ballistic blocks, and tephra accumulation were deeply described by [2,6,7,8,9,10]. The remote sensing community has demonstrated that the use of Thermal Infrared (TIR) satellite sensors represents a reliable method [11] to identify volcanic activity [12]. Despite the monitoring activity that could be performed only by using low-spatial-resolution and high-revisit-time sensors, such as Meteosat Second Generation (MSG) [13,14], Moderate Resolution Imaging Spectroradiometer (MODIS), or Sea and Land Surface Temperature Radiometer (SLSTR) on board Sentinel 3 (with 1 km ground resolution) [15,16,17], sensors with higher spatial resolution (but lower revisit time) allow for the analysis of smaller phenomena, both in terms of spatial extension and magnitude of the event [18,19,20].

In this work, an Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) sensor embarked on Terra and Thermal Infrared Sensor (TIRS), embarked on Landsat 8 (L8), time series are processed in order to verify if the surface hydrothermal crisis on the Island of Vulcano, occurring since summer 2021, can be detectable by space-based TIR sensors, contributing to the identification of possible changes in volcanic activity and following the variation trends. This identification is allowed by the availability, in the spectral range of the TIR (8–12 microns) of data acquired with Ground Sampling Distance (GSD, more generally, spatial resolution) of 90–100 m, as only TIRS L8 and ASTER have.

Starting from September 2021, the monitoring systems managed by Istituto Nazionale di Geofisica e Vulcanologia (INGV) recorded significant variation in micro-seismicity linked to the circulation of fluids and soil deformations in the La Fossa crater area [21]. Moreover, significant changes in the composition and temperature of the fumaroles on the crater rim and an increase in degassing diffused from the soil in the areas adjacent to the crater [22] were also measured. As a consequence of this unrest phase, we decided to investigate the capability of a space-based TIR sensor to support the surveillance activities on the island of Vulcano, with particular regard to surface thermal state, related to the increase in temperature in the fumarolic fields.

In [23], the authors built a complete heat loss model for a hydrothermal system feeding fumarolic emission from the Vulcano La Fossa crater, by combining ground and ASTER data and constraining the thermal flux of the Vulcano La Fossa hydrothermal system between 2000 and 2019.

In [20], we analyzed the possibility to use satellite time series to identify thermal anomalies on Vulcano, observing 2000–2017 spaceborne remotely sensed time series data of surface temperature covering the whole island and detecting hot pixels as focal points.

In this work, we enlarged the time series up to February 2022, with the aim to support the operations of monitoring and detection of changes in volcanic activity occurring in the summer 2021.

2. Study Area

The Volcano Island is one of the seven Aeolian Islands, located in the north-eastern part of Sicily (Italy) (Figure 1a). It is characterized by two main thermal areas that can be easily detected and monitored by high-spatial-resolution satellite data (TIRS and ASTER): Fangaia, located near to the Levante harbor and La Fossa, in the central part of Vulcano Island [22] (Figure 1b).

These thermal areas result from the circulation of water vapor at high temperature coupled with other gases related to volcanic activities, such as CO₂ and SO₂ [3].

In this work, we focused on the La Fossa cone site, where the fumarolic activity is highly localized (Figure 2) and it was recently characterized by a significant increase in activity. For this reason, the Italian Civil Protection Department, in line with what is also expected at the international level, changed the alert level of Vulcano Island from green to yellow [24]. That means the volcano is exhibiting signs of elevated unrest above known background activity.

Before this current one, there were already periods of increased intensity in degassing activities, particularly in 1978–1980, 1988–1991, 1996, 2004–2007, and 2009–2010 [25]. These crises were related to an unrest phase, defined as variation from the standard activity [26] and were usually recognized by monitoring geophysical and geochemical parameters.

Starting from September 2021, the local seismicity, consisting of seismic events of low energy recorded in the area of the La Fossa crater and linked to the dynamics of the circulating fluids in the superficial hydrothermal system of Vulcano, showed a significant increase in daily frequency and a positive variation in the gradient of the cumulative number, with a maximum of 59 events observed on 27 September 2021 (Figure 3) [22].

In the same period, geochemical parameters in terms of temperature of fumaroles on the crater, CO₂ and SO₂ flux in the crater area also increased (Figure 4A,B) [22].

3. Materials and Methods

The spatial scale at which thermal phenomena are observed is given by the GSD and swath, while the temporal scale is linked to the revisit time on a given area, which in turn depends on orbit characteristics and swath dimension and steering capability in the sensor. In order to observe phenomena on a “regional to local” scale, GSD of 60–100 m is required as well as a short revisit time. As of today, ASTER [27] on Terra satellite and TIRS on L8 [19,28] offer such a capability with 90–100 m of spatial resolution (swath width of 60 km and 185 km, respectively, within the 8–12 µm atmospheric window) and a revisit time of 16 days.

Table 1. ASTER and TIRS bands features; * TIRS bands are acquired at 100 m spatial resolution but are resampled to 30 m in delivered data product.

Description	ASTER	TIRS	Unit
Number of thermal bands	5	2	-
Bands centers	8.29, 8.63, 9.07, 10.66, 11.32	10.9 (Band 10), 12.0 (Band 11)	µm
FWHM of bands	0.35, 0.35, 0.35, 0.7, 0.7	0.6, 1.0	µm
Pixel size at nadir	90	100 *	M
Revisit time	16	16	Day
Swath width	60	185	Km

reports the main features of both instruments.

Thermal anomalies were studied analyzing the Land Surface Temperature (LST) estimated by ASTER and TIRS thermal channels. For retrieving LST, the Temperature and Emissivity Separation (TES) methodology [29] was applied to ASTER thermal channels, while for TIRS thermal data, due to the presence of only two bands in TIR region and considering the USGS announced caution in the use of Band 11 of TIRS due to the calibration uncertainties [30], we applied the Single-Channel Algorithm (SCA) described and tested in [31,32,33]. To match the data in terms of spatial resolution, the LST retrieved by TIRS data was resampled to GSD of ASTER data (90 m) by using the nearest neighbor method.

Cloud-free and nighttime satellite data were freely downloaded by the EarthExplorer platform [34] managed by the USGS. The choice to use nighttime data is due to the necessity to minimize the effect of solar heating in temperature estimation and to magnify the differences in temperature between the anomalies and the background. Furthermore, during the nighttime, the topographic effect is minimized and the reflectance component is negligible. Another advantage of nighttime images is that surface radiant temperature is more constant than in the daytime. For these reasons, ASTER and TIRS nighttime image series were considered in this study. A final dataset consisting of more than 300 images was collected (

Table 2. Satellite images considered for the analysis.

Sensor	Time	Number of Processed Images
ASTER	6 June 2000–4 February 2022	232
TIRS	4 December 2013–2 February 2022	76
Total	6 June 2000–4 February 2022	308

).

The availability of a substantial number of these satellite data for the Island of Vulcano offered the possibility of obtaining LST time series over a long period, thus, allowing for an accurate analysis of thermal anomalies. Due to the different methodology to derive the LST, the results were sorted according to the type of satellites and for such reason, two different LST time series were generated for ASTER and TIRS, respectively.

A first investigation was carried out by analyzing the time distribution of the surface temperature value over a specific area, covered by a single satellite pixel, using for both sensor the 90 m GSD. In both TIRS and ASTER time series, the hottest pixel was chosen by performing Pixel Purity Index (PPI) analysis. PPI analysis allows to find the extreme pixel in multispectral and hyperspectral images, which represents the most spectrally pure pixel in the image. In this work, we apply, unconventionally, this analysis on LST time-series maps to discriminate which are the pixels most frequently marked as those hottest (Figure 5). By following this approach, we identified 4 pixels that are the hottest in 284 cases on 308 total equal to 92% of the occurrence. By considering these four pixels as a loci of the thermal anomaly we depicted a 10 × 10 pixels Region Of Interest (ROI), centered on the pixels detected by the PPI analysis; therefore, the statistical analysis is applied taking into account both the area affected from the fumarolic activity (red pixels in Figure 5) and a background area (green pixels in Figure 5) that is barely affected by the activity identified as Background Reference Pixel (BRP).

One important effect that is visible working with LST time series is seasonality. This effect illustrates annual insolation variations and is typically modeled as a sinusoidal function (Figure 6). To better discriminate the intensity of the thermal anomaly, the seasonality effect has to be defined analyzing an area far from those considered thermally active. We focused our attention on two main areas; the BRP far from La Fossa crater (Figure 5 and the four red pixels, derived from PPI analysis, in La Fossa crater rim, Figure 5). Both areas analyzed are affected by the same seasonality effects (Figure 6). One of these effects has yearly wavelength and it could be modeled as a sinusoidal function (Figure 6). In order to avoid the seasonality of the data and to highlight the hydrothermal activity of La Fossa, we considered a statistical approach based on standard deviation.

The statistical analysis consists of calculating the standard deviation (σ_map) of temperatures for each LST map and studying the trend of this parameter over time. The standard deviation describes the spread of data around the mean value: low standard deviation means that temperature values are clustered around the mean value while high standard deviation means that temperature values are widely scattered around the mean value [35].

As reported in Table 2, the total dataset is composed of 308 samples covering the period from 2002 up to 2022; each sample represents the solution of the methods applied for the generation of LST map (see Section 3).

In our study, we computed the standard deviation for the EO data processed and the source of such variability was investigated.

From each LST map the parameter σ_map is computed, within the same spatial extension of the ROI (10 × 10 pixel, Figure 5), in order to obtain a historical time series. The behavior of σ_map over time is analyzed to detect extreme values, likely related to changes in thermal activity. With the aim to automatically find high values of σ_map, mean (M_ASTER, M_TIRS) and standard deviation values (σ_ASTER, σ_TIRS) of the two satellite time series are calculated. In the algorithm presented here, a value of σ_map is defined outlier if greater than M_ASTER + 2σ_ASTER and M_TIRS + 2σ_TIRS, respectively, for ASTER and TIRS datasets. Considering the Gaussian distribution of σ_map samples (Figure 7), the choice of 2σ_ASTER and 2σ_TIRS values as limits allows us to consider about 95% of the statistical population as regular and about 5% as outlier.

Figure 8 reports the scheme of the algorithm proposed in this paper.

4. Results and Discussion

In this section, the results will be shown and discussed. A first consideration will be carried out about the PPI surface temperature value trend over the time. Then, some consideration about the results derived from the statistical analysis will be pointed out, once specifying the size of the LST maps equal to the ROI.

4.1. Surface Temperature Variation vs. Time

Since ASTER and TIRS were co-registered and resized to the same pixel spatial resolution (90 m), a complete temperature time series (2002–2022) was obtained (Figure 9). In Figure 9, the plot representing the cumulative (ASTER + TIRS) LST time series, focused on the 2013–2022 period, is evidenced. The choice to set the lower time limit set on 2013 depends on the availability of TIRS data. The quasi-sinusoidal trend related to the seasonal behaviour is easily observable. Indeed, high values of temperature in summertime and low values of temperature in wintertime are evident. Therefore, concurrent with the increase in Vulcano thermal activity, a rise in temperature value is evident.

4.2. LST Map Standard Deviation over Time

The statistical approach considered in this paper allows us to neglect the seasonal effects and to highlight the role of the outliers. It is important to remark that TIRS and ASTER do not acquire on the same days. For this reason, the number and the date of the recognised outliers could be relevant to different sources.

Considering the TIRS series, in Figure 10, four outliers are pointed out.

The TIRS time series (Figure 10) shows that the 15 August 2019, 21 September 2021, 27 January, and 28 February 2022 are outliers. The first outlier (1 in Figure 10) exceeds the 2σ value and it depends on the presence of cloud coverage in the background area, which impacts the value of the standard deviation (Figure 10).

The other outliers (2, 3, and 4 in Figure 10) cover a quasi-continuous time period ranging from September 2021 up to the last data analyzed (4 February 2022) and their trend is in agreement with the INGV bulletins [21,22]. This permits these outliers to be assigned to the current crisis on Vulcano island. Whilst the TIRS data relevant to two and four outliers were acquired under optimal weather conditions, number 3 was characterized by cloud coverage of the same magnitude as number 1. However, the following considerations can be made: outlier number 1 is relative to data acquired during summer and, therefore, the low temperatures to the clouds correspond to a high environmental temperature due to the hot season. Outlier number 3 is relevant to EO data sensed during winter (27 January 2022) and its value depends on the presence of clouds that determines the lower temperature and to the hottest that cannot be referred to the environment (winter see Figure 9) but to the contemporary crisis.

By analyzing the ASTER time series (Figure 11) in the same TIRS time window (data from 2013), the three detected outliers match in time with the current crisis of Vulcano and they are consistent with the results derived from TIRS analysis, considering the good weather conditions present at the acquisition dates.

Considering the previous discussed results and joining the two datasets based on the temporal range with the available ASTER data, we obtain the largest dataset, which allows us to increase the temporal resolution in the time series and extend the thermal activity analysis of Vulcano, starting this back analysis from 2000 up to now and to recognize the anomalous events occurring in 2004–2007 and 2009–2010, as reported in [25,36] (Figure 12).

Comparing our results to 2004–2005, despite these points not being identified as outliers, they have high standard deviation value (about 2.5), as their value is close to +2σ, suggesting a relation with the occurred crises reported in [25]. For the 2009–2010 event, also described in [25,36], we obtain an increasing trend up over 2010 and also to the beginning 2012. The others match with the current crisis in Vulcano.

The application of standard deviation permits one to neglect the absolute value of the surface temperature obtained by the two different methods (TES and SCA, respectively, for ASTER and TIRS). Despite the differences between the used orbiting system, as already verified in [33], the procedures adopted led to a generation of comparable datasets where the application of standard deviation method further reduces the uncertainties due to the different datasets.

5. Conclusions

Thermal data can have an important role in the delineation of surface temperature anomalies associated with thermal activity and for the measurement of near-surface heat fluxes associated with thermal systems.

The studied case represents demonstrations of the potentiality of satellite observations in TIR for environmental applications and confirms the request to the remote sensing community for a continuous improvement in satellite sensors in the TIR, especially in terms of spatial resolution.

The use of two sensors (ASTER and TIRS), despite the differences in processing method applied for LST estimation and GSD, produced results with high correlation, offering the possibility to extend, in time, the LST time series. The use of two sensors proved to be a cost-effective technique for generating products that allow for not only the detection but also the monitoring of thermal anomalies, permitting one to analyze how the anomaly evolves over time and to highlight critical changes in thermal activities.

Despite the methodology used, this work is based on a very simple use of statistical tools and this approach permits one to highlight the outliers, which correspond to the evidence of the surface thermal state on Vulcano island.

With this work, a robust and easy-to-use procedure for the detection of thermal anomalies and their activity was applied. The proposed methodology can be easily replicated, ensuring the conceptual correctness and reproducibility of image processing.

This study will be updated by using data that will be provided by the new planned and orbiting mission with a GSD finer than what we already use (90 m), on which we could rely for the coming years.