Chronicle of Processes Leading to the 2018 Eruption at Mt. Etna As Inferred by Seismic Ambient Noise Along With Geophysical and Geochemical Observables

This work analyzes temporal variations of seismic velocities at Mt. Etna from August 2018 to February 2019. During this time period, a strong flank eruption accompanied by intense seismicity and ground deformation took place along a fracture that opened on 24 December 2018 at the base of the New South‐East summit crater. Furthermore, two moderate earthquakes—the 6 October 2018 ML 4.7 and 26 December 2018 ML 4.8—associated with the volcanic activity were recorded. In this study, we computed cross‐correlation functions (CCFs) between windows of seismic ambient noise to identify seismic velocity variations within the volcano edifice. We calculated daily CCFs at 16 stations for 120 interstation pairs using the vertical component in the 1.0–1.5 Hz frequency band. We observe that dv/v starts to decrease rapidly from the beginning of October 2018 and reaches approximately −0.45% in the pre‐eruption period. The spatio‐temporal distribution of seismic velocities shows that the reduction of dv/v mostly occurs in the vicinity of the summit and close to the flank area and is interpreted to be affected by magmatic intrusion at 0–3 km depth. To infer the source mechanism of this eruption, we compared these observations with volcano‐tectonic earthquakes, volcanic tremor, volcanic degassing, gravity, and ground deformation data. Our study suggests that a relationship between magma intrusion and associated crack opening is responsible for the decrease of dv/v.

Seismic noise interferometry was previously studied at Mt. Etna by Cannata et al. (2017) Cannata (2012) by using the coda wave interferometry applied on repeating volcano-tectonic earthquakes (VTs).
In this work, we apply the seismic ambient noise method to Mt. Etna in the period from August 2018 to February 2019 during which the aforementioned strong flank eruption occurred on 24 December 2018, accompanied by the most intense swarm of VTs recorded during the last 20 years (Figure 1 and Figure S1 in Supporting Information S1). To infer the source mechanism of this variation, we integrated seismic ambient noise results with those of VTs, focal mechanisms, volcanic tremor, gas geochemistry, gravity, and ground deformation data recorded by the multi-parametric monitoring network of Istituto Nazionale di Geofisica e Vulcanologia (INGV; Figure 1a).

Volcanological Framework
During the last decade, Mt. Etna has been characterized by frequent eruptive activities, involving most of the summit craters along with lateral fissures (Figure 1a). In particular, from 2011 to 2013, the New South-East Crater (NSEC) produced more than 40 explosive eruptions culminating in lava fountaining (e.g., Behncke et al., 2014;De Beni et al., 2015). From mid-December 2013 to the end of 2014, the lava fountains at NSEC were replaced by violent Strombolian activities mainly at NSEC (e.g., De Beni et al., 2015) as well as by an eruption from a fissure opened at 3,000 m elevation on the lower eastern flank of the North-East Crater (NEC) during July-August 2014 (e.g., Spina et al., 2017;Viccaro et al., 2016). On 28 December 2014, another lava fountain was produced by NSEC (Gambino et al., 2016). In December 2015 and July 2016, violent explosive activities affected also Voragine Crater (VOR; Cannata et al., 2018). In 2017, the activity became again more effusive with a series of weak Strombolian episodes and lava flow emissions between February and April (Viccaro et al., 2019).
During the first months of 2018, the eruptive activity of Mt. Etna was characterized at NSEC and NEC by mild ash emissions that intensified turning into a Strombolian activity in July. In August, the Strombolian activity involved the Bocca Nuova (BN), NEC, and NSEC, and on 23 August about at 18:00 UTC, the last crater started an effusive activity that gave rise to modest lava flows on the Eastern flank of the NSEC cone (Crafford & Venzke, 2018). The Strombolian activity from BN, NEC, and NSEC and the discontinuous effusive activity from NSEC carried on up to December (Borzi et al., 2020). On 24 December at 11:11 UTC, an eruptive fissure opened at the eastern base of NSEC and propagated in the direction SSE ( Figure 1a). This eruption, characterized by both explosive and effusive activities accompanied by significant ground deformation and seismic activity, ended on 27-28 December (e.g., Bonforte et al., 2019;Cannavo' et al., 2019;Carbone et al., 2019;Chauhan et al., 2020). This eruptive episode was also accompanied by significant variations in CO 2 and SO 2 fluxes and 3 He/ 4 He ratios (e.g., Gurrieri et al., 2021;Paonita et al., 2021). The following weeks were characterized by mild ash emissions from the summit craters.

Seismological Settings
Mt. Etna is located at the front of the Apennine-Maghrebian Chain on an Early-Middle Pleistocene foredeep clayey succession (e.g., Branca et al., 2011;Monaco et al., 2005). The complex interactions between volcano dynamics and tectonics gave rise to several active faults, causing moderate earthquakes with M L < 5.0 (Figure 1b; Alparone et al., 2015Alparone et al., , 2020. The well-known faults structures in the region associated with the 2018 eruption are (Figure 1b;De Novellis et al., 2019): (a) Pernicana Fault (PF) lies on the volcano's north-eastern flank with an oblique left-lateral strike-slip mechanism (e.g., Azzaro et al., 2012;Bonforte et al., 2007); (b) Ragalna Fault (RF), mainly characterized by a right-lateral strike-slip mechanism, is located on the southwest flank (e.g., Azzaro et al., 2012;Rust & Neri, 1996); (c) Fiandaca Fault (FF) is a part of the main Timpe system in the southeast and shows the dextral oblique-slip behavior similar to RF (e.g., Azzaro et al., 2012;Villani et al., 2020).
The Mt. Etna eruption in 2018 was accompanied by especially the intense seismic activity (Figure 1b; Alparone et al., 2020). The swarm of VT earthquakes could give clues about the magmatic intrusions. The two largest earthquakes that took place during the analyzed time interval occurred on 6 October 2018 with M L 4.7 and on 26 December 2018 with M L 4.8. The former occurred in the south-western flank close to the RF at 6 km depth. The latter was located in the south-eastern flank at a very shallow depth (1.2 km) along the FF (Figure 1c) and caused a significant damage along the surface rupture (Civico et al., 2019;Villani et al., 2020). To study the time-space evolution of the seismicity and its relationship with the seismic ambient noise results, we used the seismic events relocated by Alparone et al. (2020) by double-difference (Zhang et al., 2009) and the 3D velocity model by Alparone et al. (2012). This catalog contains approximately 1,700 earthquakes that took place during the interval 6 October to 26 December 2018 (Figures 1b and 1c, and Figure S1 in Supporting Information S1). Most of the seismicity occurred in a very short time interval, accompanying the 24 December intrusion and clustering in two rock volumes along the eruptive fracture with epicenters aligned in a N-S direction and slightly southeast of the summit (Figure 1b; Alparone et al., 2020). It is also worth noting that earthquakes predominantly occurred in the upper ∼6 km b.s.l. of the Mt. Etna (Figure 1c and Figure S1 in Supporting Information S1). For a better visualization of the seismicity in the region for that period, a 3D animation can be found in Movie S1 in Supporting Information S1.

Seismic Ambient Noise
In this study, we used continuous data from 16 stations equipped with broadband three-component Trillium seismometers (40s cut-off period, Nanometrics™), operated by INGV-Osservatorio Etneo (OE; Figure 1a and Figure  S2 in Supporting Information S1). We processed their vertical component using the MSNoise package (Lecocq et al., 2014). The raw data were band-pass filtered between 0.01 and 8.0 Hz, down-sampled from 100 to 20 Hz, and de-trended. Daily traces were segmented into 1-hr windows, winsorizing (clipping) at 3 Root Mean Square (RMS), and spectral whitening in the 1.0-1.5 Hz frequency band were applied before cross-correlation (Bensen et al., 2007;Lecocq et al., 2014). This frequency band was chosen after a sensitivity analysis, which showed that this band gives the most stable results in terms of both dv/v time series and sharper variations focused on the main eruptive events. In the analysis, the variations of dv/v were computed for other four frequency bands: 1.0-2.0, 1.2-1.6, 1.2-1.4, and 1.4-1.6 Hz ( Figure S3 in Supporting Information S1). Seismic velocity changes were also calculated for the lower frequency bands (0.1-1.0 and 0.5-1.0 Hz) and higher frequency bands (1.6-1.8, 1.6-2.0, 1.8-2.0, and 2.0-2.5 Hz) and further significant variations were not found. Hourly windows were then stacked The focal mechanisms are compiled from different sources as follows: orange, yellow, blue, green, and gray beach balls are taken from GCMT, Pondrelli (2002), Scarfì et al. (2013), Scognamiglio et al. (2006), and GEOFON, respectively. The black mechanisms were calculated in this study. (c) The depth cross-section of the VTs with respect to longitude.

of 16
for each day and cross-correlation quality and stability were analyzed in both the time and frequency domains by visual check ( Figure S4 in Supporting Information S1). The resulting daily CCFs were then stacked using the time-frequency phase-weighted stacking method (tf-PWS; Schimmel & Paulssen, 1997) for 5, 10, and 30 days for each station pair to increase the signal-to-noise-ratio.
Daily variations in seismic velocity ( Figure 2a) were computed using the Moving Window Cross-Spectrum analysis (MWCS; Ratdomopurbo & Poupinet, 1995) between a daily and a reference CCF, which comprises the January 2018 to July 2018 time period in our study and corresponds to the pre-eruption stage. Unlike the stretching technique of Sens-Schönfelder and Wegler (2006), the MWCS method only uses phase spectra to estimate seismic velocity changes, which reduces the contamination due to temporal variations in the source of seismic noise on the seismic velocity changes (Zhan et al., 2013). We compute time shifts, dt, in the coda of daily CCFs relative to the reference CCF. Time shifts and coherency c between the reference and daily CCFs are calculated after the 1 km/s arrival in the coda in 30 s windows and the dt measurements are kept when the coherency c is equal or higher than 0.5. In the sensitivity analysis, the coherency c was changed from 0.5 to 0.8 and any notable alteration for the further computations was not observed ( Figure S5 in Supporting Information S1). The regressing time shifts dt from each window in the causal and acausal part of the coda then provided the daily time shift dt/t. Assuming a homogeneous change throughout the sampling medium, the relative difference in travel time dt is related to the change in the seismic velocity dv as −dt/t = dv/v. The resulting time series for each possible pair of stations on Etna volcano were finally combined in a single median variation of seismic velocity for the network ( Figure 2a). Furthermore, to map the "epicentral area" showing the strongest change in seismic velocity, the method of Brenguier et al. (2014) is applied. According to this approach, a velocity change value, computed as the mean of the velocity changes measured for the pairs with which a station was involved, is assigned to each station ( Figure 3). To estimate the depth sensitivity of the analysis, we calculate the partial derivatives of the fundamental-mode Rayleigh-wave phase velocities (Herrmann, 2013). The sensitivity kernels show that dv/v measurements in our study indicate velocity variations approximately within the upper 3 km of the crust ( Figure  S6 in Supporting Information S1) considering the velocity model of Aloisi et al. (2002).

Volcano-tectonic Earthquakes
We obtained moment tensors (MTs) for the events M L > 3.9 based on a Bayesian bootstrap probabilistic joint inversion method by using the Grond software (Heimann et al., 2018;Kühn et al., 2020). The method provides the uncertainties for the model parameters of resulting moment tensor inversion (e.g., depth) and has been applied to different regions recently-for example, in Southeastern Anatolia , Zagros Simply Folded Belt , and AlpArray seismic network (Petersen et al., 2021). The MTs were calculated via waveform inversion in both time and frequency domains using three-component data (R, T, and Z). Synthetic seismograms were computed using pre-calculated Green's functions (Heimann et al., 2019) based on the combination of two velocity models for the crust (Aloisi et al., 2002) and upper mantle (AK135, Kennett et al., 1995). We used the frequency band of 0.03-0.08 and 0.02-0.05 Hz depending on the source-to-receiver distance for events smaller than M L 4.3 and the events larger than M L 4.3, respectively. We retrieved the MTs in full waveform inversion for the two largest earthquakes (M L = 4.7 and 4.8) and the double-couple solutions for the smaller earthquakes ( Figure 1b and Table S1 in Supporting Information).

Volcanic Tremor
Volcanic tremor is one of the most important and common precursory phenomena of eruptive activity because of its strict relationship with fluid flows through the volcano's feeding system (Chouet & Matoza, 2013). This is particularly true at Mt. Etna where important variations of this seismo-volcanic signal have been observed, both in the medium and short terms, before the onset of eruptions (Cannata et al., 2018;Cannavo' et al., 2019). Volcanic tremor is commonly addressed to as RMS amplitude of the seismic signal measured over a certain time window and frequency band. In this work, we calculated the daily average of the RMS amplitude of seismic signal in velocity (m/s), recorded at a seismic station, ECPN, located at about 1 km from the summit craters ( Figure 1a). We chose to filter the seismic signal in the frequency band 0.5-2.5 Hz as it contains most of the energy of volcanic tremor at Mt. Etna during upcoming unrest (Figure 2b and Figure S7 in Supporting Information S1).

Gravity Data
The study of temporal gravity changes can provide unique insight into the processes driving volcanic activity, since gravity is the only parameter that is directly influenced by the redistribution of underground mass (Carbone et al., 2017). At Mt. Etna, a mini-array of continuously running superconducting gravimeters (SGs) has been in operation since 2016 (Carbone et al., 2019). These instruments provide better quality data than the more widely used spring gravimeters, since they are not affected by environmental parameters (Van Camp et al., 2017). Indeed, thanks to their high precision and long-term stability, SGs allow to track gravity changes in the order of Figure 3. The spatio-temporal variations of dv/v and the broadband seismic stations, shown by red triangles and used in the seismic ambient noise analysis. The three columns show three different time ranges for phase 1, phase 2, and phase 3. The epicentral locations of volcano-tectonic earthquakes are shown by gray circles for each phase. The first, second, and third rows correspond to 5, 10, and 30 days moving windows, respectively. The contour lines show 1,000 m elevation interval. The light-green line shows the eruptive fissures in phase 3. The interpolation was applied by the triangulation method using the GMT software (Wessel et al., 2019).
1-2 µGal, occurring over time scales between minutes and several months. In this study, we utilize the gravity time series from iGrav#16 at SLN (Figure 1a), jointly with the other available time series (Figure 2c).

Soil CO 2 and SO 2 Flux Data
The geochemical observations from CO 2 and SO 2 flux data play a significant role in understanding the relations between degassing and volcanic/tectonic activities. We used the CO 2 data from ETNAGAS network consisting of 14 stations operated by INGV-Palermo (Figure 1a). To remove the effects of temperature, pressure, and wind variations on the geochemical parameters, specific filters were applied to the data set (Figure 2d; Gurrieri et al., 2021;Liuzzo et al., 2013;Paonita et al., 2021). Furthermore, SO 2 flux data from UV-scanner Flame-DOAS network from nine stations established by INGV-OE (Figure 1a) were used in this study (Figure 2e; Paonita et al., 2021;Salerno et al., 2009).

Ground Deformation
We used the GNSS data to evaluate the ground deformation in the considered period. The GNSS permanent network run by INGV-OE consists of 39 stations mainly equipped with Leica Geosystems instruments (Figure 1a). The GNSS raw data coming from the stations are processed on a daily basis using the GAMIT/GLOBK software (Herring et al., 2010). The processed GNSS final daily solutions are referred to the geodetic reference frame specifically designed for the area (Palano et al., 2010). To characterize the volcano deformation state, we calculated the areal strain in the volcano summit area by using the GNSS daily solutions (Figure 2f).

Volumetric Strain Change
Volumetric strain models can provide insights into plumbing system dynamics (Donaldson et al., 2017) as well as can help model seismic velocity changes. For example, seismic velocities can decrease or increase during magma uprising due to the opening and closing of cracks, respectively (Bennington et al., 2015;Brenguier et al., 2008;Hotovec-Ellis et al., 2015;Obermann et al., 2013). The inflation of a pressure source causes extensional and compressional volumetric strain above and at the sides of the source, correspondingly. Deeper sources host larger extension and vice versa (Donaldson et al., 2017). We calculated the volumetric strain due to the ellipsoidal sources proposed in Mattia et al. (2020) for the April-October and October-December pre-eruption periods, placed at ∼5 km b.s.l. and with volume variations of 12 and 21 10 6 m 3 , respectively. We calculated the volumetric strain out to 10 km from the source by computing the infinitesimal strain tensor from the Yang analytical solutions within the half-space (Yang et al., 1988) (Figure S8 in Supporting Information S1).

Results
The multidisciplinary observation of the period surrounding the December 2018 eruption of Mt. Etna volcano yielded distinct time series, which, taken together, helps unfold the events leading to the eruption (Figure 2). We jointly analyzed all the data sets and divided the period preceding and accompanying the 2018 eruption into three distinct phases as chronologically as follows (  Figure S3 in Supporting Information S1). After the eruption onset, dv/v started to increase recovering the whole variation accumulated in October, settling again on values close to the ones observed prior to the pre-eruptive period. The normalized dv/v histogram plot showing the variation of dv/v in time for all the station pairs can be found in Figure S9 in Supporting Information S1.

of 16
The VTs M L > 3.9 display mostly the strike-slip mechanism in the vicinity of Mt. Etna during the seismic swarm that took place in 2018 (Figure 1). The source parameters (strike, dip, and rake) that we calculated are in good agreement with the previous reports given by different sources (e.g., GCMT, Pondrelli, 2002;Scarfì et al., 2013;Scognamiglio et al., 2006). However, we estimated shallower centroid depths <4 km for most of the earthquakes. The two largest earthquakes, 6 October 2018 M L 4.7 and 26 December 2018 M L 4.8 indicate 8.7 ± 1.1 and 3.8 ± 0.7 km centroid depths, respectively. An example of waveform fits for the former earthquake is illustrated in Figure S10 in Supporting Information S1. The centroid depths that we obtained are compatible with the hypocentral depths as previously reported (Alparone et al., 2020). The MTs results are listed in Table S1 in Supporting Information S1.
As for the volcanic tremor, Figure 2b shows two clear increments preceding and accompanying both the Strombolian activity in August (phase 1) and the eruptive activity in December (phase 3). In particular, volcanic tremor showed a sharp amplitude increase preceding by 2 hr the onset of the activity that occurred on 23 August 2018. As for the latter, the volcanic tremor amplitude started increasing at about 08:30 on 24 December 2018 and then a couple of hours before the eruptive fissure opening. Furthermore, it is worth noting that, since the end of October, we observe a slow and steady increase in the average amplitude of the volcanic tremor until the beginning of December eruption. This could be due to the motion of magmatic masses toward more superficial portions of the volcano's plumbing system with consequent pressurization of the shallow part of the plumbing system as observed during previous unrest periods (Cannata et al., , 2018. Conversely, a decreasing trend occurs after the end of the eruption, probably caused by the magmatic pressure loss. The most striking feature in the gravity data from SLN (Figure 2c) is the increase phase, which occurred during October and early November (phase 2). Gravity increased by slightly more than 10 µGal with an average rate of 0.4 µGal/day. Since no significant elevation change took place at SLN, the free-air effect could be neglected and the resulting gravity change would be only associated with underground mass changes. It is important to note that an intense rainfall took place in mid-October 2018 ( Figure S11 in Supporting Information S1), implying that a part of the above increase in the signal is likely due to changes in groundwater mass (Carbone et al., 2019). Since we do not have precise knowledge of the amplitude of the volcano-related gravity change, nor of the depth of the volcanic mass source, we cannot provide a reliable figure for the amount of mass involved in the redistribution process.
With regard to the variations of CO 2 measured by the EtnaGAS network, the most intense oscillations can be observed between August and October 2018 (phases 1 and 2) with three increase-decrease variation cycles of CO 2 flux, reaching the highest intensity in mid-October (Figure 2d). A further significant trend is the period of decrease, which begins in mid-October and continues steadily until December, reaching minimum values that coincide with the start of eruptive activity. Previous to this work, similar CO 2 cycle variations during eruptive periods have been correlated to inputs of magmatic ascent from an intermediate plumbing system (Gurrieri et al., 2021;Liuzzo et al., 2013).
Throughout August 2018 and December 2018, the bulk daily mean plume SO 2 flux (Figure 2e) released by the summit and the eruptive fracture of the 24 December 2018 eruption was featured by a main degassing stage with values higher of those recorded over non-eruptive quiescent conditions at Mt. Etna (∼5,000 tons/day). Between November 2018 and January 2019 (phase 3), the SO 2 mass emission rate increased up to mean values between 10,000 and 12,000 tons/day. Specifically, over the eruption, the emission rates were highly variable with intraday values up to 17,000 tons/day and strictly connected to the eruptive activity observed with the start of the activity at the NSEC crater. The SO 2 emission rates reverse its trend since the second weeks of January to gradually decrease toward the typical degassing quiescent values of Mt. Etna.
The long lasting inflation period started in 2014-2015 and was shortly interrupted during the five main eruptive episodes that occurred till the December 2018 flank eruption. These episodes had not affected the inflating trend whose source was estimated to be located at a depth of about 7 km b.s.l (Cannavò, 2019). The trend was suddenly interrupted the day of eruption, when the deformation data showed a dike-like intrusion pattern, followed by a deflation period (Figure 2f).
Spatio-temporal seismic velocity changes and the corresponding seismicity were shown for 5, 10, and 30 days smoothed window in Figure 3 for different phases. The dv/v does not show any considerable variation for phase 1 although a Strombolian activity was recorded on 23 August at the summit craters (Figures 3a, 3d, and 3g). The seismic velocities indicate a strong reduction in Phase 2 down to −0.45% mostly in the vicinity of the summit and south-eastern part of the summit. The number of VTs in this period increases sharply ( Figure S1) and the hypocenters become shallower between 3 and 6 km (Figure 3b, 3e, and 3h; Alparone et al., 2020). The largest seismic event in this period, the 6 October 2018 M L 4.7, and its aftershocks are located in the southwest near RF. For Phase 3, the dv/v remains at the same level as at the end of Phase 2 (Figure 3c). A strong seismic swarm took place in the south-eastern flank, mostly shallow part of the crust, in the vicinity of the observed largest decrease of seismic velocity during this phase (Figure 3c, 3f, and 3i; Figure S1 in Supporting Information S1).
According to the strain analysis, a significant increase in both the areal and volumetric strains was observed in October-December (phases 2 and 3), thus agreeing with the observed decrease in dv/v. Moreover, the increase in the volumetric strain spatially coincides with the volume sensitive to the dv/v analysis ( Figure S8 in Supporting Information S1). Indeed, most of the seismic stations showing large dv/v variations (Figure 3c) are located in correspondence with a near-surface extensional area above a 4-km thick bulk where volumetric strain changes of ∼0.5 nanostrain are modeled.

Discussion
The period preceding the 24-27 December eruption can be divided into three distinct phases based on the important changes in key observables (Figure 4). The boundary between the first and the second phases is the 6 October 2018 M L 4.7 earthquake, coinciding with the onset of the decreasing trend of dv/v. The boundary between the second and the third phases is the trend change in the gravity signal, from increasing trends to stable values, on 31 October 2018. Finally, the end of the third phase coincides with the end of the eruption as well as with the apparent increase of the dv/v.
In the first phase, up to 6 October 2018, the relative stability of most observables does not indicate any clear precursors, even though this period includes the 23 August 2018 Strombolian eruption (Figure 4a), except for the CO 2 that shows a long-term precursor since the beginning of 2017 (please see Paonita et al., 2021). The 23 August 2018 eruption was considered weak and its source shallow (Borzi et al., 2020). Its only measurable precursory activity was an increase in volcanic tremor amplitude, which began by 2 hr before the eruption. It is worth mentioning here that the volcano displayed a resting period for 15 months before this phase since after the February-April 2017 eruption (Borzi et al., 2020). This phase was specifically associated with the deep recharge of the reservoir at 7-12 km b.s.l. as evidenced by a progressive increase of CO 2 flux (Figure 2d), of the CO 2 / SO 2 ratio and of 3 He/ 4 He ratios (Gurrieri et al., 2021;Paonita et al., 2021) in agreement with the pressure source estimated by deformation (Cannavo' et al., 2019). The seismicity in the corresponding period is scattered down to 30 km of depth (Figure 1c; Alparone et al., 2020).
During the second phase (6-31 October 2018), one of the most striking features is the seismic velocity decrease, which occurred in early October (Figure 2a). Similar short-term variations in dv/v before volcanic eruptions have been observed in many studies around the world. Brenguier et al. (2008) and Obermann et al. (2013) observed dv/v decreases preceding different eruptive periods at Piton de la Fournaise volcano. They show that seismic velocity reduction is spatially localized close to the magma reservoir and it is likely related to opening of the cracks in the volcano edifice, caused by magma pressurization. Studies at Kīlauea volcano showed that seismic velocity increases and decreases stem from continual deflations and inflations of the summit, respectively (Donaldson et al., 2017). Moreover, Olivier et al. (2019) observed a velocity decrease 10 days prior to the 2018 eruption of the Kīlauea volcano. Recently, 17 days before the 2018 Sierra Negra volcanic eruption, a velocity drop has been revealed by Ruiz et al. (2022). At Mt. Etna, the dv/v reduction during that period of phase 2 would be likely due to sudden pressurization of the deep magma batch (Figure 4b), causing shallow (within 3 km) crustal variations, which are opening of fractures just above a source that is pressurizing before migrating (see also, Donaldson et al., 2017). The strain models for such a source ( Figure S8 in Supporting Information S1) produce an extensional regime directly above and below the pressure source models. It is important to mention here that the dv/v could also be significantly impacted by variations of water storage in the top few hundred meters of the crust (e.g., Lecocq et al., 2017). However, the large rainfall only affected the region almost 10 days after (from mid-October) seismic velocity reduction started ( Figure S11 in Supporting Information S1). Our study supports the idea that reduction of medium seismic velocity is due to crack openings within the extensional regime determined from the volumetric strain analysis as was proposed in several previous studies (Brenguier et al., 2008;Donaldson et al., 2017;Hotovec-Ellis et al., 2015;Obermann et al., 2013). Similar patterns were observed at Mt. Etna as well. De Plaen et al. (2019) suggested that seismic velocity decreases result from intense pressurization of magma within the plumbing system for 2013-2014 paroxysmal eruptions. Patanè et al. (2003) indicated that magma migration within the plumbing system at Mt. Etna led to the dilatation of the volcano edifice preceding the 2002 eruption.
During phase 2, an earthquake of magnitude 4.7, the largest event since 2009, occurred on 6 October 2018 along the RF (Figure 4b). Magma ingression could trigger seismicity in the proximity of ascending magma due to changes in physical properties (e.g., stress state, elastic moduli, density, and seismic velocity) and might have triggered this earthquake. Following that, the seismicity increased in the surrounding areas and became shallower (H ≤6 km b.s.l.; Figure 1 and Figure S1 in Supporting Information S1). Alparone et al. (2020) interpreted this seismicity increase as the first clear sign of unrest, which could be linked to the progressive ascent of magma through the upper crust. Due to the growing rate of fractures, the volcanic edifice is likely subjected to the accumulation of damage, which is supported by abundant shallow seismicity as an indication of magma intrusion in this period. Spatio-temporal distributions of dv/v in the vicinity of the summit area within 10 km well coincide with the seismicity distribution and the larger changes in velocities (Figure 3). The temporal dv/v drop in this phase coincides with the M L 4.7 earthquake. To better understand the causal relationship among magma ascent, M L 4.7 earthquake, and seismic velocity drop, we calculated the dynamic stress change (Hill & Prejean, 2015) associated with the M L 4.7 earthquake and it was equal to 0.15 MPa based on the maximum peak ground velocity observations at the used seismic stations (Figure 1). Such a dynamic stress change value seems to be too low to be responsible for the observed 0.45% (1.0-1.5 Hz) and 0.20% (0.1-1.0 Hz) dv/v drop ( Figure S3) compared to previous studies (e.g., Brenguier et al., 2014;Taira & Brenguier, 2016). Hence, both the 6 October earthquake and the dv/v change are likely to be the effects of magma ascent. However, we cannot exclude that the 6 October earthquake played a role in intensifying the dv/v change.
Concurrently with dv/v drop in phase 2, there is an increase in the gravity signal by ∼10 µGal. The dv/v decrease and the gravity increase are accompanied by an increase in SO 2 and a renewed input of CO 2 . When we consider the changes of all the mentioned signals in this period, phase 2 can be interpreted as follows. The deep increase of pressure, which led to the 6 October 2018 main earthquake, caused the opening of fractures and the ascent of magma simultaneously. This synchronous phenomenon generated the decrease of dv/v but not a decrease in gravity as fracturing processes and magma uprise compensate for each other in terms of mass. Subsequently, the magma continued flowing upward and it gave rise to an increase in mass because of compression of pre-existing gas bubbles, gas dissolution (due to the pressure increase), and possibly filling of pre-existing pore space while the fracturing process ceased. Furthermore, GPS areal strain data ( Figure 2f) show a steady increase that could indicate that the volcano inflating trend did not change even in this period.
In the last phase, 1 November to 26 December, volcanic tremor amplitude, GPS areal strain, and SO 2 flux increased while the gravity signal remained stable. The increase of the first two observables could be due to the pressurization of the upper plumbing system. The growth of SO 2 flux likely points to shallow magma batch degassing. The flattening of gravity signal might be linked to no mass variation and then no important magma migration/shallow intrusion (Figure 4c). In this phase, dv/v showed slight oscillations similar to the ones observed during the second part of phase 2. We interpret this phase as characterized by a stable shallow magma batch whose exsolved gases are increasing the measured SO 2 flux. The increase in the fluid dynamics within the shallow plumbing system can also explain the intensification of the volcanic tremor amplitude, whose source is generally associated with the dynamic interactions of fluids with surrounding rocks along magma pathways (Chouet & Matoza, 2013) as well as the plumbing system pressurization suggested by the GPS areal strain. This pressurization phenomenon eventually led to a dike intrusion accompanied by a seismic swarm and culminated with the development of erupting fissures at the eastern base of the NSEC on 24 December 2018 ( Figure 1). The spatial distribution of dv/v indicates that there is a correlation between the highest drop and where the eruption took place (Figure 3). The seismicity mostly took place at a shallower depth (0-6 km) (Alparone et al., 2020) and evolved into the M L 4.8 earthquake on 26 December along the strike-slip FF (Figure 1b). Indeed, the dike intrusion along the flank also caused the acceleration of the large-scale seaward sliding of the eastern flank, which likely contributed to the occurrence of the M L 4.8 earthquake (Aloisi et al., 2020). This flank sliding occurred at rates significantly larger than past slip events (Mattia et al., 2015) and suddenly depressurized the dike intrusion at the shallow depth. This likely remobilized the magma residing at the shallow depth and contributed to ending the 4-day eruption, which otherwise had strong indicators for an event much larger in terms of erupted volumes. Eventually, dv/v starts recovering and displays an abrupt increase 0.2% after this phase. This subsequent increase of dv/v could be due to the closing of the cracks as the volcano shows the deflation behavior. However, it should be noted that the error of dv/v is high (±0.3%) in this period. This is likely due to the changes that took place in the medium during the eruption, which led to sharp modifications of the CCF compared to the reference CCF. Consequently, we do not rely on the post-eruption dv/v data since the individual station pairs show scattered seismic velocity changes ( Figure S9 in Supporting Information S1) and CCFs become unstable after the eruption.

Conclusions
We implemented a seismic ambient noise interferometry technique for the 2018 Mt. Etna flank eruption by calculating the cross-correlations of the continuous signal recorded by a dense network. We used the moving-window cross-spectral analysis on phase-weighted daily stacked vertical components to successfully identify seismic velocity changes in the pre-eruption period mostly localized around the summit and flank areas. The observed temporal variation of seismic velocities (dv/v) is associated with the pressurization of the plumbing system within the volcano edifice during the period preceding the eruption. We integrated dv/v results with VTs, volcanic tremor, volcanic degassing, gravity, and ground deformation to provide more perspective on how it relates to the dynamics of the magma ascent.
Combining all these observables confirmed the relationship between the magma intrusion and the observed reduction in seismic velocities and provided an additional perspective on the impact of two significant earthquakes on the development of magma transport and how they related to the changes in the local stress field and ground deformation.
Ultimately, using dv/v helped provide a better context to understand the lead-up to an eruption that had the potential of becoming much larger but was aborted after the stress release that followed the shallow earthquake caused the intrusion to stop.

Data Availability Statement
The seismic ambient noise analysis in this study was applied via MSNoise python package (Lecocq et al., 2014).

Acknowledgments
The research has benefited from funding provided by the Italian Presidenza del Consiglio dei Ministri-Dipartimento della Protezione Civile and Istituto Nazionale di Geofisica e Vulcanologia, DPC-INGV Allegato B2 2019-2021-Task 10 "TRUST" (grant number: OB.FU.9996.030) (Scientific Responsibility: F.C.). This paper does not necessarily represent DPC's official opinions and policies. P.B. has also received funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Project 407141557. A.C. thanks CHANCE project, II Edition, Università degli Studi di Catania (principal investigator A. Cannata) and grant PIACERI, 2020-2022 programme (PAROSSISMA project, code 22722132140; principal investigator Marco Viccaro). The study was partially funded by the H2020 NEWTON-g project (Grant Agreement 801221). Open Access Funding provided by Istituto Nazionale di Geofisica e Vulcanologia within the CRUI-CARE Agreement. We are very grateful to Mehdi Nikkhoo for his constructive discussion and comments on this work. We are indebted to the technicians and technologists of the INGV Osservatorio Etneo and Sezione di Palermo for enabling and improving the acquisition of data. We also thank the editor and the reviewers for their valuable comments and suggestions that helped us to improve the article.