The 2019 Eruptive Activity at Stromboli Volcano: A Multidisciplinary Approach to Reveal Hidden Features of the “Unexpected” 3 July Paroxysm

Abstract

In July and August 2019, Stromboli volcano underwent two dangerous paroxysms previously considered “unexpected” because of the absence of significant changes in usually monitored parameters. We applied a multidisciplinary approach to search for signals able to indicate the possibility of larger explosive activity and to devise a model to explain the observed variations. We analysed geodetic data, satellite thermal data, images from remote cameras and seismic data in a timespan crossing the eruptive period of 2019 to identify precursors of the two paroxysms on a medium-term time span (months) and to perform an in-depth analysis of the signals recorded on a short time scale (hours, minutes) before the paroxysm. We developed a model that explains the observations. We call the model “push and go” where the uppermost feeding system of Stromboli is made up of a lower section occupied by a low viscosity, low density magma that is largely composed of gases and a shallower section occupied by the accumulated melt. We hypothesize that the paroxysms are triggered when an overpressure in the lower section is built up; the explosion will occur at the very moment such overpressure overcomes the confining pressure of the highly viscous magma above it.

1. Introduction

The Stromboli volcano activity is characterized by continuous degassing accompanied with explosive transients of variable intensity. Different eruptive vents can be the source of the explosions and they often change in shape and number. Occasionally, more violent explosions named “Strombolian paroxysms” interrupt its ordinary activity, causing considerable apprehension among those living in the villages located at the base of the subaerial part of the volcano. These events involve several almost simultaneous craters and erupt volumes of materials much larger and with considerably higher energy than during ordinary activity. Bevilacqua et al. [1] compiled a new historical catalogue of paroxysms that have occurred in the last ca. 140 years. During a paroxysm, two different types of bombs and scoriae can be emitted, filling the entire summit area of Stromboli: the highly crystalline “brown” and the poorly crystalline “blonde”. The two facies, though texturally very different, have similar compositions with both being basaltic and reflecting a peculiar feature often associated with large explosive events [2].

The violent outbursts occurring at Stromboli on 3 July and 28 August 2019 are analysed here as a case study of volcanic activity without the occurrence of evident long- and medium-term precursors. Viccaro et al. [3] consider them exceptional examples of how the volcanic behaviour at open-conduit basic systems may change drastically without apparent notice. Moreover, Giordano and De Astis [4] underline how the 3 July paroxysm was unexpected and preceded only by significant ground deformation precursory phenomena, examined by Giudicepietro et al. [5] and Di Lieto et al. [6], which started just 10 min before the eruption. Andronico et al. [7] evidence the lack of ‘immediate’ precursors, pointing to some longer-term changes, and suggest an urgent implementation of volcanic monitoring at Stromboli to detect such long-term precursors.

Giudicepietro et al. [5] began highlighting the presence of medium-term seismic precursors of paroxysmal activity in seismic data and provided valuable evidence for the development of an early warning system for paroxysmal explosions based on strainmeter measurements. Therefore, starting from this insight and extending the analysed temporal period, the signals recorded by INGV monitoring networks have been re-analysed with the aim of finding more hidden peculiarities in the time series that are able to reveal medium- and short-term precursors of the impending outbursts. Our attention focused on the more sensitive instruments (tiltmeters, deep borehole strainmeters and seismometers) to investigate the days, hours and minutes before the paroxysms and on GNSS data to identify eventual episodes of magmatic intrusion at shallow or deep levels. We also analysed the thermal radiance as observed through satellite images and we performed an analytical approach to study the images collected by the cameras deployed in the summit area of the volcano. Our findings confirm that Stromboli volcano offers only barely detectable signals. The proposed analysis of variations in shape and frequency content of VLP (Very Long Period event), and the trends in tilt and strainmeters time series coupled with the thermal radiance in the summit area of the volcano and a “numerical” approach to visual monitoring of the volcano, represent a real possibility to improve our ability to forecast such dangerous events.

Finally, we performed an analytical inversion of tilt, strain and displacement data variations recorded during the five minutes prior to the explosion to estimate the small volume responsible for triggering the July paroxysm. This source model, along with multiparametric observations, gave us the possibility to define a conceptual model of the behaviour of some paroxysms occurring at Stromboli.

The 3 July and 28 August 2019 Paroxystic Explosions

The paroxysm of 3 July took place at 16:46 (local time). The violent explosion occurred at the southwestern edge of the summit crater area, known as the Ginostra vent. It was preceded a couple of minutes earlier by minor lava overflows at all the vents of the crateric terrace. The eruptive column rose more than 2 km above the summit area and the cloud had a south-west dispersal direction. A close analysis of the explosive column (possibly due to a fortuitous picture in Figure 1) shows that it was composed of two distinct sections: the lower one, dark in colour in which the parabolic trajectories of the big ejecta directed them to fall in the summit area; above this first rose-shaped amount of ejecta it was possible to observe a second one, which surmounted the first one and was characterized by a grey to blonde colour and by higher and more acute parabolic trajectories of the ejecta. The eruptive column underwent a partial collapse towards the Sciara del Fuoco on the north-west side of the volcanic edifice, which produced a pyroclastic flow directed towards the sea.

The paroxysm of 28 August took place at 12:17 (local time). The event featured three distinct explosions; the first two were localized at the southwestern vent (Ginostra vent) and the third one, of minor intensity, occurred at the north-eastern vent. The eruptive column rose about 4 km above the summit craters, then collapsed to generate a pyroclastic flow directed to the sea along the Sciara del Fuoco (Figure 1b).

2. Materials and Methods

A multidisciplinary approach, comprising geological, geochemical and geophysical analyses, is necessary to obtain comprehensive knowledge of the complex phenomena at Stromboli where volcanic dynamics are characterized by small changes hardly detectable by using a monodisciplinary approach. We compared the results obtained from different analyses using most available datasets: the images from remote cameras with an innovative method, the thermal data acquired by the SEVIRI and MODIS satellite sensors and the time series from geodetic data (GNSS, tilt and strain). In order to strengthen the new findings, we compared the seismic data in the VLP frequency range with strainmeter data. Figure 2 shows the monitoring devices used in this work. A full description of data and methods, analysis techniques and additional figures are reported in Appendix A.

2.1. Thermal Observations from Fixed Camera and Satellite

Data acquired by the Flir A320 Thermal Camera installed at the INGV-OE Monitoring Station SPT (Stromboli Pizzo Thermal) located at Pizzo sopra La Fossa, 918 m above sea level (ASL), were processed. We considered the hourly frequency of the explosion and the dispersions of the material involved as follows (see Appendix A for details):

$$C_D={\mathit{F}}_{\mathit{h}}·{\mathit{D}}_{\mathit{h}}$$

(1)

where ${\mathit{F}}_{\mathit{h}}$ is the hourly frequency of the explosions and ${\mathit{D}}_{\mathit{h}}$ is the dispersion (calculated as number of pixels of the thermal anomaly) and C_D is defined as the cumulative dispersion.

The choice is justified by the observation during the long-time intervals of ordinary explosions: average or low-energy events make ${\mathit{F}}_{\mathit{h}}$ and ${\mathit{D}}_{\mathit{h}}$ entirely uncorrelated, while during bigger ordinary events, the system responds by lowering the frequency to maintain its balance.

Accordingly, we plotted the C_D for the period 3 April 2019 to 3 July 2019 (when the camera was destroyed by the paroxysm) related to each area filtered by a simple 49-sample mean average (Figure 3a). Another interesting result obtained from the automatic detection of explosions from videos is the attribution to the different crater area of every single explosion. If we plot both the Cd related to each area of emission (north and south), the results can be summarized as in Figure 3b.

We processed multi-spectral data acquired by MODIS and SEVIRI satellite sensors by using the HOTSAT thermal monitoring system [8,9]. This system allows for the automatic detection of thermally anomalous pixels (hotspot) and quantifies the thermal activity by computing the radiant heat flux. During the first six months of 2019, we find hotspots mainly in MODIS images, since they are more sensitive to lower thermal energy (see Figure A2a, Appendix A), but after 3 July we followed the more energetic phases of the eruption thanks to SEVIRI (see Figure A2b and Figure A3, Appendix A). Figure 3a shows good agreement between the MODIS-derived radiant heat flux during April–3 July 2019 and the cumulative dispersion computed from the ground-based thermal camera. Indeed, the higher values registered by MODIS from June measure the continuous spattering activity, which provides a persistent thermal anomaly in the vent area. The MODIS signals are thus a consequence of the increasing explosive activity that is assessed by the cumulative dispersion.

2.2. Ground Deformation and Seismic Data

Long-term tilt signals (see Appendix A for details) did not show any variation until the end of May when the N275°E component of TDF began a lowering trend (tilt down toward summit area). No evident changes were recorded in the other TDF component or in PLB signals (Figure A4—Appendix A).

Tiltmeters recorded changes before the 3 July 2019 paroxysmic explosion, indicating an uplift of the summit area (Figure 4a,b). Variations are small (about 0.1 microrads) and detected about 2–3 min before with PLB and 3–4 min with TDF. During the two explosive episodes of July and August 2019, tilt showed variations of 1–2 microrads.

No meaningful changes in trend or abrupt variations are evident in the daily time series of global navigation satellite system (GNSS) data (Figure A7, Appendix A), excluding any possibility of deep or shallow intrusion of large dimensions preceding the two paroxysms of July–August 2019. We estimate velocities, uncertainties and noise properties of the four GNSS stations using HECTOR software [10]. A clear difference in terms of velocities is shown by the southern sector of the island (Figure 4d), and this feature is clearer if a “SVIN-fixed” reference system is applied (Figure 4e). The residual velocities of STDF and SPLN station show a contraction between the northern and southern sectors of the island (Figure 4e).

The high rate (1 Hz) GNSS data show large displacements related to the 3 July 2019 paroxysm (Figure A8a, Appendix A). However, they are of unclear origin and possibly due to the large dispersion of ash and hot water vapor during the event that influenced the propagation of the GNSS signal. The 30 s-GNSS data recorded at SPLN station in June 2019 show a clear uplift between days 23 and 28 (Figure 4f), indicating a possible shallow propagation of volcanic fluids very close to the southern edge of the island just a few days before the 3 July explosion. This variation, although of limited magnitude, is an anomaly in the high frequency time series of the SPLN station. Moreover, SPLN is the only GNSS station that recorded a variation in this period, probably due to its position very close to sea level and in the southern sector of the volcano.

Before the 3 July paroxysm, the SVO strainmeter recorded intriguing signals (Figure 4c): comparing the precursory effusive phase recorded by the monitoring camera installed near the southern vents of the main crater, the minimum strain value (frame 14:34:32) occurs about 11 min before the paroxysmal explosion (frame 14:45:40) [5].

The two paroxysms that occurred in 2019 have a similar strain time history: Di Lieto et al. [6] show the details of the strain recorded a few minutes before both events and evidence a correlation among them and the paroxysm that occurred in 2007, suggesting a common source mechanism. During the months preceding the July paroxysm, no significant changes in the unfiltered strain data were observed.

Fast rising within the conduit of gas bubbles generates seismic signals in the shallow part of the feeding system [11] in the frequency band 2–50 s, which are also visible on tiltmeters and strainmeter data, as a transient of about 50–200 nrad [12] and 2–4 nε, respectively. Hence, we focused on the VLP frequency band. We will refer to strain transients recorded in this band as “strain VLP”. From a statistical analysis conducted on almost 100,000 strain VLPs recorded in the period June 2018–December 2019 by the SVO strainmeter, we found that all the examined strain VLP parameters routinely analysed by the scientific community on seismic signals change slightly about 30 days before the 3 July event, showing a weak correlation once they are plotted together (Figure 5b–f). Indeed, clearer variations occur in strain VLP shape, which abruptly changes on 3 June and then again on 25 June (Figure 5a). The significant number of strain VLPs belong to two different families (Figure 6a): the first family contains events characterized by a pronounced quasi-symmetric positive, dominant, bell-shaped peak with an oscillation before and after the main peak, clearly visible from 3 June, in contrast to the second family, whose predominant characteristic is a longer oscillation without a prevailing peak, which appears from 25 June as the first one suddenly disappears. The first family consists of a dilatational strain transient with durations of 30 s, while in the second family dilatational and compressional phases are more balanced, suggesting a variation in the visco-elastic medium characteristics, with the same values for their amplitude and duration.

The 28 August paroxysm was not preceded by changes in any parameter of the strain VLPs. It occurred during an active phase, with all the monitored parameters already modified with respect to the normal activity and with a very high amplitude of the volcanic tremor [5].

2.3. Seismic and Strain Data Comparison

In Giudicepietro et al. [5], seismic data (with a focus on volcanic tremor, explosion quakes, volcano–tectonic and VLP events) from November 2018 to September 2019 are analysed, showing the presence of significant variations through the whole month preceding the paroxysm of 3 July. Introducing a new parameter called “VLP size”, differences in VLP shape recorded by the seismic network were observed by Giudicepietro et al. [5], who linked temporal evolution of explosive source to a higher gas content in the Strombolian explosive activity, and also by Chouet et al. [11], who found a correlation with activity at the different vents.

In line with Giudicepietro et al. [5], we performed a cross-correlation analysis of a selected dataset of VLP recorded by the seismometer installed at STRA station (hereinafter “seismic VLP”), which is located 500 m away from the vents at about 850 m ASL, to verify the changes in strain VLP shape. Based on strain VLP event arrivals, seismic waveforms were triggered and extracted from the vertical component of seismic recordings in the time interval 1 June–3 July 2019. We collected about 2400 seismic VLPs filtered in the frequency band 2.5–20 s. With the aim of investigating similarity and waveform variation over time, seismic traces were cross-correlated with each other. Choosing 40-s-long windows and a cross-correlation coefficient threshold of 0.8, seismic traces were grouped into 11 families. Only families consisting of at least 15 events, whose stacked and normalized waveforms are shown in Figure 6b, were selected. The time variation of the daily occurrence rate of events of each of the six families selected has also been carried out (Figure 6c). Results show an increase in the event number starting from 25 June and indicate that family 1 is the most populated one, and seismic VLP events belonging to this family were recorded during the whole period (Figure 6c). The main variation associated with this event type is a decrease in number occurring around 22 June, which was then followed (on 27 June) by an increase until the paroxysm of 3 July. Furthermore, after the family 1 occurrence rate decreases, seismic VLP waveforms were grouped in several families (Figure 6c). After around 22 June, all the different seismic VLP shapes found by cross-correlation were observed. We highlight that seismic VLP types belonging to families 3 and 4 have similar waveforms with respect to the most frequent ones (Figure 6b). After 25 June, beside seismic VLP belonging to family 1, only those of families 2, 5 and 6 are observed (Figure 6c). Visual inspection of waveforms allowed us to identify similarities between the seismic VLP event shape of family 1 (Figure 6b) and the shape of the first kind of strain VLP (Figure 6a) and between the seismic VLP event shape of families 2 and 5 (Figure 6b) and the shape of the second kind of strain VLP (Figure 6a). It should be noted that while in seismic traces seismic VLP events belonging to family 1 have been pinpointed throughout the period, strain VLP events of kind 1 disappear after 25 June. The slight differences between seismic families contrast with the clear strain VLPs’ shape changes. Therefore, we focused on seismic VLP events belonging to family 1, split the time interval into two periods, 1–24 June and 25 June–3 July, and computed stacks. Stacked and normalized waveforms of the two periods, reported in Figure 6d, highlight a slight modification of seismic VLP: in the second period, seismic VLP event-stacking is characterized by a longer tail, suggesting decreased damping of seismic VLP events, whereas the frequency tends to remain constant.

The comparison between results obtained by investigating strain VLPs with seismic VLPs confirms a change in VLP shape, and as a consequence in the source mechanism and/or in the fluid properties after 25 June. In addition, seismic VLP component analyses evidenced that this process started a few days earlier and was probably complex, indeed, starting from 22 June and lasting until 25 June, seismic VLP event types belonging to families 3 and 4, which exhibit waveforms similar to those of the seismic VLP of family 1 (Figure 6b–d), were also recognized.

3. Modelling the Volcanic Source

We performed an analytical inversion of tilt, strain and displacement data variations recorded during the five minutes before the explosion (see Appendix A.4. for more details). We used a pressurizing source embedded in a homogeneous half-space described as a simple spherical source [13]. In particular, all available data have been inverted simultaneously using the software developed by Cannavò [14]. The effects of the topography were taken into account by using the method of Williams and Wadge [15]. We obtained a pressure source located under the craters area at about one hundred meters above sea level (Figure 7). The source is characterized by a volumetric expansion of about 10,500 m³. This value found for the source of the explosion is much lower than that evidenced by Di Lieto et al. [6]. This is due to the fact that the model presented in the current manuscript takes into account all available ground deformation signal changes occurring in the five minutes prior the paroxysmal event to estimate the small volume responsible for triggering the explosion. Instead, Di Lieto et al. [6] consider the strain signal variation recorded during the phenomenological phase of the paroxysmal activity that involves both the middle (about 1.5 km b.s.l) and shallower storages. The volume of the pressure source (about 10⁶ m³) is calculated by considering the peak to peak strain variation over the entire signal duration. The different modelling approach used in Di Lieto et al. [6] also conditions the source depth estimate, which has been found to be consequently lower in the present work. Moreover, it is noteworthy that the Mogi’s model [13] does not provide a correct solution for the pressure source located above the elevation of recording stations and, therefore, we limited the maximum values that the “Z” coordinate can assume to the minimum elevation of our stations (about +80 m above sea level). Our solution achieves this value. Despite this limitation, our approach data fit is good (Figure 7). In fact, we obtained a reduced chi-squared χ² equal to ~1.0 considering an a posteriori standard deviation of 0.01 m and 0.02 m for the horizontal and vertical displacements, an uncertainty of 1.0 μrad for the tilt data and of 0.02 με for the strain data. Seeking to overcome the intrinsic limitation of the Mogi [13] approach, if we exclusively invert the tilt variations deduced by the three summit seismic stations (STRA, STRE, STR1, Figure 2), we obtain a pressure source that is not well-constrained located at an altitude of (+480 ± 620) m with respect to sea level. Therefore, the good result obtained using all available data leads us to believe that the real pressure source representing the five minutes preceding the eruptive explosion is probably not much above the obtained level of +80 m. The estimated small volume variation does not affect the signals recorded at the high-frequency GNSS network. Conversely, it was visible to the borehole instruments due to their high sensitivity. Moreover, it is well known that geodesy alone cannot constrain important parameters such as chamber volume or pressure and that deformation studies resolve exclusively for a volume variation [16] that cannot be compared to the erupted volume. However, the small value that we obtained for the volume variation preceding the paroxysm could suggest a key role played by the continuum gas phase with respect to a supply of new magma, as explained in the following paragraphs.

4. Discussion

The intermittency of Strombolian explosions has been interpreted by Jaupart and Vergniolle [17,18] as a consequence of the periodic collapses of a foam layer formed by the accumulation of bubbles that remain trapped below the magma surface. According to this view, the maximum packing condition is reached when the contact area between bubbles becomes equal to the cross-sectional area of the bubble. Similarly, Chouet et al. [11] proposed a source mechanism on the basis of seismic broadband data and the observed intermittency of explosive eruptions, which point to a degassing process in which gas bubbles rising in the inclined liquid-filled dike tend to concentrate against the upper wall of the dike, thus enhancing the collisions between bubbles and their coalescence. As demonstrated in laboratory experiments with analogue fluids [19], localized bubble concentration against inclined surfaces makes slug flow unavoidable at all but the lowest gas-flow rates for non-vertical conduits. Therefore, a crucial element in this degassing process is the inclination and the shape of the conduit.

Recently, Oppenheimer et al. [20] have proposed, on the basis of analogue laboratory models, the important role of near-surface (down to 800 m) crystallization and the variations of the crystallinity and the interactions of crystal-bubbles in regulating the intensity of degassing and explosive activity. These authors proposed a “weak plug” model for Strombolian explosions, evolving from low viscosity style towards more crystalline, a stronger and less permeable plug corresponding to larger events. These last events are characterized by an increase in the crystallinity, determining a different conduit condition and larger size explosions. This model predicts some features of the Strombolian explosions, such as the variability of their sizes, duration, pulsation and fountaining according to the degree of near-surface crystallization, but suffers from a lack of a quantitative explanation of geophysical signals. Moreover, in Caracciolo et al. [21] petrological evidence for a ‘soft plug’ model based on observations of earlier explosions are discussed.

To find an eruptive mechanism that could take into account the large spectrum of geophysical data, the direct observation of the paroxysm and the nature of the erupted material, we decided to adopt the “push and go” model proposed for the feeding system of Mount Etna volcano in Ferlito [22]. The magma is described not as molten rock residing at a certain depth within the volcanic edifice or in the crust, but, on the contrary, as a flux of fluid continuously moving upwards, in which two vertically distributed portions can be recognised. Below is a solution made of ~70% in volume by a continuum gas phase (mostly H₂O) in a supercritical state (density 360 kg/m³) and 30% in volume by basaltic components dissolved in it, named water melt solution (WMS); the overall density is low (1140 kg/m³). At the top of the WMS, the continuous loss of gas, characterising all persistently active volcanoes, will leave a highly dense (2800 kg/m³) basaltic melt, defined as a continuous melt phase (CMP). The CMP has a low enough viscosity to permit the continuous passage of the gas bubbles coming from the degassing of the WMS (Figure 10a). The transition between the WMS and the CMP can be erratic within the plumbing system. In this new paradigm of magma, the eruptions can be considered to occur due to a significant pressure increase within the WMS in the deep plumbing system, which cannot be kept confined by the weight of the overlying basaltic melt (CMP). This process can be viewed like a piston acting within a cylinder, in which the overpressure and expansion of the gas-rich WMS appears crucial, able to push the high-density basalt of the CMP up to the surface effusion. There are two main causes for the overpressure within the WMS: an increase in the gas flux from a deeper source or a decrease in the permeability of the CMP, which in turn could be caused by a viscosity increase in the CMP due to the cooling of the melt and incipient crystallization (see [20]). The eruption will be effusive, with lava flows, if the CMP is able to dynamically contain the expansion of the WMS below; whereas, if the CMP is unable to confine the WMS expansion, the latter will be violently emitted and the resulting eruption will be an explosive paroxysm.

From a joint comparison among all data, we identified two different volcanic state phases: a precursory phase, starting in early June 2019 (yellow and blue filled rectangles in Figure 8) that consists of a medium and a short pre-paroxysmal time, and a paroxysmal phase, starting on 3 July (red filled rectangle in Figure 8), that includes all the eruptive volcanic activity, characterized by effusive activity and by a second paroxysm on 28 August.

4.1. The Precursory Phase

From June 2018 to the end of April 2019, we do not observe any significant variation in the measured parameters. At the beginning of June there is a significant increase in the cumulative dispersion parameter (considering the medium value of the preceding months), which is directly connected to the energy of the explosive events, and of the number of strain VLPs, whereas GNSS does not show any variation. MODIS data revealed almost continuous thermal anomalies from June 6, showing a radiant heat flux curve with an increasing trend reaching a maximum value of 0.45 GW on 2 July at 12:00 GMT. Tilt evidences a trend change on the TDF N275°E component, which, however, has no counterpart on the other components or the PLB data (Figure 8), suggesting a local source for this variation. Such a lack of “depth” related signals indicate that this phase of unrest must be associated with the uppermost portion of the feeding system, specifically the CMP. Considering that the cumulative dispersion parameter indicates an increase in the magnitude, not in the number, of the explosive events that are caused by the breaking of gas bubbles at the surface of the CMP, it must be deduced that more energetic explosions must be generated by larger bubbles or by bubbles that have a higher gas overpressure. The large rose-shaped incandescent ejecta observed during this phase (Figure 1) would suggest the first hypothesis is more plausible, and large gas bubbles are more likely. Since the most efficient mechanism to form bubbles within a basaltic melt is coalescence, it must be inferred that the viscosity of the melt had to be sufficiently low to allow congruous bubbles mobility (Figure 10a). It can also be deduced that the flux of hot gases passing through the CMP was high enough to maintain a thermal state able to compensate the heat dispersal and therefore hinder crystal nucleation and growth, thus retarding an increase in the CMP viscosity (see Ferlito [22], p. 19). In this first phase, we can reasonably say that the gas flux and therefore the WMS flux in the Strombolian feeding system was high; unfortunately, the lack of a gas sensor network, such as the one set up at Mount Etna, which has allowed us to explain the contribution of the gas flux to the 28 December 2014 eruption (Ferlito et al. [23]), prevents us from making a similar quantification for Stromboli. Indirect confirmation is reported by Inguaggiato et al. [24,25], who found an increase in CO₂ soil emissions starting from June, weeks before the paroxysmal event. These data also show higher variance with respect to previous recordings, which index the instability of the volcanic activity.

A similar indication is provided by the increased number of strain VLPs belonging to two main families (Figure 8 and Figure 9). It should be noted that in the first phase, until the first days of June, the feeding system was in an ordinary state and the boost of the gas flux was somehow accommodated by a corresponding increase in the gas flux from the summit craters without any evident change within the feeding system. This first period was followed by a significant variation in the radial tilt and by a sudden appearance of family 1 strain VLPs, which persists until 25 June when family 2 becomes predominant. Waveform analysis of seismic VLP events confirmed the changes occurring in this period (Figure 6). Seismic VLPs at Stromboli are interpreted as being due to pressure change caused by the generation or the ascent of a gas slug in the volcanic conduit [11,26]. Their waveform similarity has been studied in the literature in order to investigate their characteristics in terms of steady activity and changes related to eruptive activity [27,28,29]. After 25 June, several families exhibiting different waveforms appeared (Figure 9) and waveforms belonging to family 1, recognized throughout the period (1 June–3 July), underwent modification in oscillation duration (Figure 6b), suggesting changes in the source mechanism or in the damping property of the medium.

As is known from simplified models of seismo-volcanic sources [30,31,32], the resonance frequency and damping of the system is strongly influenced by the nature of liquid and gas content. Moreover, since seismic VLPs are linked to inertial displacement of material in which they propagate [33,34], the different characteristic waveform of two types suggests two origins: a variation in mechanical magma properties or a different location. In Giudicepietro et al. [5], all seismic VLPs show little variation in incidence angle and azimuth and, starting from 25 June, become very localized. The locations do not show remarkable variations before or during the eruptive phase of summer 2019. These observations seem to reinforce the first hypothesis. In contrast, an explanation of strain VLP shape characteristic is not available in the literature: this work is a starting point in understanding the volcanic processes on which these observations are based.

The variation of these parameters indicates that during this phase, starting on 25 June, the volcanic edifice was reacting to an overpressure of the WMS that the system was unable to discharge through the CMP and the summit craters (Figure 10b). A question that needs to be addressed is whether this different behaviour can be related to a decreased permeability of the system or to an increase in the WMS flux. There is no easy answer, but the products emitted during the short paroxysm of 3 July, in spite of a similar chemical composition, had two distinct textures: the partly crystallized “brown” (Figures S1 and S2) and the totally glassy “blonde” (Figures S3 and S4). The crystal content of the brown indicates that before the paroxysm the CMP must have reached a considerable degree of viscosity, which in turn could have decreased the permeability. Moreover, during the increase in the WMS flux, occurring in the first phase, only the gaseous component was emitted and therefore discharged from the feeding system, whereas the basaltic (melt) component was mainly retained within the uppermost portion of the feeding system, increasing the size of the CMP up to the point in which the gas flux was unable to provide the heat sufficient to maintain low viscosity to hinder the crystallization of the CMP. The decreased CMP permeability delayed the flux of the gas bubbles through it and must have therefore led to an overpressure of the WMS. Such an isotropic pressure must have acted against the entire upper portion of the feeding conduit, boosting a deformation that can be considered responsible for the radial tilt variation observed at the TDF station (see Figure 8). Now, any volcano being a complex dynamic system, it appears obvious that such a condition of instability is temporary and must evolve somehow. Unmistakable signs of the dramatic evolution began to show from 25 June, one week before the paroxysm. The number of VLPs increased significantly and so did their intensity. At the same time, the energy of the explosions had a conspicuous boost. The radial tilting (TDF N275°E component) increased notably and the SPLN GPS station showed an uplift (Figure 4f). This further rise of the physical signals was the prelude of the paroxysm and can be interpreted as the last bearable increase in WMS overpressure before the break-up of the CMP.

The 3 July paroxysm that occurred at 14:45:45 UTC was the sudden and violent emission of the CMP unable to contain the dramatic pressure increase within the underlying WMS, which produced an inflation of 10,500 m³ and was preceded, about one minute earlier, by a shallow (800 m deep [5]) volcano tectonic event (VT), suggesting a fracture in the base layer of the CMP column. According to the proposed “push and go” model, the explosion would occur at the very moment at which the overpressure in the WMS overcomes the confining pressure of the CMP above it (Figure 10c). We maintain that our modelled pressure source is roughly located at the interface between the WMS and the CMP. Exclusively for speculative reasoning, we would like to verify how large the CMP would need to be to keep the underlying layer (WMS) confined by its own weight. Considering that the pressure source is located at about 700 m below the surface, a hypothetical cylindric conduit representing the CMP and having a volume equal to the estimated “ΔV”, to roughly balance the pressure of the modelled source, should have a radius of about 2 m. This value is compatible with the dimensions of the open vents usually visible at the summit craters of Stromboli.

4.2. The Paroxysmal Phase

What might be of further interest is examining the 75 min preceding the paroxysm. In Figure A6, from minute −90 to minute −75, the only notable variation is given by the strain VLPs occurring at about a 5 min intervals. From minute −75 to minute −60, the positive strain has an increase that reaches a peak of 0.55 × 10⁻⁸. This positive strain could have been caused by an overpressure likely due to a boost in the WMS flux. In the following 50 min, the excess flux was discharged, passing through the overlying CMP and thus significantly decreasing the positive strain. Two distinct patterns of strain release are notable (Figure A6): at the beginning the rate was fast, in the first 5 min it decreased by 0.2 × 10⁻⁸ with no strain VLPs associated; whereas from minute −55 to minute −10, the strain decreased with a slower rate from 0.2 to −0.25 × 10⁻⁸ and in a discontinuous “step way” fashion, in which the occurrence of each strain VLP event is associated with an immediate step of strain decrease. This information suggests that the strain VLPs generate at the interface between the WMS and the CMP are probably associated with the transition of WMS from a continuous gas + melt flux to a discrete flux of gas (bubbles) within the melt. Two such patterns of pressure release can be interpreted as the contemporaneous existence of two distinct depressurization modes: a “continuous” mode, in which the pressure decrease is due to lava emission through the summit vent and/or gas discharge through the articulated network of fractures that opens at the summit of Stromboli; and a “discontinuous” mode, in which the pressure in the WMS is released through the sudden passage of the gas through the interface with the CMP.

Finally, the last 10 min before the paroxysm are particularly interesting: the strain increased one order of magnitude in less than 10 min; such a sudden and fast increase is not easily explainable as the result of an augmented viscosity of the CMP. What can be hypothesized instead is a dramatic increment in the WMS flux, which was too fast to be released through the CMP. In fact, the travelling time for the gas bubbles to pass through the CMP is in the order of a few minutes, but if the flux increase is faster the system cannot release the gas overpressure in the WMS. The lack of strain VLPs just before the paroxysm is therefore indicative that the gas of the WMS was not able to pass through the CMP, but instead pushed its entire mass (the brown) upwards that was entirely erupted in a few seconds along with the WMS that constituted the glassy texture (the blonde).

The violent emission of the CMP and the WMS occurring on 3 July has strongly modified the uppermost portion of the feeding system. A consideration needs to be made: the very fact that after the sudden emission of CMP and WMS there was no collapse of the summit edifice indicates that the material that was expelled was immediately reintegrated by the rise of deeper melt. Indeed, the eruption continued after the paroxysm with the emission of a lava flow and a regular explosive activity from the summit craters, indicating that the conduit was filled up again with the CMP. A final consideration needs to be made on the paroxysm that occurred on 28 August (Figure A6): the strain curve in the 90 min preceding this second paroxysm appears to be slightly similar to the one of 3 July, but is much smoother and more importantly there is a lack of significant strain VLPs.

5. Conclusions

The paroxysm of 3 July 2019 prompts a general reflection: although an impressive number of monitoring devices have been installed for many years on the island, such a paroxysm occurred without any meaningful variation of the geophysical parameters used for the evaluation of the activity of the volcano (see, e.g., [4,5,7,24]). Volcanic tremor, number of VLP events, ground deformations and earthquakes all showed small unclear variations: we have attempted here a careful revision of “hidden” characteristics of these signals to propose a set of parameters able to detect a change in the state of the volcano, which might lead to unrest or even a paroxysm, so as to help prevent dangerous consequences. However, the parameters considered during the monitoring of the activities of Stromboli volcano is not the only aim of our research: indeed, a step beyond needs to be taken in order to understand why the monitoring of major or paroxysmal events at Stromboli volcano is so elusive and sometimes difficult to interpret. Small variations observed in the geophysical signals recorded in Stromboli could be related to two factors: the small quantities of volcanic fluids involved in Stromboli and the role of gas in their evolution. Hence, we apply a conceptual model previously proposed for Mount Etna volcano (the so called “push and go”) that fits the observed variations occurring at medium and short term before the paroxysm. In our opinion, a high rate of (N > 100/day) and predominance of family 2 strain VLPs over family 1 strain events indicates a high likelihood for an explosive event within the next week to 10 days. In this sense, this information could improve our capability to forecast dramatic events such as the one that occurred in July 2019. Regarding the model, we are aware that, from now on, a more detailed investigation must be carried out at Stromboli and at other basaltic volcanoes to verify its limits, its possible applications, and to investigate some aspects that are still unclear.

Performed data analyses and integration of proposed parameters, estimable/retrievable from real- or near real-time data, helped us shed light on the processes leading to the 3 July Stromboli paroxysm. The following parameters and their variation over time, could be helpful to evaluate the volcanic hazard from explosive activity at Stromboli:

• Continuous and automatic evaluation of the cumulative dispersion (Figure 8) from video recordings of the summit explosive activity in Stromboli, because the mere observation of the number of explosions occurring from the Stromboli craters is not sufficient and a proxy of the energy involved in the explosive activity is necessary. This value should be coupled with the satellite-derived radiant heat flux to measure the thermal energy due to the extension and persistence of both the vent area and hot deposits as viewed from space.
• Classification of strain and seismic VLP waveforms in different “families”. The transition from one VLP family to another, or the superposition of several seismic VLP families, can be an indicator of changes in the fluid properties, such as the change in permeability of the higher portion of the magma in the main conduit and this can be considered an alteration of the normal condition leading to the mild ordinary explosive activity.
• Small variations in long/medium term of tilt and GNSS time series, coupled with thermal observations from satellite, must be considered the complementary data whose co-variation can confirm the occurrence of a medium term (one–two weeks) alert of potentially impending large explosive activity.

For these reasons, the present work introduces a new approach to joint multiparametric analysis of monitoring system data. Soon, it could be considered a prototype of an early warning system useful for civil protection purposes on a time scale in the order of weeks/days before a sudden change in dangerous volcanic phenomenology.