Data Inventory and Location of seismic Signals Recorded During the 2021 Unrest on the Island of Vulcano, Italy

1. Introduction

In September 2021 geochemical and geophysical parameters routinely monitored on the island of Vulcano, Italy, showed relevant anomalies, causing major concern among the population and Civil Protection authorities.

Vulcano is part of the Aeolian Archipelago, an arc-shaped volcanic chain of seven islands located in the Southern part of the Tyrrhenian Sea, Italy (Figure 1). The tectonic framework of this region is primarily controlled by the complex N–S to NNW–SSE trending convergence between the African and Eurasian plates ([1] and references therein). At Vulcano Island, the tectonic setting is characterized by the presence of en-echelon fault segments trending NNW–SSE, with oblique (right-lateral) transtensional kinematics referred to the so-called Aeolian–Tindari–Letojanni Fault System (ATLFS in Figure 1). In the northern part of the island, this fault system interacts with the WNW–ESE transpressional fault belt known as the Sisifo–Alicudi Fault System (SAFS in Figure 1; [1,2]).

Fumarolic activity affects in- and near-shore of the island since its last eruption in 1888-1890. While no further eruptions have occurred, the area has experienced several periods of volcanic unrest, which are typically defined by rising fumarole temperatures, shifts in gas chemistry, heightened seismic activity, and measurable areal dilatation ([3] and references therein).

Geophysical studies document processes related to either regional tectonics [4,5,6,7] or changes in the shallow hydrothermal system [8], also favoring slope failures [9].

As with most active volcanoes worldwide, Vulcano is characterized by several types of seismic signals. They are related to i) the dynamics of the shallow hydrothermal system, and ii) fracturing processes caused by shear failure and slip along a fault plane with the so-called volcano tectonic (VT) earthquakes e.g. [10], generally associated with the dynamics of the regional fault systems [5,6,7,11]. Regarding the seismicity related to hydrothermal activity inland, it has been studied since the 1970s [12], and up to the 2021 unrest seismic events were classified into different categories based on their waveforms and frequency content, ranging from Long Period (LP; frequency content 0.5 - 5 Hz) to High Frequency (HF; frequency content 5 - 25 Hz) events, including monochromatic events and tornillos [13]. Their foci were generally located within ~2 km underneath La Fossa Crater area [14,15].

The 2021 unrest was marked by a drastic increase in hydrothermal seismicity (Figure 2), dominated by an unprecedented occurrence rate of Very Long Period (VLP) events never before recorded in such high numbers at Vulcano. For these events, Federico et al. [3] reported spectral peaks of ~0.3 Hz and a centroid location North of La Fossa Crater at depth of ~1 km. The source location of these events suggested a connection between pressurized fluids and geological structures [16,17]. In literature, VLP events in volcanoes capable of producing phreatic eruptions appear indeed to be reliable indicators of pressurization of the shallow hydrothermal system, e.g. [18].

In this study we present an application of machine learning to seismic transients recorded on Vulcano between September 2021 and December 2022 to create an inventory that includes various signals, from high-frequency (including anthropic noise) to VLP events. This application was also tested when a resumption of microseismicity occurred between May and June 2024. In addition, we located selected events in order to explore: i) whether or not fluids acted in various sectors of the island, tracing the diverse ‘signature’ of the clusters from pattern classification; ii) the focal depths involved; and iii) how the location of the seismic source/s can be related to other changes in the volcanic system.

2. Seismic Network

The Istituto Nazionale di Geofisica e Vulcanologia, Osservatorio Etneo (INGV-OE) runs the 24/7 seismic monitoring of Vulcano and the other Aeolian Islands with digital stations equipped with three-component broadband seismometers (Figure 1). In particular, IVCR is the reference station for seismic surveillance of Vulcano due to its proximity to La Fossa Crater.

The signals are sampled at 100 Hz and telemetered at the INGV-OE data acquisition center in Catania, where they are acquired in real time and stored in a digital repository. Soon after the 2021 unrest began, the permanent network was boosted with the installation of a mobile seismic network consisting of five stations. Three stations were installed on the island of Vulcano and two in the southern part of Lipari (Figure 1). The mobile stations were also digital and equipped with three-component broadband velocimeters.

3. Tectonic Earthquakes

A moderate increase in volcano tectonic seismicity occurred in late October 2021 [7], well after the onset of the 2021 unrest. Cocina et al. [19] documented three time intervals with a relatively high occurrence rate of VT earthquakes: October 30 to December 31, 2021 (42 earthquakes; 0 ≤ ML ≥2.6); March 31 to April 30, 2022 (66 earthquakes; 0 ≤ ML ≥2.6); and December 4 to 31, 2022 (34 earthquakes; 0 ≤ ML ≥4.6). In the third time interval, an earthquake reached the highest magnitude value during the entire unrest (ML4.6). The hypocenters were located in the northern part of the island (with depth < 5km b.s.l.), in the eastern part (with depths mostly down to 10 km b.s.l.), and offshore, in the southwestern part (with depths up to 15 km and more). These findings confirmed the existence of seismic volumes, which were also active during the previous 34 years [20].

4. Application of Unsupervised Learning

An abrupt change in seismic activity started on September 7, 2021, with the occurrence of a variety of signals, an example of which is given in Figure 3. The seismic radiation was rather complex, characterized by the intermingling of waveforms and frequency content, likely related to various sources.

An important aspect, which hindered the definition of a reliable monitoring concept, was the presence of noise, mostly originating from human activities, such as cars, hydrofoils, or helicopters. Standard monitoring concepts, such as the well-known RSAM [21], or spectrogram-based methods, do an excellent job on volcanoes like Etna, where the persistent background radiation (so-called volcanic tremor) provides important information about the state of the volcano, e.g. [22]. Conversely, the application of a priori fixed criteria would have been of questionable value on Vulcano due to the variety of exogenous and endogenous origin of the signals.

Given the multitude of recorded seismic signals and the limited a-priori consolidated knowledge about their typology and origin, the application of a machine learning technique with unsupervised learning was a convenient choice to screen the data. In unsupervised learning we seek indeed to discover structures in a dataset, such as groups of similar patterns, by separating the elements that belong to these groups without the supervision of an expert. Our analysis was based on the creation of the so-called Self Organizing Maps or SOM [23]. SOM have been successfully applied to seismic signals on various volcanoes, e.g. [24,25,26,27,28]. An in-depth description of SOM is given in the section 9 ‘Methods’ (subsection 9.1 Self-Organizing Maps).

We considered ~6600 seismic transients recorded at the IVCR station (Figure 1) from September 1, 2021 until December 31, 2022. The transients were extracted from the continuous data stream by applying a simple short-term average/long-term average (STA/LTA) peak detector and a minimum amplitude threshold of 12500 counts, corresponding to 20.12 µm/s. These trigger conditions allowed us sampling not only events, but also anthropic noise. Each transient was time windowed with a length of 20 s (2000 points). Instead of using waveforms directly, we derived feature vectors by calculating amplitude spectra and, alternatively, the auto-correlation function (ACF) of the signals. In this choice, we followed Falsaperla et al. [29], who successfully used ACF in the context of supervised classification of explosion quakes recorded at Stromboli. In our case study, ACF was obtained by calculating power spectra using an FFT length of 8192 points and taking its inverse. We considered 819 points of the ACF to create the feature vectors, and reduced the number of components for the spectra averaging the power spectral values of 10 single density values. In Table 1 we summarize the SOM quality parameters obtained by applying the KKAnalysis software [30,31] (see section 9. ‘Methods’) to both ACF and spectra.

Ideally, the dispersion of feature vectors in the nodes is equal to the total dispersion measured over all feature vectors. In reality, the dispersion in the nodes is slightly less, i.e., the nodes do not fully represent the total dispersion of all feature vectors. The smaller the Quantization Error (QE), the better the representation of patterns given by the nodes. The Topological Error (TE) describes how well our map represents the neighboring relationship among the patterns. It is obtained by counting how often patterns close to each other on the map are not close in the original feature space. We tested different normalization options, which can actually affect the quality of our results. From the results reported in Table 1, we may prefer the ACF applying a “range normalization” (the values are normalized so that they fall within a range between 0 and 1). This option performs at best with TE, and very good with QE. A valid alternative is the option autocorrelation with logarithmic or logistic normalization. In particular, logistic normalization mitigates problems arising from outliers, which can affect range-based normalization. In the following, we will focus on this last option.

Figure 4a shows the SOM obtained by projecting the prototypes of the 6600 seismic signals onto a 2D map using the ‘small map’ option of the KKAnalysis software. Each prototype identifies a family of signals and is represented as a colored hexagon called Best Matching Unit (BMU).

Neighboring families are assigned similar colors. The numbers in the hexagons indicate the number of patterns belonging to each BMU. We obtained 102 BMUs from the original 6600 input data. In the 2D map of Figure 4a each BMU is numbered by counting the hexagons from the bottom row, with a zigzag motion from left to right towards the top. Also, the frequency content of the signals increases on the map from left to right. The waveforms shown in Figure 4b for five BMUs (briefly described below) provide examples of the variety of signals recorded. BMU19 (the hexagon filled with a bluish color in Figure 4a) forms a family of 72 events; the waveforms in Figure 4b document that they are multiplets characterized by prevailing very low frequencies.

BMU10 is a hexagon marked with light-green colors at the center of the SOM in Figure 4a. Looking at the waveforms, we again notice that very low frequencies prevail. They differ from those of BMU19 due to a high-frequency signal that starts shortly after the peak amplitude (Figure 4b). The waveforms of the BMU94 still have a well-developed low-frequency content, but the high-frequency part is significantly increased (Figure 4b). BMU17 and BMU99 are marked by hexagons filled in red colors in Figure 4a, and encompass waveforms with high-frequency content (Figure 4b).

We tested the spatial separation of the foci for the events belonging to different BMUs. We calculated the differences between the centroids of the foci and applied the Hotelling’s T² statistics for the so-called distinction of ensembles, e.g. [31,32]. The T² test follows the F distribution, and allows us to assess the significance of the separation of two or more ensembles, the patterns of the BMUs in our case study. The results obtained are reported in Table S1 and Figure S1 in the Supplementary Data.

Figure 5 depicts the distribution of all the 102 BMUs over time, showing the presence of a prototype with its colored symbol from September 1, 2021 to December 31, 2022.

Some clues on the origin of the signals can be obtained by looking at the occurrence times. For example, while the signals of BMU19 occur any time, those forming BMU99 are absent in the late evening and early morning, and are related to anthropic noise partly generated by boat traffic, such as hydrofoils (for a comparison, see Figure S2 and Tables S2- S3 in the Supplementary Data). Yellow to light-green symbols, representing events such as the multiplets of BMU19 (Figure 4a, b), are visible in Figure 5 only between September and October 2021. The ‘burst’ of seismic activity to which they belonged was only temporary, and they disappeared in the following months. In December 2021, the colors switch to blue and blue-green, and are associated with a ‘burst’ of seismic activity with the simultaneous presence of numerous BMUs (Figure 5). In the next months, blue symbols prevail and, from May 2022 on, a color change with blue-green and dark green tones is observed again. The grey bands in Figure 5 highlight the three most relevant time spans in terms of occurrence rate and magnitude of the VT earthquakes according to the results from [19].

5. Hypocentral Location

Moving from the evidence of many families of seismic signals, whose occurrence rate changed over time, we searched for their locus and depth. Firstly, we considered 500 seismic signals with a good signal-to-noise ratio, belonging to different BMUs, to manually pick phase arrivals and subsequently calculate their hypocentral location over the whole time span analyzed. From the data sample of 500 signals, we obtained 2950 P pickings, which allowed us to locate 257 events between September and December 2021, and 142 events between January and December 2022, for a total of 399 events. We discarded 101 data because: i) some events had not reliable location (for example, they had too low energy to appear above the background noise or there was the temporary missing acquisition of some station/s); and ii) a manual screening of the waveforms highlighted that the transients were caused by noise. We applied the following criteria to consider a reliable location: number of stations ≥ 6, azimuthal gap ≤ 200°, root mean square error (RMS) ≤ 0.2 s, and horizontal and vertical error ≤ 0.5 km. The stations considered for the picking collection were those marked in Figure 1, and encompassed both permanent and mobile sensors set up on the islands of Lipari and Vulcano.

Figure 6 depicts the daily occurrence rate and the cumulative number of events with frequency peak < 1Hz. In addition, in the figure we show the epicentral location and foci distribution in N-S and W-E cross-sections for some BMUs with the highest daily occurrence rate of events in the following time spans: September 6 – October 29, 2021, December 1 – 20, 2021, May 22 – June 12, 2022, July 9 – September 3, 2022, and May 11 - June 24, 2024.

The inspection of Figure 6 reveals that there were changes in the epicentral location and focal depth of the seismic events, which became more clustered and deeper over time.

The plot of all localized events from September 2021 to December 2022 is given in Figure 7. Here, the hypocentral locations refer to 56 BMUs. The epicenters are clustered in La Fossa Crater, with only a few scattered around it. The maximum depth attained in 13 months was 2.2 km (see Figure S3 in the Supplementary Data). Conversely, we clearly observe two main focal volumes separated by a discontinuity at about 600 m b.s.l. We will return to this finding in the ‘Discussion’.

Table 2 summarizes some statistics of the obtained results, concerning mean and maximum depth, mean gap and RMS, mean errors in horizontal (seh) and vertical (sez) direction.

Secondly, we analyzed the two ‘bursts’ of activity highlighted in Figure 5, which could shed light on the questions: ‘Where was the source of seismic activity between September and October 2021, namely at the beginning of the unrest?’; and ‘Where were the seismic sources in December 2021, when numerous BMUs occurred simultaneously?’.

For the first burst, we focused in particular on the days between September 30 and October 1, 2021, for which we obtained 11 BMUs (blue dashed hexagons in Figure S4 in Supplementary Data), the highest number of concurrent BMUs in this time span. The hypocentral locations refer to 62 events, which occurred mainly North and West of La Fossa Crater at very shallow depths (Figure 8). The statistical analysis in Table 3 highlights a mean depth of 0.2 km b.s.l., the lowest of the entire dataset (Table 2).

The second burst of seismic activity of 2021 was in December. We obtained the hypocentral location of 123 events, which occurred in the days between 7 and 9. It is remarkable that these events belonged to 37 BMUs (red dashed hexagons in Figure S4 in Supplementary Data), the maximum of concurrent BMUs during the entire unrest. Also remarkable is that their epicenters clustered inside La Fossa Crater and clearly involved two focal volumes: from the surface to ~500 m b.s.l. and from 600 m to 1 km b.s.l. (Figure 9).

The maximum daily occurrence rate reached during the two afore-mentioned bursts was 67 and 58 events, respectively. From January to December 2022 two more bursts of activity (from May 22 to June 12 and from July 9 to September 3) exceeded the daily occurrence rate of 30 events, with a maximum of 42 and 32 events, respectively. The hypocentral maps in Figure 6 document that also these events remained clustered inside La Fossa Crater, but with relatively deeper foci up to 2.2 km b.s.l. (see also Table 2). By the end of November 2022, the occurrence rate decreased sharply, remaining well below 10 events per day over the next 14 months. This low rate of seismic activity marked a trend towards background conditions, also highlighted by the diagram in Figure S5 (in Supplementary Data), which depicts the cumulative energy of the seismic signal calculated from the maximum squared amplitude of the transients. However, a sudden, temporary increase (with a maximum occurrence rate of 23 events) occurred again between May 11 and June 24, 2024 (Figure 6). This new burst of seismic activity led us to examine its locus and verify whether the new transients had seismic signature comparable to the previously identified patterns, namely whether they belonged to the BMUs of Figure 4a. Overall, 148 events out of 180 transients had a reliable location according to our criteria. They shared similar hypocentral parameters (and corresponding errors) of those recorded during the unrest between 2021 and 2022 (Figure 6 and Table 2). However, their cumulative energy remained well below the peaks reached in 2021 and 2022 (Figure S5 in Supplementary Data).

To verify whether the new signals belonged to the BMUs of Figure 4a, we analyzed these data using the signal inventory created with the transients recorded from September 2021 to December 2022 as a ‘reference’ dataset. The strategy we followed for this test is explained in the section 9. ‘Methods’ (subsection 9.2 ‘Test of online application’). We observed that the new transients shared the same BMUs as the signals of our previous inventory, confirming the hypocentral parameters results.

6. 3D Analysis of the Hypocenter Distribution

An inspection of Figure 10, in which the 547 localized events are color-coded by year, reveals temporal variations in the hypocentral distribution across the N-S and W-E cross-sections. Overall, the seismic volumes involved from 2021 to 2024 were common, despite a greater scatter affected the 2021 events, which had a relatively wider aerial distribution (red dots in Figure 10). Shallow foci prevailed in the whole dataset, with a mean depth of 0.5 km b.s.l. (see Table 3). Notably, less than 3% (16 events) of the foci were deeper than 1.5 km b.s.l.

To shed light on finer details, the hypocentral locations for 2021, 2022, and 2024 are displayed in separate 3D plots in Figure 11, providing a clearer view of the seismic volume involved each year as described below.

2021. (September – December)

The 3D diagram for the year 2021 documents the activation of a shallow source volume (Figure 11a). The involved region has an extension of ~3 km on the map, from West of the Baia di Levante to the southern flank of La Fossa Crater (see dashed line in Figure 10). ~76% of the events (195 out of 257) had foci ≤ 0.5 km b.s.l. 59 events (23%) belonged to another main source volume located below La Fossa Crater between 0.6 and 1.4 km b.s.l. (Figure 11a). There were only three events with depth ≥ 1.5 km, the deepest of which (~1.8 km b.s.l.) occurred on December 12.

2022. (January – December)

All events in 2022 had epicenters inside La Fossa Crater. 68 (48%) out of 142 had depths ≤0.5 km b.s.l.; 8 events (5%) had foci between 1.5 and 2.2 km b.s.l., while the remnant 66 events (47%) formed a well-clustered group with intermediate depths (Figure 11b). The deepest hypocenter (2.2 km b.s.l.) was attained on July 11 and belongs to the time span from July 9 to September 3, 2022 (forth box in Figure 6) during which the foci became deeper on average with respect to the previous months of the unrest.

2024. (May – June)

On the whole, the 148 events of this data sample showed similar characteristics to those of the years 2021 and 2022 (Figure 11c). Most of the epicenters were located within La Fossa Crater; however, some events also occurred on the southern and western flanks of the crater. 76 events (51%) had shallow foci (≤0.5 km b.s.l.). A source volume between 0.6 and 1.4 km b.s.l. accounted for 67 events (45%). Only 5 events had depths from 1.8 to 1.9 km b.s.l.

7. 3D Analysis of Adaptive and Gaussian Clustering

To explore the level of clustering of the hypocenters of the entire dataset, we applied adaptive and Gaussian clustering to the latitude, longitude, and depth of the 547 events recorded in 2021, 2022, and 2024. A description of these clustering methods is provided in the section 9. ‘Methods’ (subsection 9.3 ‘Adaptive and Gaussian clustering’). The results obtained were depicted in a 3D plot for each technique (Figure 12).

The results of the adaptive clustering indicate partitions on both horizontal and vertical scales (Figure 12a). The horizontal partitions highlight: i) a group of very superficial events (green) and ii) the critical depth around 600 m, which marks the separation of the three shallowest clusters (green, dark blue, and red) from the deepest ones (light blue, and pink).

The vertical partition refers to three clusters: red, pink, and light blue. The red cluster is the only one with an almost (apparently) continuous distribution from depth (~1 km b.s.l.) to the surface. The pink and light-blue clusters have the deepest foci, from about 600 m b.s.l. down to ~2.2 km b.s.l. Unlike their scattered distribution from 1.3 km b.s.l. to greater depths, they form a dense focal volume at about 1 km.

The results of the Gaussian clustering mainly depict partitions on horizontal scale. Again, we recognize a discontinuity at a depth ~600-700 m, above which foci of the green and pink clusters prevail. The separation between the shallowest (green/pink) and the deepest (light blue) clusters is evident, even though a few events of the shallowest cluster reach the depth of about 1 km b.s.l. Vertical partitions affect only the first kilometer b.s.l.

8. Discussion

Seismic signals alone do not tell much about a volcano unrest. Rather, we understand this information as an indirect hint, such as a high body temperature is associated with some pathology. In the context of the Vulcano unrest started in 2021, Federico et al. [3] discussed multidisciplinary data recorded until April 2022 by dense monitoring networks of sensors run by INGV. Of particular concerns were parameters related to thermal ground waters (e.g., the increase in temperature, conductivity, salinity, CL/SO₄ ratios), which showed significant anomalies, and the high emission of CO₂, which entailed a considerable risk for humans and animals [33]. These signs of an impending crisis, affecting the hydrothermal system, were also supported by the inflation of the cone of La Fossa Crater and by the occurrence of the variety of seismic events with low- and very low-frequency component (frequency peak ~0.3 Hz) from September 2021 on [3].

From the seismic viewpoint, VT seismicity neither heralded nor played a major role at the onset of the unrest. Earthquakes occurred indeed in three time intervals starting from October 30, 2021 (Figure 5), and the maximum magnitude (ML4.6) was reached in December 2022, more than one year after the unrest began. Conversely, there was the need to classify and locate the varying typologies of seismic transients (many of which with frequency content below 1 Hz) to shed light on the location of their source volume/s and improve online monitoring using a data inventory as a reference. Note that seismic events with origins related to fluid dynamic were numerous during past unrests, such as in 1985, e.g. [34]. However, the long- and very-long period component could not be recorded due to the use of technical equipment with short-period seismometers until the 1990s. Only recently Lupi et al. [35] found in hindsight the temporary occurrence of few VLP events in June 2018 preceding an increase in volcanic gas emissions recorded in September 2018 [36].

In the following, we discuss the results obtained from our analyses to answer the three questions: i) whether or not fluids acted in various sectors of the island, tracing the diverse ‘signature’ of the clusters from pattern classification; ii) the focal depths involved; and iii) how the location of the seismic source/s can be related to other changes in the volcanic system.

8.1. Did Fluids Act in Various Sectors of the Island?

Variety and number of the seismic signals recorded at the beginning of the seismic unrest in September 2021 moved us to create a data inventory to distinguish the diverse typologies of transients. For this purpose, we opted for the SOM, a machine learning technique for which no a priori knowledge is required (the so-called ‘unsupervised learning’). We obtained 102 BMUs, each one encompassing a typology of transients from very long-period (VLP) events to high frequency signals (Figure 4a). Useful was the identification of transients that do not represent a matter of alert, such as those related to anthropic noise (e.g., boats, hydrofoils, and cars; for example, BMU99 in Figure 4b and pink symbols in Figure 5), with a typical daytime occurrence (Figure S2 in Supplementary Data). Conversely, the occurrence of VLP events, such as those of BMU19 (Figure 4b), along with other transients with low- and very low-frequency content throughout their duration (Figure 4b, Figure 8 and Figure 9), documents the involvement of fluids and can be indicative of potentially hazardous pressurization conditions of the hydrothermal system. The hypocentral analysis of these events highlights that they remained confined mostly in the region of La Fossa Crater over time (Figure 6 and Figure 7, and 9). A slightly greater epicentral scatter occurred during the first burst of seismic activity in September-October 2021, namely, at the start of the seismic unrest (Figure 5 and Figure 6). In that time span, we found shallow events (on average 200 m b.s.l.) affecting also the northern and western part of the island (Figure 8, Figure 10 and Figure 11), where Federico et al. [3] measured a huge increase in gas and vapor emission. The waveforms of these seismic events had a mingling of high- and very low-frequencies in varying parts of their signature (Figs 3, 4b). The light-green and yellow BMUs, which identify these events on the SOM, disappeared in the next months (Figure 4 and Figure 5).

8.2. Focal Depths Involved

The analysis of the BMUs revealed the occurrence of several families of multiplets in bursts (Figure 5, Figure 8 and Figure 9). This type of seismic activity, with fluids and fracturing interaction, is typical in volcanic and hydrothermal systems [37]. We calculated the hypocentral location of 399 events, which occurred between September 2021 and December 2022 and belonging to 56 BMUs (Figure 7). The results highlighted that they were mostly below La Fossa Crater within two main shallow focal volumes (from the surface to ~500 m b.s.l., and from 600 to ~1 km b.s.l.). As previously stated, the average depth was shallow (200 m b.s.l.) during the first burst (between September and October 2021), while relatively scattered foci reached the maximum depth of 2.2 km b.s.l. in 2022 (Table 3).

After more than 14 months of low occurrence rate, there was a sudden, temporary increase in seismic activity between May and June 2024. Consequently, we analyzed 180 additional seismic transients recorded in that time span. A test application of the SOM to the new data confirmed the non-destructive source of the new recorded signals, which shared families and, for 148 locatable events, location, and depths similar to our previous inventory (Figure 10 and Figure 11). This evidence pinpoints the importance of the 2021 unrest in terms of duration with respect to previous episodes, such as those in 1985 and 1988 [20].

The 547 events of our dataset were temporally distributed as follows: 257 from September to December, 2021, 142 from January to December, 2022, and 148 from May to June, 2024. Overall, we hypothesize that their seismogenic sources were confined to a fracture network where the build-up of overpressure of fluids and vapor of the hydrothermal system acted. The resulted forced unsealing was promoted by fracture-gas transfer and/or gas-fracture transfer, as revealed by the seismic event signature, with varying intermingling of the very low-frequency content of the signal, in agreement with the conceptual model proposed by Konstantinou [38] for hydrothermal systems.

8.3. How can the Location of the Seismic Source/s be Related to Other Changes in the Volcanic System?

The peculiarities observed in the time and space distribution of the seismic signals mirror also in data regarding other parameters monitored at Vulcano Island. We mentioned the sharp increase of soil CO₂ fluxes from September 23 on, which reached its climax in October 2021 [3]. A similar development occurred in the water temperature of two wells, as well as salinity and conductivity. Besides, enhanced activity in terms of VLP events from the second half of May 2022, was accompanied by the emission of whitish liquids along the beach, observed on May 22 - 23, 2022 [3].

The 2021 Vulcano unrest attracted the interest of the remote sensing community, which acquired a wealth of data, furnishing further evidences of the state of the volcano. For example, Campus et al. [39] applied MIROVA, which is a satellite-based near real-time method for the detection of volcanic hotspots. The MIROVA algorithm, originally developed for the “Moderate Resolution Imaging Spectro-radiometer” (MODIS) data, has been amplified through recent years allowing to handle date from different sensors, such as the “Visible Infrared Imaging Radiometer Suite” (VIIRS) [40,41]. Figure 13a shows the time series of the number of alerted pixels (Npix) within a distance of 1 km radius from La Fossa Crater from May 2021 to December 2022 [39]. The increase of Npix clearly highlights two main phases of the unrest, which are concurrent with the trends of the Volcanic Relative Power (VRP) in Figure 13b and the time development of microseismicity with frequency content < 1 Hz (Figure 13c).

Localized deformation at volcanoes, such as those measured at La Fossa Crater, are indicators of pressurization inside a hydrothermal system, thus potential precursors of phreatic explosions [42,43]. Di Traglia et al. [17] obtained a deformation source located at La Fossa Crater with a prolate spheroid geometry, modelling geodetic (MT-InSAR and GNSS) data. This deformation source and the hypocenters of the seismic events of our data inventory over the same time span (September 2021 – December 2022) share the same location (Figure 7). The location of a resistivity anomaly is also consistent with the zone separating the two main focal volumes visible in Figure 7. This anomaly is on the order of hundreds of Ω m and was detected beneath La Fossa Crater at approximately 600 m depth using 3D magnetotelluric modeling [16]. A discontinuity at that depth was also identified by Fulignati et al. [44] and interpreted as the transition between two zones with different permeability and geothermal gradient. Accordingly, we assume that unsealing processes promoted the response of a fracture network above and/or below the afore-mentioned discontinuity due to the action of the hydrothermal system.

9. Methods

9.1. Self-Organizing Maps (SOM)

SOM are made up of a number of nodes, each of which represents a number of patterns. The nodes are small clusters of their own, and can be understood as prototypes of patterns they represent. Throughout an iterative training process, the node weights are adjusted at every time step so that the sum of the distances between original data and their representing prototype nodes converges to a minimum. Formally, during the training of SOM, one can minimize the sum of the distances with the relationship

$$D_{ij}=\sqrt{{(W_i-V_j)}^T(W_i-V_j)}$$

(1)

where V_j is the normalized input feature vector and W_i the weights stored in the nodes. V and W have the same dimension. A core step is the identification of the closest node to the actual input vector, i.e., the best matching unit (BMU) for the j_th pattern. Neighboring nodes lying within a certain radius of influence are considered as well. Once BMU and the nodes falling within the area of influence are identified, the weights are gradually adjusted according to the so-called learning rate λ, which decreases with time t. A second parameter φ - called influence radius - describes the dependence on distance ∆ of the upgrade of a node. φ (∆,t) is maximum for the BMU, whereas nodes outside the radius of influence are not upgraded at all. Similar to the learning rate λ, also φ (∆,t) decreases with time. Various options are available for the shape of φ, such as the Gaussian one, which was adopted in our application. The upgrade of weights follows the relationship

$$W_i(t+1)=W_i(t)+\varphi (\Delta ,t)\cdot \lambda (t)\cdot D_{ij}(t)$$

(2)

An appealing property of this scheme is that topological relationships among the original data also hold for the BMUs. Consequently, patterns represented by neighboring nodes in the SOM are truly close to each other also in the original data space. This characteristic of topological fidelity is indeed achieved by using the influence radius φ (∆,t). An in-depth description of SOM is reported in the most recent edition of Kohonen’s classical work [45]. SOM can be particularly instructive when the weights take the form of a color code (Figure 4a). For this purpose, Principal Component Analysis is applied, projecting the multi-dimensional weights vector in a 2D representation space made up of the principal axes z1 and z2. These axes are normalized in the range [0….1]. A third value z3 is obtained either as 1–z1 or 1–z2.

9.2. Test of Online Application

Stability of preprocessing and affordable computational effort allow us to envisage an online application for the pattern classification of seismic transients recorded on Vulcano using SOM. The application of such a routine to the data stream continuously acquired requires the use of a representative ‘reference’ dataset, which is processed together with the new patterns. The reference dataset is chosen in the light of the expert knowledge and should be representative for all relevant classes of signals. We tested the possibility of a similar application considering the inventory of the signals created from September 2021 to December 2022 as a reference dataset, and the burst of seismic activity occurred between May and June 2024 as the new data to classify. In our application, these new data were stored in a ring buffer whose size had to be small compared to the reference dataset. This was a necessary condition to avoid that the basic characteristics of the SOM (in particular the number of nodes and the shape of the map) did change when the ring buffer was added to the reference dataset. We adopted a ring buffer with 45 events - less than 1% of the 6600 patterns in the reference dataset. For the sake of convenience, we updated the ring buffer four times, passing successively 45 out of the 180 transients recorded between May and June 2024. In this way, the essential characteristics of the SOM (102 BMUs) remained stable, proving that: i) the new signals shared the same characteristics of the transients forming the previous inventory; and ii) the effectiveness of the classification strategy and its possible online application.

9.3. Adaptive and Gaussian Clustering

Cluster analysis based on unsupervised learning establishes similarities of patterns using an a-priori defined metrics. The choice of this metrics depends on the purposes of the users applying the analysis. A further issue is the choice of the clustering strategy. In Partitioning Clustering, we have to fix the number of clusters a priori. Typically, one starts with a small number of clusters, and augments this number in the following steps, until a solution (a partition in agreement with some formal criteria or that is feasible for the user’s goals) is found.

In adaptive clustering schemes, we consider the (squared) sum of the distances (${\mathit{D}}_{\mathit{t}\mathit{o}\mathit{t}}$) measured within the clusters, i.e.

$${\mathit{D}}_{\mathit{t}\mathit{o}\mathit{t}}=\sum _{j=1}^k{\sum }_{i=1}^{m_j}{(\mathit{x}}_{\mathit{i}}-{\overline{\mathit{x}}}_j{)}^T{{\mathit{G}}_{\mathit{j}}}^{-1}\left({\mathit{x}}_{\mathit{i}}-{\overline{\mathit{x}}}_j\right)$$

(3)

which we want to be minimum. Here ${\mathit{x}}_{\mathit{i}}$ is the feature vector of the i_th pattern (the coordinates of the hypocenter), ${\overline{\mathit{x}}}_j$ is the centroid of the j_th cluster, and ${\mathit{G}}_{\mathit{j}}$ the dispersion matrix of this cluster.

The metric resembles the Mahalanobis distance with the difference that we use the elements of the dispersion matrix instead of the covariance matrix. Cluster hulls form elongated ellipsoids, accounting for possible correlations between the features. Note that ${\mathit{D}}_{\mathit{t}\mathit{o}\mathit{t}}$ will always decrease augmenting the number of clusters. However, the rate of ${\mathit{D}}_{\mathit{t}\mathit{o}\mathit{t}}$ decrease flattens when a suitable choice for the number of clusters is achieved. This was the case when we chose five clusters (Figure 12a).

Besides, we can apply clustering focusing on the probabilities that a pattern belongs to a certain cluster. The Expectation Maximization (EM) algorithm is based on the a priori assumption that the distribution of our samples is obtained from the sum of multivariate normal distributions. The resulting partition is referred to as “Gaussian Mixture Model”. Each cluster is described by the parameters ${\overline{\mathit{x}}}_j$ and ${\mathit{C}}_j$(i.e., their centroid vectors and covariance matrices). As for the adaptive clustering, we adopted 5 clusters (Figure 5b). Given cluster j, the probability that a pattern x_i belongs to this cluster follows from

$$p\left({\mathit{x}}_i\right|j)={\mathit{C}}_j^{-1/2}\mathrm{e}\mathrm{x}\mathrm{p}(-({\mathit{x}}_i{-\overline{\mathit{x}}}_j\left){\mathit{C}}^{-1}\right({\mathit{x}}_i{-\overline{\mathit{x}}}_j{)}^T)$$

(4)

At the same time, we should be aware that clusters j and k themselves are not equally likely to occur. Therefore

$$p\left(\mathit{x}\right)={\sum }_{j=1}^{M}p\left(\mathit{x}\right|j)P_j$$

(5)

where P_j is the probability a priori to see cluster j. In other words, all M distributions contribute in a weighted way to $p\left(\mathit{x}\right)$. Typically, we do not know a priori the parameters of the distribution, and our task is now the blind identification of all the parameters, i.e., the centroid vectors and covariance matrices of the single cluster, ${\overline{\mathit{x}}}_j$and ${\mathit{C}}_j$, as well as the a priori probabilities $P_j$. For this purpose, we can apply the “Generalized Mixture Decomposition Algorithmic Scheme” (GMDAS, for more details see [46]).

10. Conclusions

The application of machine learning techniques to the seismic transients recorded during the unrest that began at Vulcano in September 2021, allowed us to identify several families of signals. Specifically, we chose SOM to create the data inventory due to: i) the complexity of seismic records, which hindered standard monitoring (e.g., RSAM application); and ii) the difficulty in defining a-priori classifications for techniques based on supervised learning.

Moving on from the SOM results, we found that the families of seismic events in the time span September - October 2021 were very shallow (on average ~0.2 km b.s.l.) and scattered around La Fossa Crater (Figure 8, Figure 10 and Figure 11). Many of them appeared as events with a mingling of frequencies throughout their duration (Figure 3). From mid-October 2021 on, the families of the seismic events clustered below La Fossa Crater (Figure 7 and Figure 10). In the three years analyzed (2021, 2022, and 2024), the 547 hypocenters localized were well confined (mean vertical error 0.4 km), mostly within two volumes up to a depth of ~1 km b.s.l. Only less than 3% (16 out of 547 events) of the foci were deeper than 1.5 km b.s.l.

We surmise the events of September – early October 2021 opened the paths to gas and vapor leading to their huge measured release. Overall, we interpret the analyzed seismic activity in form of bursts as the results of episodes of unsealing followed by temporary sealing of the ground. Strong variations, such as those reported by Federico et al. [3], in temperature and gas/vapor pressure could allowed such processes in the highly fractured rock of the volcano.

Future analyses will consider precision location of the most energetic events, although the common micro-seismicity of the region limits the possibility of application to the majority of the seismic events considered in this study. Finally, our analyses showcased that SOM can be an effective tool to support monitoring and warning of future unrest at Vulcano. We also recommend their application to other volcanoes worldwide for which there is a limited a-priori consolidated knowledge about the typology and origin of the recorded seismic signals.