Anatomy of thermal unrest at a hydrothermal system: case study of the 2021–2022 crisis at Vulcano

Hydrothermal systems can generate phreatic and/or phreatomagmatic explosions with little warning. Understanding the temporal and spatial evolution of geophysical and geochemical signals at hydrothermal systems is crucial for detecting precursory signs to unrest and to inform on hazard. Thermal signatures of such systems are poorly defined because data records are often too short or discrete compared to activity timescales, which can be decadal. La Fossa system of Vulcano has been monitored since the 1980s and entered a period of unrest in 2021. We assessed the thermal signature of La Fossa using ground- and satellite-based data with various temporal and spatial scales. While continuously-recording stations provided continuous but point-based measurements, fumarole field vent surveys and infrared images obtained from satellite-flown sensors (ASTER and VIIRS) allowed lower temporal resolution but synoptic records to be built. By integrating this multi-resolution data set, precursory signs of unrest could retrospectively be detected from February to June 2021. Intensity of all unrest metrics increased during summer 2021, with an onset over a few days in September 2021. By September, seismic, CO2, SO2 and other geochemical metrics also indicated unrest, leading Civil Protection to raise the alert level to yellow on October 1. Heat flux, having been 4 MW in May 2019, increasing to 90 MW by September, and peaking at 120 MW in March 2022. We convolved our thermal data sets with all other monitoring data to validate a Vulcano Fossa Unrest Index (VFUI), framework of which can be potentially applied to any hydrothermal system. The VFUI highlighted four stages of unrest, none of which were clear in any single data set: background, precursory, onset and unrest. Onset was characterized by sudden release of fluids, likely caused by failure of sealed zones that had become pressurized during the precursory phase that began possibly as early as February 2021. Unrest has been ongoing for more than 18 months, and may continue for several more years. Our understanding of this system behavior has been due to hindsight, but demonstrates how multiparametric surveys can track and forecast unrest.

Whether unrest leads to explosive activity or not, geophysical signals recorded for restless systems can provide insights into physical and dynamic processes in operation, in and above the hydrothermal system (Montalto 1994;Mannen et al. 2018;Moretti et al. 2020;Girona et al. 2021).In addition, even unrest that does not lead to explosive activity presents hazards as sometimes testified by high contents of toxic gases in the ambient air, as for example at Mammoth Mountain (USA) during 1989 (Sorey et al. 2000).
However, some metrics or signs of unrest are poorly defined, and this includes the thermal signature.This is mostly because temperature measurements tend to be punctual and so define the thermal character of the system only for a particular point in a cycle of unrest, and lack the temporal resolution and duration to capture a complete cycle (Aubert 1999;Matsushima et al. 2003).As a result, continuous, high-temporal resolution thermal data sets recording the onset of unrest are lacking.
The ground temperature for points within the fumarole field of Vulcano's La Fossa system (Aeolian Islands, Italy) have been monitored since 1984 (Diliberto 2011).In addition, surveys have been repeated annually since 1995 to detect changes of vent temperatures across the field (Harris and Maciejewski 2000;Harris et al. 2012), and ASTER and Landsat Thematic Mapper imagery have been used to constrain the thermal anomaly associated with the fumarolic activity and derive heat flux (Gaonac'h et al. 1994;Harris and Stevenson 1997;Silvestri et al. 2019;Mannini et al. 2019).Since January 2020 onwards, we have also been running up to five continuously recording temperature sensor stations in and around the fumarole field.In parallel, we have completed fumarole temperature up to three times a year, coupled with infrared camera imaging.Over the summer of 2021 Vulcano's hydrothermal system entered a new phase of unrest, which allowed us to record and define the thermal signature of unrest.The northern hemisphere, Mediterranean summer generally runs from June through August, and the winter from December through February; where the summer solstice in 2021 was on 21 June, and the winter solstice was on 21 December.
We thus use the unrest that began in the second half of 2021 to define its signature in different thermal data sets, and examine how the temperature trends couple with other geochemical and geophysical data sets.This allows us to develop a model whereby unrest is preceeded by three-to-four months of slowly increasing heat flux, before an abrupt increase at the onset of unrest, with heat flux increasing by 10 s of MW over just one or two days.Pailot-Bonnétat et al. Earth, Planets and Space (2023)

Geological setting and context
Vulcano is located in the southern part of the Eolian Archipelago (Fig. 1A) and is the emerged summit of a volcanic edifice that extends c. 1000 m b.s.l.(below sea level) (De Astis et al. 2013).It is a composite volcanic system built during a series of eruptive phases, with first activity being dated at 130 ka (Keller 1980;De Astis et al. 2013).Structures and vents are aligned NNW-SSE, an alignment which is controlled by the Tindari-Letojanni strike-slip fault system that passes beneath the island (Ventura 1994;Mazzuoli et al. 1995;De Astis et al. 2013).The associated NE-SW extensional stress field is linked to the two major calderas on the island, the Piano and the Fossa, which have been interpreted as pull-apart basins (Ventura 1994;Mazzuoli et al. 1995), as well as the alignment of zones of fumarolic activity.

Thermal signature associated with fluid ascent at the Vulcano hydrothermal system
The Fossa, Faraglione and Vulcanello eruptive centers (Fig. 1B) have been active since 5.5 ka (De Astis et al. 2013).The last eruption occurred in AD 1888-1890, and its description by Mercalli and Silvestri (1891) became the archetype for Vulcanian eruptions (Mercalli 1907).Mercalli (1883) noted that emissions during late 1879 were sufficiently high that access to the crater for mining became impossible.De Fiore (1922) noted the establishment of a new fumarole field (maximum temperature 110 °C) in the Fossa crater by 1913, andSicardi (1940) recorded increasing then decreasing temperatures between 1923 and 1937.The maximum recorded was 615 °C in 1924(De Fiore 1924).
Fumarolic activity is mainly located in La Fossa crater and on the Spiaggia di Levante beach (Tedesco et al. 1995).While Ferrucci et al. (1991) found a magma body at 2-3 km depth, the magmatic contribution to the fumarolic gases has been inferred from chemical composition, helium isotope ratios and 3He output (Nuccio et al. 1999).Many authors have interpretated the gas composition at Vulcano as resulting from mixing of magmatic and deep hydrothermal fluids (marine and meteorological sources), with the hydrothermal system being a biphase (water-vapor) boiling saline solution (e.g., Carapezza et al. 1981;Chiodini and Cioni 1995;Nuccio et al. 1999).Carapezza et al. (1981) calculated that this saline solution was boiling at a temperature of 330-340 °C at a pressure of 13-17 MPa; Nuccio et al. (1999) assumed these P-T conditions in late 1970's, whereas P-T values around 400 °C and 20 MPa were hypothesized during the 1990s.
Mixing models have shown that unrest is associated with increased CO 2 and He fluxes, representing gas liberated by vesiculated magma (Nuccio et al. 1999).
Background values of these gases are interpreted as being derived mainly from diffuse degassing of volatile impoverished melt (Nuccio et al. 1999).The unrest of 1988-96 was characterized by homogenization of the mixing fractions.This effect was interpreted as vaporization due to boiling, transforming the system from biphase to vapormonophase; able to supply a uniform composition to all fumaroles (Nuccio et al. 1999).Fumarole temperatures also reached a maximum of 690 °C (Capasso et al. 1994;Chiodini and Cioni 1995), the area of fumarolic emission increased (Bukumirovic et al. 1997;Italiano et al. 1998), and the isotopic signature evolved toward a clear magmatic component (Tedesco et al. 1991;Capasso et al. 1997).
Since the late 1990s, fumarole temperatures have generally decreased (Diliberto 2011(Diliberto , 2017;;Harris et al. 2012), and the biphase system has re-established.Gas expansion during ascent, and the subsequent gas temperature on emission is, however, highly dependent on rock permeability over the last few hundred meters (Nuccio et al. 1999).Changes in shallow permeability can trigger heating phases (increases in fumarole temperature), where increases in permeability associated with fracturing of sealed zones allow more efficient heat flux without new magmatic input (Diliberto 2011(Diliberto , 2017;;Harris et al. 2012).Increased permeability has also been shown to relate to deflation, of a source at sea level, due to enhanced steam emission and a consequent reduction in volume of water in the hydrothermal system (Gambino and Guglielmino 2008;Harris et al. 2012).

The new (2021-present) phase of unrest
After around two decades of declining levels of thermal activity (Harris et al. 2012;Mannini et al. 2019), a new period of unrest was registered during the summer of 2021.CO 2 , SO 2 fluxes and temperature data from the INGV monitoring network revealed slowly increasing emissions starting in June 2021, with a sharp increase in September 2021 (Inguaggiato et al. 2022a).The number of Very Long Period (VLP) seismic events also increased in September 2021, and was accompanied by 2 cm of inflation, indicating fluid pressurization within fractures (Federico et al. 2023).The temperature of the fumaroles increased, but the maximum temperature (372 °C) remained around the critical temperature of the Fossa crater boiling saline solution (~ 400 °C) and was very localized (Federico et al. 2023).This meant that, in contrast to the 1988-96 unrest, the hydrothermal system underwent only partial vaporization, so the system did not completely dry out and the biphase state was maintained (Federico et al. 2023).
Re-mobilization of sulfur stored in the system was also observed in direct measurements of SO 2 flux and plume chemistry (Aiuppa et al. 2022).It also resulted in the eruption and deposition of sulfur (Additional file 1: Appendix A).Sulfur melting and emission of molten sulfur is not rare during unrest at Vulcano, where sulfur flow activity was common during the 1988-96 unrest (Harris et al. 2000(Harris et al. , 2004)).Sulfur melts between 112 and 120 °C (Meyer 1976), so heating of sulfur deposits in the shallow system above this temperature range during unrest will cause them to fluidize and erupt.However, between 1998 and 2020, no new eruptive sulfur structures, such as the hornitos, flows or air fall were observed, and the degree of sulfur deposition was markedly reduced.

Thermal character of the fumarole field
Heat flux at a fumarole field is partitioned between diffuse soil emissions and emission at fumarole vents (Sekioka and Yuhara 1974;Chiodini et al. 2005;Harris 2013).Unlike at some systems, e.g., La Soufriere (Jessop et al. 2021), the heat flux at Vulcano has been proven to be dominated by the diffuse component, as shown by independent studies based on the water flux, satellite remote sensing, thermal camera, and fumarole temperature measurements (Chiodini et al. 2005;Harris et al. 2009;Mannini et al. 2019).This is also the case at Nisyros (Greece) and Campi Flegrei's Solfatara (Chiodini et al. 2005).At such systems, the source of the heat has been ascribed to condensation of ascending fluids (Chiodini et al. 2005), where heat is then transferred to the surface by permeable convection and, over the last meter, conduction (Aubert 1999).Heat is then lost from the surface by radiation and convection (Sekioka and Yuhara 1974;Harris 2013).
At Vulcano, the zone of fumaroles is thus surrounded within a broader zone of soil emission.Heat flux from the zone of soil emission dominates the energy budget, accounting for 93 ± 2% of the total budget, where total heat fluxes were in the range 5-13 MW between 2000 and 2019.Reducing conditions within the soil emission zone changed the surface color from red to gray (Fig. 1D).We here label the former zone hot (grey) zone, and the latter cold (red).However, although accounting for less than 8% of the total flux, changes in the thermal character (number and temperature) of the fumaroles provides information on whether the system is heating or cooling (Harris and Maciejewski 2000;Matsushima et al. 2003), as well as on localized changes in shallow system permeability (Harris et al. 2012).

Methods
Between January 2020 and January 2022, we carried out six field campaigns, spaced by three to eight months.We carried out thermal measurements in January and October 2020, March, June and September 2021, and January 2022.This sampling frequency was limited by travel restrictions and lockdowns due to the Covid pandemic.The full list of measurements and their locations as well as instruments used is given in Additional file 1: Appendix B.

Surface temperature measurements
Two temperature stations separated by 50 m were installed in January 2020 at the west edge of the fumarole field: RED_station and GRAY_station (Fig. 1C).Each station records temperature at 15 cm depth (T 15 ), surface temperature and air temperature 5 cm above ground level, every 5 or 10 min.Two other temperature stations (MZ1 and MZ2) were installed inside the hot zone at the end of June 2021 (Fig. 1C).In addition, we used data from the permanent station maintained by INGV-PA (DDS station, Fig. 1C) on the east side of the fumarole field and HTF (high temperature fumarole) FA, F5 and F5AT (Fig. 1C).DDS station provides temperature at 15, 30, 45 and 60 cm depths.During each campaign, T 15 and surface temperature were measured every 5 m along a 50 m-long profile (Profile1, Fig. 1C) between RED_station and GRAY_station (Fig. 1D), every 1 m within a 5 × 5 m grid (Grid1, Fig. 1C), and every 15 m on 250 m and 300 m long profiles traversing the entire field (Profiles 2, 3 and 4, Fig. 1C).In March and June 2021, we also made T 15 measurements for every 5 m across a 110 × 50 m grid (FA_grid, Fig. 1C).Once the 2021 unrest began, we could not resume the FA_grid measurements due to dangerous conditions.However, simultaneous measurements at Grid1 in March and June allowed us to assess the representativeness of the 5 × 5 m area.Detailed descriptions of the temperature measurements at FA_grid are given in Additional file 1: Appendix B.

Soil CO 2 flux measurements
Soil CO 2 flux (g m −2 d −1 ) was measured using the accumulation chamber method of Chiodini et al. (1998).Data was processed with the "Flux Revision" software (https:// www.wests ystems.com/ downl oad/) using the raw soil CO 2 content regression line measured by the instrument and ambient pressure-temperature (acquired from a Kestrel 5500 hand held weather station) to obtain the flux (ppm s −1 ).Conversion in (g m −2 d −1 ) is obtained by considering the volume and area of the chamber together with air temperature and pressure.

Fumarole surveys
Once a year since 1994, and two to three times a year since 2020, a survey using Thermal InfraRed (TIR) thermometers was conducted to measure vent temperatures.The survey was always conducted following the same path and pattern to assure its statistical consistency.On the basis of historical evolution and thermal characteristics, the fumarole field was divided into five zones by Harris and Maciejewski (2000).These are: (1) the rim rifts which are two fissures that cut across the rim on the west edge of the field, (2) fumaroles in fractures aligned along the rim, the (3) upper, (4) middle and (5) lower zones (Fig. 1D).

Thermal images
Night time thermal images were acquired coincident with scheduled satellite overpasses (Ramsey 2016).The camera was placed opposite the fumarole field, on south rim of the crater (Fig. 1C).This position allows a synoptic view of the five fumarole zones from the rim to the bottom of the crater.For a line-of-sight distance of 310 m to the bottom of the crater and 460 m to the rim, the pixel size ranged from 0.4 m to 0.6 m.To image the entire field, a panel of three images was taken as a panorama.

ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer
) is a sensor flown on the TERRA NASA's low earth orbit satellite, which measures spectral radiance at 14 wavebands between the visible and thermal infrared.Given that we are interested in kinetic surface temperature, we here use data from the five Thermal Infrared (TIR) bands between 8 and 12 µm.Specifically, we used the AST_08 image product with 90-m pixel resolution, in which pixel temperatures are automatically corrected for emissivity and atmospheric effects.Nighttime scenes were selected to avoid solar heating and thermal inertia effects.Masking allowed extraction of temperature data for each pixel in two regions of interest (ROI): the Fossa crater, where the anomaly is visible, and Vulcanello, a cone north of the Fossa which currently has no anomalous heat and gas release (Fig. 1B).Maximum and mean temperatures from the ambient control (Vulcanello) were compared to those of the active target (i.e., the Fossa).
Within the Fossa ROI, the number of anomalous pixels was calculated using the temperature distribution (Additional file 1: Appendix D).Almost all distributions are bimodal: the first mode being the background temperature distribution and the second the anomaly (Fig. 2a).If heating occurs, the hot part of the distribution becomes more robust making it more skewed to hot temperatures.Pixel location in terms of ground target and viewing angles are variable from overpass to overpass, which might affect the clarity of the bimodality.However, the pixels detected in this study are coherent with thermal anomaly location and evolution found in situ.The number of pixels within the second mode was extracted and multiplied by the pixel area to obtain the thermal anomaly area (Fig. 2b; Additional file 1: Appendix D).This method was validated by cross-comparison with an independent deep learning approach (supervised UNET network) (Corradino et al. 2023).

Conversion to heat flux
Following Sekioka and Yuhara (1974) and Matsushima et al. (2003), radiative (Mrad) and convective (Mconv) heat fluxes were calculated on pixel-by-pixel basis in ASTER and thermal camera images, as follows: where σ is the Stefan-Boltzmann constant, ε is emissivity, h c is the convective heat transfer coefficient, T a is the temperature of the anomalous pixel, and T b is background temperature.Multiplying by pixel area converts these values from flux densities (M in W m −2 ) to heat fluxes (Φ in J s −1 ).For emissivity, we used the mean of our measured values (see Additional file 1: Appendix B) and for h c , we used 24 J K −1 after Sekioka and Yuhara (1974) and Matsushima et al. (2003).For the background temperatures, we used the minimum and mean temperatures of non-anomalous pixels.Values were summed for all pixels to obtain total Φrad and Φconv and the totals were added to obtain total heat flux [Φtot = Φrad + Φconv].

VIIRS processing
VIIRS (Visible Infrared Imaging Radiometer) is a sensor flown on the NASA/NOAA Suomi National Polarorbiting Partnership and NOAA-20 low earth orbit satellites.VIIRS measures spectral radiance at 22 wavebands between the visible and thermal infrared, where we used data from two of the mid-infrared band I4 (3.55-3.93µm) and the TIR band I5 (10.5-12.4µm) which have pixel sizes of 375 m.
VIIRS data were analyzed using a specialized version of CL-HOTSAT (Ganci et al. 2016).CL-HOTSAT is a satellite data processing system used both for the nearreal-time monitoring of high-temperature volcanic features (Rogic et al. 2019;Ganci et al. 2023) and for the extraction of the input parameters of numerical models used for lava flow hazard and risk assessment (Cappello et al. 2016;Zuccarello et al. 2022).The VIIRS data were downloaded as Level 1 Sensor Data Record (SDR) format including L1B calibrated radiance products VNP02IMG and VJ102IMG, for SNPP and JPSS-1/NOAA20 platforms respectively, and VNP03IMG and VJ103IMG VIIRS L1 terrain-corrected geolocation products (1) (2) containing the derived line-of-sight (LOS) vectors for each of the 375-m image-resolution or I-bands.The geolocation product was used for reprojecting and georeferencing the calibrated radiances that were then corrected for the atmospheric effect.Nighttime images were selected to avoid possible effects due to solar reflection.Moreover, in order to locate the thermally anomalous pixels a new hotspot detection algorithm was developed based on the Normalized Thermal Index (NTI) computed using the I4 and I5 bands (Wright et al. 2002).
This hotspot algorithm is a contextual algorithm in CL-HOTSAT where the volcanic area includes the whole La Fossa crater area and the non-volcanic area the rest of the island.We flag a pixel in the volcanic area as a hotspot if its value for NTI (NTIHS) is greater than the normal variation of NTI in a non-volcanic area (NTINVA) according to: For each pixel flagged as a hotspot the radiant heat flux is computed by assuming the pixel as composed by two thermal components.A dual band technique is then applied to the hotspot pixel in order to infer their thermal structures and the Stefan-Boltzmann law is applied to compute the associated radiant heat flux.

Punctual measurements of temperature and soil CO 2 flux
Since 2020, the location of the transition between the cold and hot zones on Profile1 has been stable and located at the 20 m mark (Fig. 3).Between September 2015 and March 2021, T 15 for the cold zone was stable, the temperature range for any given measurement time varying by no more than 3 °C (Additional file 1: Appendix (4) NTIHS > Mean NTINVA + 3 * SD NTINVA .F). T 15 for the hot zone were up to 20 °C higher than the cold zone, with hot spots at 80-90 °C.These hotspots were highly localized, but spatially and temporarily variable.However, in September 2021, the hotspot expanded to include the entire hot segment of the line, along which temperatures were > 80 °C (Fig. 3).At the same time the hot-cold transition moved 5 m into the cold zone, with temperatures at the 20 m mark increasing from 25 °C in January 2020 to 90 °C in September 2021.By January 2022, the 20 m mark temperature had decreased by 40 °C, as had the extent of the hotspot.However, T 15 levels remained higher than before June 2021 (Fig. 3).This implies an expansion of the hot zone, as well as enhanced heating of the hot zone from June 2021 onwards, before stabilization at a new level associated with unrest.
In March 2021, T 15 for Grid1 was between 30 and 40 °C with a hotspot at 60 °C in the SE corner and a second at 46 °C in the NE corner (Fig. 4a).In June 2021, hotspots were located in the same two locations, but the SE hotspot had increased in temperature and the whole area had heated by 20 °C (Fig. 4b).By this time, temperatures at the hotspots had reached boiling point temperatures.In September 2021, almost the whole area was at boiling point (Fig. 4c).By January 2022, boiling point temperatures had focused on a NE-SW hot fissure along which temperatures were 50-60 °C (Fig. 4d).This zone was already apparent in the September grid (Fig. 4c) and persisted through May 2022 (Fig. 4e).The linear feature likely represents a line of preferential fluid circulation, such as a buried fracture.This NE-SW alignment is the direction of the rim rifts and all regional structures as controlled by the extensional stress field associated with the Tindari-Letojanni fault system (Ventura 1994;Mazzuoli et al. 1995).Temperature distributions across the FA_grid (Additional file 1: Appendix C, surface = 5500 m 2 ) showed similar behavior between March and June 2021.
Soil CO 2 flux was also measured along Profiles 1-4 and Grid1 (Fig. 1C).An increase in temperature between June and September 2021 was matched by an increase in soil CO 2 flux at all sites.At Grid1, for example, the average flux increased from 900 g m −2 d −1 in June 2021 to 4522 g m −2 d −1 by September.Locations of peaks in T 15 and soil CO 2 flux on Profile2 were correlated in June and September 2021 (Fig. 5), marking high permeability locations.However, the correlation declines where surface temperatures are closer to background, an effect of lower permeability so that soil CO 2 flux is almost zero.

Fumarole field survey
Vent temperatures show a general decline following the 1988-96 unrest, with lowest levels being recorded in March 2021 (Additional file 1: Appendix G).Thereafter there was an increase in the mean temperature which, in 14 months recovered to levels higher than in 1994.To clean possible anti-correlations (due for example to very low number of fumaroles measured, cf.Additional file 1: Appendix G) we developed a thermal index (FTI): where T max and T mean are the maximum and mean temperatures for each survey, SD is the standard deviation and Nb v is the number of vents measured.This index is designed to use T max as a proxy for the absolute thermal state of the fumarole field, T mean for the heat flux, Nb v for the size of the fumarole field, and standard deviation for the thermal behavior, where small standard deviation has been shown to be associated with cooling, and high standard deviation with heating (Harris and Maciejewski 2000).
The FTI shows a general decline through 2021, with a recovery to around the 1996 levels by June 2022 (Fig. 6).This approach minimizes the influence of changes in permeability.Increased fracturing and fluid flow have been shown to result in increased heat flux without any change in magmatic input, as occurred between 2002 and 2010 (Gambino and Guglielmino 2008; Alparone et al. 2010;Harris et al. 2012).It thus confirms the current trend as not related to such a process.We see the same thermal reactivation from June 2021 in all zones (Additional file 1: Appendix G).Reactivation was also apparent in the field, where we observed (1) activation of a fissure network with low temperature fumaroles outside of the main fumarole field area, (2) increased sulfur deposition, (3) increased high temperature sulfur mineralization (liquid orange sulfur and black sublimates), (4) eruption of sulfur flows and sulfur pyroclasts (Additional file 1: Appendix A), and (5) melting of the sulfur component of the substrate to leave a fine-grained grey silicate powder (Harris and Maciejewski 2000).

In-situ temperature in the hot zone
We assess the magnitude of heating using the average daily surface temperature minus air temperature at GRAY_station (= ΔT).In 2020 ΔT was negligible and stable (ΔT ± 1 °C, Fig. 7).However, ΔT increased to 2-3 °C in April 2021, reaching 8 °C by the end of 2021 (Fig. 7) as the system went into unrest.At depth of 15 cm, rainfall causes perturbations of ground temperature, where T 15 decreased by up to 20 °C during rainfall events, and taking one day to two weeks to recover (Fig. 8).However, around these events, there was a systematic increase in T 15 by 0.167 °C per day until August 23 at MZ1, at which point the rate increased to 0.7 °C per day through mid-September.Thereafter we record three pulses during which T 15 reached boiling temperature.The three pulses had durations of around half a month and were centered on late-September, late-October and late-November, the later pulses marked by discrete peaks with T 15 up to 110 °C.The time series collected at GRAY_station and, to a lesser extent, DDS mimic the behavior of MZ1, but with lower T 15 and less marked pulses and peaks (Fig. 8).This is consistent with DDS being further from the fumarole field, where the lower temperature and response at DDS suggests lower permeability conditions than at MZ1.

Thermal camera images
Thermal images were collected on 25 June and 29 September 2021 (Fig. 9a, b) during nighttime satellite overpasses.These images showed expansion of the thermal anomaly, especially in the west and south-west sectors of the crater, where the anomaly spread into previously cold zones.It also shows reactivation of areas within the Fig. 6 Long term thermal trend of the Fossa fumarole vents.Gaps are when the period of time between two surveys is greater than 1 year.The point for January 2022 was likely underestimated due to poor visibility (plume condensation) fumarole field, especially along the rim rifts, in the upper zone and areas of the middle zone.
Surface temperature distributions are corrected for air temperature to take into account seasonal cycles, meaning that the first bins of the temperature distributions are collocated (Fig. 9c, d).The skewness in June is higher because background temperatures dominate more than in September (Fig. 9c; Table 1).September's distribution also has a more robust high temperature tail, resulting in a lower skewness and highly positive kurtosis (Table 1), indicative of heating (Harris et al. 2009).The September distribution can be separated into three populations: (1) background temperatures (normal distribution), (2) heated ground (normal distribution), and (3) pixels containing fumaroles and subsequent heated ground (skewed distribution) (Fig. 9d).

ASTER images
There was no correlation between the temporal variation of mean or maximum pixel temperatures and the number of anomalous pixels (Additional file 1: Appendix H).However, there was plainly an increase in the size and magnitude of the anomaly (Additional file 1: Appendix D, E).We thus developed a method to highlight the thermal anomaly and extract any temporal trends using the number of anomalous pixels to take into account changes in size of the anomaly (Fig. 10).The resulting Satellite Thermal Index (STI) is given by:   where Tmean crater is the average temperature from the anomaly, Tmean control is the average temperature of the ambient control area on Vulcanello, SD crater is the stand- ard deviation for the anomalous temperature distribution, SD control is the standard deviation on Vulcanello and Nb p is the number of anomalous pixels.The role of the numerator is similar to that of Eq. ( 5) and minimizes the influence of seasonal and diurnal cycles.The index values are underestimated in case of partial cloud cover or rain fall (such as the mid-July 2021 storm event) (Fig. 10), but we see an increase possibly beginning in mid-May, to reach a first peak at the end of September 2021 (which might be underestimated because of the partial cloudiness).This is consistent with the increased number of anomalous pixels from mid-May.From September onwards the STI remained at a high level (Fig. 15).Note that trends in STI (Fig. 10) mimics those in the in-situ temperature data, providing ground truth validation.
The VIIRS-derived radiative heat fluxes are more scattered due to the greater pixel size (375 m) and subpixel cloud contamination, but are generally within range of ASTER-derived fluxes (Fig. 11a).
Thermal camera only targets the fumarole field on the NE inner flank, whereas ASTER and VIIRS target the cone-wide heat flux (the entire crater and all flanks).Thermal camera derived radiative heat fluxes are thus lower, being 1-2 MW (Fig. 11a), and are consistent with the VIIRS-derived radiative heat flux obtained for the same portion of the system by Coppola et al. (2022), i.e., 1.2 MW.This means that, in terms of radiative heat flux, the fumarole field contributes to 5-10% of the cone-wide flux, which is consistent with Mannini et al. (2019).
In terms of total heat flux, ASTER recorded 30-60 MW until October 2021, with a possible ramp up beginning in June.Total heat flux then increased abruptly to 80-90 MW during October 2021, with a peak at 120 MW in March 2022 (Fig. 11b).These values are higher than the flux of 21 MW calculated on the basis of steam flux in 1998 by Chiodini et al. (2005), which was at the end of the previous unrest.These total heat fluxes are also higher than the pre-unrest (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) baseline of 5-14 MW as given by Mannini et al. (2019) and as obtained from the same data (ASTER) using the same approach as applied here.The change from baseline fluxes is a result of the both the size (number of pixels) and magnitude Total heat fluxes from thermal imagery are 18 and 32 MW (Fig. 11b), being higher in September than July.This is consistent with the heating trend, and shows that the percentage contribution of the NE inner flank to the total heat flux decreased between the two acquisitions (Fig. 11b), consistent with spreading of heated zones and soil degassing chimneys onto the outer flanks.

Gas and heat flux partitioning
There is no correlation between T 15 and soil CO 2 flux if we consider the entire data set.One possible explanation for lack of correlation is that it is hard to seal a high-volume accumulation chamber on a rough and hard ground surface, so CO 2 can escape around the edges of the chamber.However, the primary cause is the casehardened crust (cf.Malin et al. 1983), a 10 cm thick cemented surface layer of ash and blocks produced by alteration of hydrothermal salts.Very low permeability of this crust hinders CO 2 escape, and causes gas to flow laterally to cracks in the surface.At the same time the hardened crust will have a higher thermal conductivity than a granular material and is more efficient at transferring heat to the surface (see Tables in Clark 1966).Hot spots have a better correlation with peaks in soil CO 2 flux because they are at permeable, non-crusted locations, although we find dislocation between heat and CO 2 emission points at a meter scale (Fig. 5).This is especially true for measurements made within one meter of each other, as for example at Grid 1, where data are scattered and uncorrelated (Fig. 12b).Instead, measurements separated by several meters are beyond the spatial scale of the dislocation, improving correlation (Fig. 12a).The correlation is also better for Profile 1 in September (R2 = 0.80) than in June (R2 = 0.73) (Fig. 12a) because the number of hot spots increased (Fig. 5).

Soil temperature pulses
Three hypotheses have been used to explain the link between increased seismicity and fumarole temperatures: (1) an external cause, such as rainfall (Nakano and Kumagai 2005), (2) a transient permeability increase leading more efficient heat flux (Harris et al. 2012), and/or (3) increased release of hot gases from the deep magma body (Alparone et al. 2010;Cannata et al. 2012).We see all three effects in our data.
The temperature pulses of late-September and late-October 2021 are synchronous with the first two Very Long Period (VLP) swarms, and the VLP swarm of mid-December is followed by a temperature peak at GRAY_ station (Fig. 8).VLP events registered at Vulcano's hydrothermal system have been interpretated as linked to increases in the gas volume fraction in the fluid, sudden release of pressured gas, and generation of small cracks (Alparone et al. 2010;Cannata et al. 2012).In such case, fluid pressure increases to exceed the mechanical resistance of any obstruction inside hydrothermal fluid-filled cracks (= increased fracturing).The fluid flow then resonates inside the newly opened cracks (= increase of fluid flow) (Alparone et al. 2010).
Rainfall events have also been shown to trigger both VLPs and temperature changes (Nakano and Kumagai 2005), but rainfall events in our data are not systematically related to VLP swarms (Fig. 8).Intense rainfall in July was not associated with any VLP events, but caused T 15 to decrease before recovering to the level recorded before the rain event.Instead, the first VLP events of September occurred a few days after the first rain since July.This could be interpreted as a water influx which destabilized a fully-charged (inflated) and highly unstable system, which was not the case in July when the system was still charging and inflating.Discharge (out flow) then followed in October and December 2021 with the release of predominantly magma-derived fluids (Aiuppa et al. 2022) and elevated levels of seismicity and temperature (Fig. 8).
The depth of the VLP events, and thus the hydrothermal system (cf.Alparone et al. 2010;Cannata et al. 2012), at Vulcano has been estimated to be between 0.5-1.1 km b.s.l. by Alparone et al. (2010) and 0.6-1.2km b.s.l. by Gangemi et al. (2022) in September-October 2021.The synchronicity between the VLP swarms and peaks in T 15 at the periphery of the fumarole field (Fig. 8) suggests that fluids released from depth contributed to near-immediate surface heating.In addition, the high permeability extended over a greater area so that the heat flux became widely distributed, rather than being focused at primary fractures feeding high temperature fumaroles.
In December, fumarole temperatures plateaued at a maximum of 360 °C (Federico et al. 2023, Additional file 1: Appendix I), implying channelization towards the primary fractures.At the same time, the surface temperature response to the VLP swarm was only recorded at GRAY_station, and lagged by roughly two weeks.We interpret this as being due to less efficient heat transfer to the surface due to decreased permeability conditions.This would result from self-sealing and clogging pore space, which can occur in a relatively short timescale at a hydrothermal system (Stix and de Moor 2018;Montanaro et al. 2022;Ardid et al. 2022).By December, sealing had thus already begun at the system periphery, with sealing and channelization thus developing over the timescale of a 4 ± 2 weeks (Fig. 8).

Time-line of unrest
A summary time-line of the evolution of all data sets as the system went into unrest is given in Fig. 13.Onset of unrest is clear in all data sets, but the point in time at which unrest becomes evident is variable.In hindsight, precursors to the unrest were found at HTF, where temperatures began to increase as early as February 2021 (T3, Additional file 1: Appendix I).Increases were detected at GRAY_station during April, and in the Satellite Thermal Index in late-May, and at the fumarole on the crater rim (T2, Additional file 1: Appendix I) in June 2021.During the summer, increases in temperature were also recorded in the diffuse degassing zones (DDZ), by the FTI, at the T 15 line and Grid 1.However, the sampling rate at these locations (every 3-4 months) limits temporal precision as to the exact month of the onset.What we do see, though, is a building thermal onset to the unrest as the system charged.
An increasing trend in soil CO 2 flux began in mid-2020.However, baseline levels were not exceeded until mid-June 2021 (Inguaggiato et al. 2022b).SO 2 flux began increasing in April 2021, exceeding the baseline sporadically in June, and persistently from August 2021 (Inguaggiato et al. 2022b).VLP events increased from zero to several hundred events per day during September (Gangemi et al. 2022), at which time CO 2 metrics also picked up (Inguaggiato et al. 2022b).During the by Inguaggiato et al. (2022b) and the monthly sum of VLP events from Federico et al. (2023).Values used are given in Table 2.The denominator for each grouping are the baseline from years of relative quiescence so as to ensure that each grouping is both weighted and unitless.Here, y RED is the average surface temperature at RED_station during 2020 (4.5 °C), i.e., year of baseline data, and y CO 2 and y SO 2 are the average soil CO 2 and SO 2 fluxes over 2016-2021 as taken from Inguaggiato et al. (2022b).These were respectively 1753 g m −2 d −1 and 29 t d −1 .The sum of the three normalized groupings is then divided by three so that the final value for the VFUI is on a scale between zero and one.et al. ( 2015) developed a Volcanic Unrest Index (VUI) to communicate level of unrest to civil protection and the public, by rating unrest intensity.The approach of Potter et al. (2015) is a semi-quantitative index which uses a greater range of data types than those available here to set an unrest level between 0 and 4. Our index is likewise a semi-quantitative measure of unrest, but which uses those metrics that are available to the case in hand to detect onset of unrest, rather than assign a intensity level to ongoing unrest.In this regard the two approaches are complimentary.Although initialized with three data types (thermal, geochemical and seismic) for Vulcano, our index is designed for potential use at any hydrothermal system entering unrest and where similar monitoring data sets are available.
Our VFUI is plotted in Fig. 14 to track and categorize unrest phases based using monitoring data.The VFUI reveals four phases (Fig. 14): (1) background levels with VFUI < 5, before May 2021, (2) a precursory period through mid-August, when VFUI systematically increased to 20, (3) the onset between mid-August and mid-October, when the VFUI increased to a peak of 80, before (4) declining but remaining at elevated levels (> 50).

Thermal crisis as a sign of unrest: a model
Based on the sequence of thermal events, we can develop a four-step model for the precursory to, and onset of, unrest.While Fig. 15a describes the state of the system before unrest, Fig. 15b-d link the evolution of the thermal crisis to the phases of unrest.
1.The initial state of the system, as given in Fig. 15a, is based on the model of Nuccio et al. (1999), with a scale added for components of the system following Harris and Maciejewski (2000), Gambino and Guglielmino (2008) and Aiuppa et al. (2022).It depicts the situation described in our geological setting section whereby magmatic fluids ascend from a body at around 4 km depth into a mixing zone at 0.5-1 km b.s.l.Above the mixing zone, heat is channelized to the surface through a "chimney" (Harris and Maciejewski 2000), which is a highly permeable zone of fluid, heat and gas ascent surrounded by low permeability clays and self-sealing.2. The thermal crisis was triggered by increased supply of magmatic fluids to the shallow hydrothermal system (Fig. 15b).This occurred between January and June 2021, as indicated by the change in geochemistry of fumarole gas samples (Aiuppa et al. 2022;Federico et al. 2023).3. The increased supply of magmatic fluids caused an increase in heat flux, which was recorded at the sur-  , 15b). 4. Initially, increased fluxes were focused in the heat chimney defined by self-sealing around the fractured conduit of heat and fluid ascent that had developed below the Fossa fumarole field since the end of the 1988-96 unrest (Fig. 15a), but the zone of heat emission began to spread during Spring 2021 indicating expansion of the chimney to peripheral zones of the fumarole field (Fig. 15b). 5.The impermeable seal surrounding the heat chimney failed early in September and was associated with seismic swarms, causing the area of thermal and gas emission to expand suddenly and extensively.At the same time, fluid ascent became much more efficient due to a widespread reduction in permeability so that heat and gas fluxes increased abruptly (Fig. 15c).
This 2021-2022 unrest contrasts with that of 1988-96, when expansion of the heated area was less, but during which there was activation of areas within the existing fumarole field (Bukumirovic et al. 1996(Bukumirovic et al. , 1997)).In addition, the 1988-96 unrest was characterized by a marked increase in maximum temperature (Capasso et al. 1994;Chiodini and Cioni 1995), whereas the 2021-2022 unrest was not.This suggests that the self-sealed heat chimney remained intact during the 1988-96 unrest so as to focus an increased flux in a smaller area; to be reflected in an increase in the maxima 690 °C.Instead, failure of self-sealed zones in September 2021 caused expansion of the emission zone, resulting in a widespread, but less intense, heating of the ground.Indeed, although CO 2 fluxes increased well above those of the 1988-96 unrest, maximum temperatures remained around 380 °C.

Conclusion
After two decades of relative quiescence, Vulcano entered a phase of unrest in late summer 2021 following a four to five month-long period of precursory signals.Precursory signals are mostly marked by increased heat and fluid flow to cause temperatures, gas fluxes and compositions to change.Unrest onset was then marked by an acceleration in the increase of thermal and gas emission, as well as the onset of seismicity and deformation.We show here how ground-and satellite-based thermal parameters can be used to track and characterize the phases of unrest from baseline, through precursory signs to onset.Our understanding of the 2021-2022 unrest was only possible with hindsight.However, results from this multiscale thermal approach, involving a multiparameter data comparison, will be useful for interpreting future events at Vulcano and evolving thermal conditions at less well monitored hydrothermal systems around the world.
The main difference between the current unrest and that of 1988-96 is that a much higher heat and gas flux was involved, spread over a much greater area.This was due to failure of self-sealed zones in September 2021, which did not appear to happen in the 1990s when an increased heat flux instead remained channelized in the self-sealed chimney.Aiuppa et al. (2022) explained the SO 2 degassing during October-December 2021 by the presence, at a depth of 4-5 km, of a relatively small volume (~ 3 × 10 6 m 3 ) of mafic magma.We see the charging of this injection as subtle increases in heat and gas flux, and the seal failure by an abrupt increase in these fluxes.Given the duration of previous unrest at Vulcano, this may continue for several more years.Following Selva et al. (2020), the question is: will the unrest end passively or explosively, when and how? design of the project, acquisition and interpretation of the data, participated in the the original draft of the work and substantial revisions.ISD contributed to the acquisition and interpretation of the data and made substantial revisions.GG contributed to the acquisition, analysis and interpretation of the data, participated in the writing of the draft and made substantial revisions.AC contributed to the acquisition, analysis and interpretation of the data and participated in the writing of the draft.GBillota contributed to acquisition and analysis of the data with software processing.GBoudoire contributed to the acquisition, analysis and interpretation of the data and made substantial revisions.AG contributed to the acquisition and analysis of the data.FG contributed to the acquisition, analysis of the data.MR contributed to the conception of the project and made revisions.

Fig. 1 A
Fig. 1 A Location of Vulcano island in the Eolian Archipelago, Sicily.B Areas of interest in Vulcano island.C Location of permanent stations and punctual measurements.D Location of the limit between hot (GRAY) and cold (RED) zones, and the five zones of intense fumaroles release.HTF High Temperature Fumarole being F5AT, F5 and FA

Fig. 2 a
Fig. 2 a Examples of surface temperature frequency distribution on ASTER images.The distribution starts weakly bimodal with a positive tail and the bimodality increases along with the thermal anomaly.Mode 1 is interpreted as the background reference temperature and mode 2 as the anomalous temperature in the ROI area.b Examples of anomalous pixel locations selected following the temperature frequency distribution

Fig. 3
Fig. 3 Temperature at 15 cm depth on Profile1 west of the fumarole field (2019-2022) (see Fig. 1C for location).The large black arrow marks the hot zone offset into the cold zone after September 2021

Fig. 7
Fig. 7 Surface temperature minus air temperature at GRAY station for 2020 and 2021

Fig. 9
Fig. 9 Expansion of the thermal anomaly associated with the fumarole field between (a) June and (b) September 2021.The rim zone is not visible from this viewpoint.c Comparison of surface temperature frequency distributions from thermal camera images obtained in June and September 2021.d Surface temperature frequency distribution in September 2021

Fig. 10
Fig. 10 ASTER-derived STI.Three scenes (09/29/21, 01/26/22 and 03/14/22) were still selected even though the cloud cover was partially covering the crater (colored star mark).The July 2021 storm events are marked by a blue triangle.The STI and NBp baselines are respectively placed at 40 and 10 pixels Tmax crater * Tmean crater * SD crater * Nb p Tmax control * Tmean control * SD control , (temperature) of the an anomaly increasing (Additional file 1: Appendix D, E).

Fig. 11 a
Fig. 11 a Comparison of Radiative Heat Flux (RHF) from ASTER, VIIRS and ground-based thermal imagery.Minimum bound for ASTER heat flux is using the first mode of temperature distribution for pixels in the Fossa ROI as background temperature (BG1).Maximum bound uses the minimum temperature from the same ROI (BG2).b Total heat flux time series as calculated from ASTER on the whole exposed anomaly and compared to the heat flux from ground-based thermal imagery (targeting only the inner flank of the northern slope)

Fig. 12 a
Fig. 12 a Comparison of T 15 and soil CO 2 flux along Profile 1 and (b) across Grid 1

Table 1
Temperature statistics for the camera thermal images

Table 2
(Inguaggiato et al. 2022b)te VFUIParameters from left to right: monthly average temperature of fumaroles T2 and T3 (Appendix I), monthly average soil CO 2 flux and SO 2 flux(Inguaggiato et al. 2022b), monthly sum of VLP events (Instituto Nazionale Di Geofisica e Vulcanologia, 2021b, 2021a), STI and monthly average of surface temperature minus air temperature at GRAY station.The parameters generated in the course of this study are bold.STI's value in italics are estimates to fill in the gaps in ASTER images