Seismic Source Migration During Strombolian Eruptions Inferred by Very‐Near‐Field Broadband Seismic Network

We analyze seismic waves excited by small Strombolian explosions to understand the source process of volcanic explosions. We deployed five broadband seismometers at only 100–300 m away from the active craters of Stromboli volcano, Italy. Moment tensor inversion of the entire seismic signals in the 0.05–0.2 Hz band locates the source at a depth of 170 and 150–200 m west/southwest of the crater where acoustic waves are excited. Contrary, the sources of seismic waves in the 0.2–0.5 and 0.5–1.0 Hz bands are excited almost at the explosion onset and are located close to the crater. We show for the first time that explosions are preceded of about 10–20 s by a small amplitude seismic phase. Semblance analysis shows that this phase is radiated from a depth of 170 m beneath the western part of the crater area. Our analysis indicates that the source moves about 50 m toward the active crater 10–20 s before the explosion occurs at the surface. At the explosion onset, the source moves back to the same location of the small preceding phase. These lateral migrations of the seismic source are estimated by moment tensor inversion and semblance analysis. We suggest that migration reflects the bending of the shallow feeding system toward northeast. Seismic waves are thus reflecting the history pressure generated by the rising of a gas‐rich pocket in the very shallow portion of a magma mush and by the following restoring force occurring after the explosion.

and James et al. (2008) numerically investigate the process of gas slug ascent in the conduit using the equations of motion of liquid magma and equations of state of gas in the slug. James et al. (2008) further include the motion of magma surrounding the gas slug during the ascent toward the surface. Gas slug ascent velocity and associated pressure changes measured in several laboratory experiments (e.g., Llewellin et al., 2012;Seyfried & Freundt, 2000;Vergniolle & Jaupart, 1990) are well matched with the model presented by James et al. (2008). Hence, the results of the numerical simulations and laboratory experiments have been used to quantitatively interpret the gas overpressure and the associated infrasound in terms of slug-driven Strombolian explosion James et al., 2009). This gas slug dynamics has been considered to be an origin of the source of very-long-period (VLP) seismic waves. At Stromboli volcano, seismic VLP signals are often observed from few seconds before to tens of seconds after each explosion (e.g., Gurioli et al., 2014;Harris & Ripepe, 2007;Ripepe et al., 2001). Moment tensor inversion analyses (Auger et al., 2006;Chouet et al., 2003) and polarization analyses (Giudicepietro et al., 2009;Marchetti & Ripepe, 2005;Martini et al., 2007) indicate that the source of VLP is located at about 500 m a.s.l. below, but outside of the crater rim in the Sciara del Fuoco slope. The estimated moment tensor solutions are explained by the opening/closing of a crack embedded in the volcanic medium (Chouet et al., 2003). Tilt motions recorded at Stromboli volcano revealed inflation-deflation cycles of ground deformation located close to the summit area (Genco & Ripepe, 2010;Ripepe, Delle Donne, et al., 2021). These inflation-deflation are repeatedly observed every a few to tens of minutes. Each explosion is preceded by a slow inflation lasting for about 200 s and the inflation rate increases about 10 times in the last 10-20 s that may indicate a drastic acceleration of the gas dynamics before the explosion (Genco & Ripepe, 2010;Ripepe, Delle Donne, et al., 2021). Ground tilt then shows a sharp deflation for ∼<30 s during and after the explosions.
Laboratory experiment on the acoustic signal associated with the growth, flow and burst of gas bubbles moving in a cylindrical reservoir with a narrow pipe (Ripepe et al., 2001) shows that a low-frequency signal is observed as soon as a gas slug starts rising in the pipe and that a high-frequency signal is excited when the gas slug breaks at the water surface in the pipe. This is similar to what we observe during explosive eruptions at Stromboli. To quantitatively compare the experimental results with the observation, the slug ascent velocity is inferred from the geophysical data. Considering the depth of the VLP source as the centroid of the reservoir, the slug ascent velocity is estimated to be 10-70 m/s (Gurioli et al., 2014;Harris & Ripepe, 2007) from the time delay between VLP and explosive onset visible on the thermal camera. Similarly, at Aso volcano, Japan, the slug ascent velocity is estimated to be 1-160 m/s for Strombolian eruptions (Ishii et al., 2019). However, these observed velocities are much faster than the theoretically predicted velocity of 1.5 and 3.4 m/s for the conduit radius of 1 and 5 m (Batchelor, 1967). Kawaguchi and Nishimura (2015) investigated the spatio-temporal changes of volcanic deformation due to gas slug ascent in an open conduit (James et al., 2008). The result shows that the spatio-temporal characteristics and amplitudes of tilt motion observed at Stromboli (Genco & Ripepe, 2010) cannot be explained by the gas slug ascent in a conduit. These discrepancies suggest the necessity to rethink of the gas slug ascent model in favor of other mechanisms of Strombolian eruptions. More recently, a new model based on the gas flow through a shallow (<300 m) crystal mush (Barth et al., 2019) has been used to explain VLPs at Stromboli as the final stage of a ∼200 s long magma pressurization/depressurization cycle (Ripepe, Delle Donne, et al., 2021).
The objective of this study is to clarify the magma/gas dynamics associated with the origin of the seismic VLP signal at Stromboli. We deploy five seismic stations close to the craters to estimate the seismic sources with the as high resolution as possible. The stations were located at a distance of 100-300 m from the active craters, although steep topographic shape of the volcano, especially the Sciara del Fuoco, made it difficult. We apply two source location methods to clarify spatio-temporal changes of the seismic source locations before, during and after the explosive eruption. We first apply a moment tensor inversion method to the full-waveform and determine the seismic source locations (centroids). We then locate seismic sources with elapsed time based on particle motion 10.1029/2021JB022623 3 of 23 analyses. Finally, we compare these seismological results with a conduit model presented in previous studies and discuss the source process.

Temporary Geophysical Observation at the End of September 2016
The volcanic activity at Stromboli volcano have been continuously monitored by an integrated geophysical network of the University of Florence. The permanent seismic network (Figure 1) consists of four broadband seismic stations (PZZ, ROC, SCI, and STR) equipped with broadband seismometers (Güralp CMG-40T). A permanent five-element infrasonic array (EAR) allows to determine the location of the active crater. This array has an L-shape geometry with an internal spacing of ∼100 m to record coherent infrasonic waves in the 1-10 Hz frequency band . Acoustic data are recorded by a 16 bits acquisition system with a sampling frequency of 54.2 Hz. These seismic and acoustic data are radio-transmitted to the monitoring center of the Department of the Civil Protection (COA) on the island, and these data are collected, processed, and published in real time on a Web site (Valade et al., 2016).
At the end of September 2016, a temporary seismic and acoustic network of five broadband stations (ST1, ST2, ST3, ST4, and ST5) was installed ( Figure 1a). These stations were equipped with Güralp CMG-40T broadband seismometers and the data were recorded by 24 bits Güralp CMG24 digitizers with a sampling rate of 100 Hz. The seismometers were deployed at distances of only 100-300 m from the crater and with the aim to have the better azimuthal coverage. The acoustic signals were also recorded by infrasound microphones at the temporary stations (Lacanna & Ripepe, 2020). These seismo-acoustic stations were installed on 23 September and all the temporary stations except ST4 were withdrawn on 2 October.

Data Characteristics
We analyze seismic data recorded at the eight stations located within 1 km from Central (C) crater ( Figure 1b). We analyze the data relative to about 17 h (from 07:05:00 to 23:59:59 UTC on 26 September 2016) which have high signal-to-noise ratio and no data gap at none of the recording stations. We select 43 seismic events associated with explosions that at ST1 show amplitudes of acoustic waves above 10 Pa. Analysis of the infrasonic array indicates that the acoustic signals are all associated with explosions occurring at the northeast (NE) crater ( Figure 2). An example of seismic and acoustic signals associated with an explosion recorded at ST1 is shown in Figure 3. This station records acoustic signals with the largest amplitudes among all the network. Ground displacement derived from the ground velocity seismogram shows that upward (compression) VLP signal is followed by large downward (rarefaction) amplitude of VLP seismic signal that starts at the onsets of high frequency seismic and acoustic signals (see Ripepe, Delle Donne, et al., 2021 for details). Spectral analysis of the ground velocity seismogram indicates that there are dominant spectral peaks around 0.1 and 10 Hz (Figure 3b). However, peaks are not the same at all the stations. For example, large spectral amplitude around 10 Hz at ST1 may be attributed to site effects, because at other stations (ST3 and PZZ) dominant frequencies are instead at around 0.1 and 2-5 Hz, which are typical frequency ranges at the permanent stations of Stromboli volcano (Ripepe, 1996).
Seismic data are analyzed in three frequency bands: (1) 0.05-0.2 Hz (periods of 5-20 s, VLP1), (2) 0.2-0.5 Hz (2-5 s, VLP2), and (3) 0.5-1.0 Hz (1-2 s, LP). The VLP1 band is the lower frequency part of the VLP (0.05-0.5 Hz, 2-20 s; e.g., Chouet et al., 2003), and represents the dominant frequency content in the seismograms. VLP2 band is the higher frequency part of the VLP, in which small spectral peaks are found at about 0.2 Hz (5 s) and 0.27 Hz (3.7 s) ( Figure 3b). The LP frequency band represents the component between the two typical dominant frequencies: around 0.1 and 2-5 Hz (Ripepe, 1996). Figure 4 shows causal-filtered (fourth order) seismograms in the three VLP1, VLP2 and LP bands. Initial motions with large amplitudes are almost simultaneously detected for the three different frequency bands. Since the seismic network was deployed within a distance of about 1 km from the craters, these phases mainly consist of near-field term of the seismic waves.
Besides, the analysis of the VLP1 band shows that 21 events (out of the 43 events) have a small-amplitude positive signal that precedes by 10-20 s the onset of the explosion (Figure 5a). The amplitude of this phase is only about 2-3 times larger than the ambient noise ( Figure 5b). This preceding phase is more evident when the VLP1 phase is large (Figure 5c), and it is probably masked by the ambient noise or permanent tremor when the amplitude of VLP1 band is small. This preceding phase is visible only in the VLP1 frequency bands and seems to coincide with the drastic acceleration of the ground inflation detected by tilt sensors (see Figure 6 in Ripepe, Delle Donne, et al., 2021, for details). Therefore, the preceding phase is an important signal that may represent the motions of magma and/or gas in the magma reservoir or conduits in the last tens of seconds before the explosive onset.

Moment Tensor Inversion for Different Frequency Band
We first determine source locations (centroids) and mechanisms by using the entire seismic signal in the three VLP1, VLP2, and LP bands. We apply the moment tensor inversion method of Maeda et al. (2011). This method is developed to deal with horizontal seismograms strongly contaminated by tilt motions. The equations used for the method and the computation system are described in detail in Appendix A.
We compute the Green's functions by using the open-source software package of Seismic Wave Propagation Code (OpenSWPC, Maeda et al., 2017). This numerical simulation code of seismic wave propagation is based on the staggered-grid finite difference method with the fourth-order accuracy in space  and second-order accuracy in time (Levander, 1988). The topography of Stromboli is obtained from a digital elevation map with a resolution of 0.5 m (Lacanna & Ripepe, 2020) and is decimated to 10 m for this computation. We assume P-wave velocity p E V of 3.5 km/s, S-wave velocity s E V of 2.0 km/s, and the medium density  E of 2650 kg/m 3 (Chouet et al., , 2003. The node spacing of the point sources is set to be 50 m. Considering the topographic condition at Stromboli, the volcanic edifice is described by 637 point sources. The optimal source location and mechanism is defined by grid searching the minimum misfit between the observed and the synthetic seismograms. The misfit is calculated as (e.g., Ohminato et al., 1998): where U k t n obs ( )  and U k t E M , and zz E M ) with roughly the same polarity.
The estimated source locations for the 43 events in the VLP1, VLP2, and LP frequency bands are represented in Figure 9. The gray color bars indicate the numbers of the optimal source locations estimated for the 43 events. For the VLP1 band, we have that 29 out of the 43 events are located around the rim of the southwest (SW) crater (location A) (Figure 9a). These locations are distributed around 600 m a.s.l. and 150-200 m southwest of the NE   crater. This depth is about 100 m shallower than the previous values derived by moment tensor inversion (Chouet et al., 2003) and ground deformation (Ripepe, Delle Donne, et al., 2021). Misfits 1 E E are estimated to range be-  It is noted that the moment tensor inversion method determines the centroid of the source that explains the entire waveforms for about 30s. VLP1 signal begin shortly before the onset of an explosion ( Figure 3a), but the larger amplitude is recorded only after the explosion occurs ( Figure 4). Hence, a volumetric change at a depth of 170 m beneath the western edge of the crater area mainly occurs during and/or after the explosion (see Ripepe, Delle Donne, et al., 2021 for details). Besides, seismic waves in the VLP2 and LP frequency bands have the onset which mainly coincides with the explosion so that these waves can be associated with the triggering process of the eruption itself.

Source Locations Estimated by Semblance Analysis
We determine source location of the eruption seismic signals by particle motions and semblance analysis. Contrary to the moment tensor inversion, the semblance analysis has the merit to improve the location due to the time resolution. This enables us to examine the spatio-temporal changes before, during and after the explosion, although the assumption that the source is simply isotropic may not be realistic.
As shown in Figure  Particle motions of these three phases at the very-near-field stations (ST1-ST5) show strong rectilinearity (Figures 10b and 10c), While particle motion for the Preceding phase is pointing to the western area of the craters (Figure 10b), for the Main phase1 is directed more to the east toward the craters. In this case, particle motion coincides with the location of the VLP2 and LP frequency bands as defined by the moment tensor inversion (Figures 10b and 10c). Note that particle motion of Main phase2 points again back to the western portion of the crater area ( Figure 10c). These temporal changes suggest that the source of the VLP1 migrates from west to east toward the crater before an explosion and goes back to the west during the explosion itself ( Figure S2 in Supporting Informarion S1).
To better estimate this migration of the seismic source in the VLP1 frequency band, a semblance analysis following the method of Kawakatsu et al. (2000) is applied to the three different phases of the VLP1 frequency band. This method measures waveform coherency among the stations where high coherence indicates that the seismic energy is isotropically radiated from the source to the stations as compressive body waves. Semblance can be defined in order to incorporate also the information contained in the rectilinearity of the particle motions by emphasizing the radial component   E R and subtracting the other components ( , ) E V H as follows: (2) Before calculating Equation 2, each seismogram is normalized to give equal weight to each station. The rootmean-square amplitude of signal of each station is normalized as follows, We calculate the semblance values at each node of the grid. Grid interval is set to be 25 m, a P-wave velocity of 3.5 km/s is assumed (Chouet et al., , 2003 and eight stations are considered to locate 21 events (

Resolution and Accuracy of the Results
Our moment tensor analysis for the three VLP1, VLP2, and LP frequency band shows a significant volumetric component at the source but with clear differences in the distributions of the xx E M , yy E M and zz E M components. The best-fit of the source locations (centroid) for the three VLP1, VLP2, and LP frequency range are located about 100-300 m in the east-west direction. Semblance analysis shows that the source of the VLP1 migrates in the eastward direction before and after an eruption. Such differences in source mechanism and location of the different frequency components (VLP1, VLP2, and LP) of the seismic signal recorded at Stromboli may reflect spatio and temporal changes of magma dynamics during the explosive process. Here, we carefully examine these differences by estimating the errors and performing the inversion using synthetic tests.

Moment Tensors
To verify the reliability in the observed difference of the source mechanism solutions, we calculate model resolution matrices (Menke, 1984;Stein & Wysession, 2003) to understand how well parameters can be resolved by our inversion. The model resolution matrices are calculated from the singular-value decomposition of the Green's functions used in Section 4.  ). Hence, we conclude that the geometry of the source mechanism such as crack or spherical is not well constrained. This means that we are not able to discuss the differences in the source mechanisms between the different seismic VLP1, VLP2, and LP spectral component.

Source Locations Determined by Moment Tensor Inversion
To verify the reliability of our source location, we perform four different tests to estimate possible errors in our location procedure. The first analysis is based on the statistics of waveform residuals and the spatial distribution of the misfit 1 E E (Equation 1). The second one is a jackknife test which allow to examine the effect of station distribution on the source locations. Third and fourth analyses consist of two synthetic tests aimed to investigate the effects of noise and network configuration on the moment tensor solutions. E E distribution covers an area which is  E 50-100 m large in each direction. As shown in Figure 11, the estimated source locations within one standard deviations for the VLP1, VLP2, and LP frequency bands do not significantly overlap to each other. Therefore, we conclude that the VLP1, VLP2, and LP source locations are distributed in the east-west directions beneath the craters and are about 100-300 m apart from each other. Further observation and analysis of more seismic events will probably improve this result and will provide more reliable solutions.
A jackknife test on the moment tensor inversions for the three frequency bands is performed using four different station geometries: (1) only the five temporary stations (ST1-ST5) deployed in the very-near-field, (2) all the stations, including the permanent seismic network, but without ST3, (3) all the stations but without ST4, and (4) all the stations without ST5. We analyze the same event shown in Figures 6-8. The result of the test shows that when the record at one station is not used, the best-fit sources are shifted between 0 and 158 m (three grid points in one direction with one grid point in another direction) from the optimal solution estimated by using all the stations ( Figure S3 in Supporting Informarion S1). These differences in the location are almost the same as the confidence intervals estimated from the probability density distributions of the waveform residuals ( Figure 11). When we use only the five very-near-field stations, VLP1 and VLP2 are located very far (255 and 292 m) from the optimal source. The jackknife test shows that when the number of stations and the azimuthal coverage are insufficient or when we remove the stations with the larger seismic amplitudes, source location becomes unstable.
In addition, we performed synthetic test to evaluate the effect of noise on the location. Synthetic waveforms were calculated for two different source models. The first one is an isotropic expansion/contraction source represented by a moment tensor having only diagonal components ( xx E M , yy E M , and zz E M ) with a ratio of 1:1:1. The second one is a horizontal opening/closing crack similar to the best-fit solution previously found and with a moment tensor having diagonal components of 1:1:3. In both cases, source location is set to be in location A (Figure 9a).
The Green's functions calculated for each station are then convolved with a source time function derived by the vertical component of the VLP1 frequency band recorded at ST1. We then add random noise with amplitude different for each station, according to the signal-to-noise ratio associated with the 43 events considered in our analysis. We thus calculate 43 synthetic pseudo-VLP1 signals by randomly changing the noise. Moment tensor inversion was then applied to each pseudo-VLP1 to calculate the effect of the noise on the source location. This test shows that source location distribution for the two assumed source models ( Figure S4 in Supporting Informarion S1) mainly coincides with the location A at least for 40 out of 43 of the pseudo-VLP1 signals calculated for the isotropic source and 38 out of 43 for the horizontal crack, respectively. In the case of the isotropic source, the misfit ( 1 ) E E of the 40 events located in A is around 0.4%-4.3%, while the other three events have misfit as large as 3.9%-11.1%. When the horizontal crack is considered, the misfit of the 38 events in A is 2.5%-10.7%, while the other 5 events range between 0.9% and 9.5%. We note that residuals of these synthetic tests are smaller than those estimated for the observed data. We suggest that this probably reflects the migration of the source position with time in the recorded VLP1 waveform which is instead assumed stable for the synthetic pseudo-VLP1 calculated for the two source models. Our synthetic tests demonstrate that we are not able to discriminate among location difference spaced out of 1 or 2 grids (50 or 100 m) using only the misfit, which remains almost in the same range of confidence ( Figure 11). We do not exclude that some noise level may introduce a larger error also beyond the one or 2 grids distance. Location D estimated from real data (Figure 9a) may be explained by such noise contamination.
Finally, our second synthetic test examines how the solutions of the moment tensor inversion are affected by the network configuration. We compare our very-near-field network with a network of 18 pseudo stations located between 500 and 1000 m away from the central crater C ( Figure S5 in Supporting Informarion S1) which reproduces the network geometry of previous experiments (e.g., Chouet et al., 2003). We generate pseudo-VLP1 by assuming three source models: (i) isotropic expansion/contraction source (ISO), (ii) opening/closing crack with N45 E striking direction, that is, parallel to the long axis of the crater terrace, and with a ratio between the principal axes of the moment tensor of 1:1:3 (Crack1), and (iii) a horizontal sill-like crack opening/closing source with a moment tensor having diagonal components of 1:1:3 (Crack2). Source location is always set to be in location A (Figure 9a). We perform moment tensor inversion by using pseudo-VLP1 calculated for the two networks. Distribution of 1 E E shows that the area with the small misfits estimated for our very-near-field seismic network is between 30% and 50% smaller than the assumed pseudo network of 18 seismic stations ( Figure S6 in Supporting Informarion S1). Besides, the distribution of 1 E E around the optimal source location A depends on the source mechanism. The misfit distribution expanding toward north-west direction in our result (Figure 6a) may reflect the source mechanism.
We thus compare the source mechanisms expressed as eigenvectors of the moment tensors between the input source models and the output solutions ( Figure S7 in Supporting Informarion S1). Results show that source mechanisms are not so well constrained in the case of our very-near-field seismic network. Comparing the principal axes and the eigenvalues of the solutions, the very-near-field network may not be more appropriate to determine the source mechanisms than distant networks, as has been already discussed by the resolution matrix in Section 5.1. The resolution of the moment tensors for the pseudo network is slightly better, but there is still a trade-off especially for the zz E M component. Hence, we conclude that topography at Stromboli is limiting the geometry of the seismic network and therefore the moment tensor solution are in general not well constrained and may not be used for a detailed discussion on source mechanisms ( Figure S8 in Supporting Informarion S1 and Appendix C).

Semblance Analysis
Reliability of the source location derived by semblance analysis is a priority to understand how realistic the migration of the source inferred for the different phases identified in the VLP1 frequency band before (Preeding phase), during (Main phase1) and after (Main phase2) the explosive onset. Since seismic wavelength is in the order of tens of kilometers, semblance, in theory, does not seem to have the required spatial resolution to constrain the source location. Therefore, we perform semblance analysis using synthetic seismograms by assuming an isotropic expansion/contraction source with two different locations at 50 m away from each other.
The spatial distribution of 3 E S obtained by grid search (Figure 12) indicates how well the source location is constrained. Figure 12a shows the contour lines of 3 E S around the optimal source locations of Preceding phase and Main phase1 for a typical VLP1 event. The contour lines are interpolated at the same intervals in the regions smaller than source grid interval (25 m). The semblance values gradually decrease with distances from the optimal source locations because of the long-period nature of the analyzed seismic waves. As a consequence, semblance close to the optimal source location shows high values and the difference between the source locations for the three phases is not so large. The semblance values for Main phase1 and Main phase2 show also similar characteristics.
We then apply the same analysis to the pseudo VLP1 phase to understand the effects of network configuration, travel times and wave propagation on the semblance analysis. The result of this test shows semblance ( 3 E S ) of 0.928 and 0.932, for both sources located only 50 m apart, respectively. When we overlap the 3 E S contour lines obtained for the two source locations we get a distribution similar to the real VLP. Our test shows (Figure 12c) similar results both when synthetic full-waveforms or in several time windows are used. Contrary to what expected, this strongly suggests that semblance analysis can well determine source location difference as small as 50 m and that the obtained source migration of several tens of meters is plausible. The absolute locations are estimated to be one grid (25 m) west of and above the modeled sources (Figure 12c), and such difference may be caused by the station coverage around the source and by the effect of topography on wave propagation.
Finally, we verify the reliability of the eastward source migration before an explosion (Figure 10a) by measuring the arrival time difference between the Preceding phase ( 0 E t ) and Main phase1 ( 1 E t ) at the ST4 and ST5 seismic stations. Time differences were estimated by cross-correlation using a 4-s-long time window centered around the maximum amplitude recorded in the two phases (Figure 13a). Considering the location of these two stations, the time difference ( ) ( ) t t t t 1 0 ST is expected to be positive when the source migrates eastward. Although the wavelength is so large to make difficult to measure time delays with high accuracy, we found that, as predicted, time delays are positive (Figure 13b).
In conclusion, we believe that the reliability of the source location using six different methods has been largely demonstrated. All the different methods show that different frequency band (VLP1, VLP2, and LP) have different locations and that the seismic source of the VLP1 migrates with time in the east-west direction. Even if absolute location strongly depends on focal mechanisms, the relative locations are very reliable (<50 m) for a given source and thus temporal migration of the source using the semblance method is reliable and significant.

Source Migration
Since the moment tensor inversion strongly reflects the source location and the mechanism of the large seismic amplitudes, the centroid of VLP1 mainly represents the location where seismic waves are excited after the explosion onset. On the other hand, VLP2 and LP frequency band are related to the seismic waves excited almost at the same time of the explosion. This means that the source location of VLP1 coincides with the Main phase2 while source locations of VLP2 and LP are associated with the Main phase1. Eastward migration of the source is revealed by the moment tensor inversion and by the semblance analysis. Relative distances between the position of the seismic source responsible for the VLP1, VLP2, and LP frequency bands are almost the same and very close to the Preceding phase, Main phase1 and Main phase2. Since semblance analysis assumes an isotropic source mechanism, source location derived by the moment tensor inversions which is independent by the source mechanisms has to be considered as more reliable. We correct the source location of the Preceding phase by considering the offsets derived by the two different methodology, and we infer that the seismic source migrates mostly in the east-west 10-20 s before an explosion. The first source (Preceding phase) is located about 150 m west of the craters, just before or at the explosion onset, source is located beneath the craters and after the eruption onset, the source moves back to the west of the craters in the same location of the Preceding phase.
Downward migration of the VLP source was described also by Rowe et al. (1998) as associated with gas bubble bursts in basaltic magma at Erebus volcano. They interpret spectral peaks at periods of 7, 10, and 20 s of VLP seismic signal as a fundamental mode with the two overtones caused by the resonance of magma reservoir/conduit or nonlinear fluid-flow excitation.
Seismic signal at Stromboli also shows spectral peaks at periods of 3.7, 5, and 10 s (Figure 3b). Since these can be interpreted as a fundamental mode (10 s) with two overtones (5 and 3.7 s), the westward source migration following the explosions may be interpreted as the resonance process due to the propagation of pressure waves and/ or fluid flow in the magma reservoir/conduit system. Although the physical parameters to represent the resonance process such as sound velocity and conduit dimensions are not estimated in this study, our observation suggests the existence of a finite-size pressure source.

Amplitude and Temporal Characteristics of VLP Seismic Phases
To understand the origin of the Preceding phase and its link to eruption dynamics, we consider the peak amplitude and its occurrence in time ( 0 E a , 0 E t ) of the 21 seismic signals with the best signal-to-noise ratio and we then calculate the cross correlation with amplitude of the Main phase1 ( 1 E a ), of the Main phase2 ( 2 E a ), of the raw seismic signal ( 3 E a ) and of the acoustic signal ( 4 E a ) recorded at the ST1 station (Table 1 and Figure S9 in Supporting Informarion S1). Besides, amplitudes of this different phases were compared with the time delay ( ) t t 1 0  between the maximum amplitudes of the Main phase1 ( 1 E t ) and of the Preceding phase ( 0 E t ).
We found that the higher correlation (R = 0.8) is between the amplitude of the raw seismic signal ( 3 E a ) and the amplitude of the acoustic signal ( 4 E a ). This correlation supports the evidence (Figure 14e) that amplitude of the seismic signal is highly controlled by the flux of gas and lapilli outside the vent during the explosive mechanism.
In addition, amplitude of the Preceding phase ( 0 E a ) well correlates (R = 0.68) with the amplitude of the Main phase1 ( 1 E a ) (Figure 14c), which seems indicating that the small seismic signal generated 10-20 s before the explosion (Preceding phase), probably by the magma/gas movement inside the conduit (see Ripepe, Delle Donne,  of the Preceding phase and the considered amplitudes (Table 1).

Conceptual Model of the Explosion at Stromboli
Based on these results, we suggest here a possible magma dynamics acting for several tens of seconds before and after the explosive eruptions at Stromboli volcano. Geometry of the assumed feeding system is based on previous models derived to explain magma degassing (Suckale et al., 2016), effusive eruptions (Ripepe et al., 2015(Ripepe et al., , 2017Valade et al., 2016) and ground deformation .
The west-to-east migration of the seismic source before (Preceding phase), during (Main phase1) and after (Main phase2) the explosive eruption is suggesting that the feeding system is bending toward northeast in the last few hundreds of meters before the surface. This result is in agreement with the migration of the VLP seismicity observed during the effusive eruptions (Giudicepietro et al., 2009;Ripepe et al., 2015;Valade et al., 2016) when the back azimuth of the VLPs moves deeper and westward following the lowering of the magma level in the shallow feeding system.
Our conceptual model is considering seismic VLP as part of the slow ground inflation (Genco & Ripepe, 2010;Ripepe, Delle Donne, et al., 2021) that starts about 150-200 s before the explosion onset (Figure 15a). This inflation is explained as generated by the pressure increase due to the accumulation of gas-rich magma below a crystal-rich and dense magma mush (Suckale et al., 2016) acting as a "cap" for the gas (Ripepe, Delle Donne, et al., 2021) and pushing up the last (150-200 m) of magma column.
About 10-20 s before the onset of explosion, the ground inflation rate increases ∼10 times as the gas gets closer to the surface (Genco & Ripepe, 2010). Reaching the surface gas flux accelerates, generating a pressure source that excites small seismic waves (Preceding phase, Figure 15b). About 5 s before the explosion onset (Figure 15c), the pressure source moves eastward to the crater, exciting the Main phase1 in the VLP1 frequency band and the LP seismic waves. The eastward migration of the source is most probably reflecting the geometry of the shallow feeding system which is bending toward northeast in the last 120-170 m below the surface.
When the explosion occurs (Figure 15d), large seismic waves in the VLP2 and LP frequency component of the seismic spectrum are excited below the crater and associated with acoustic waves (e.g., Chouet et al., 1997;Harris & Ripepe, 2007;Ripepe et al., 2001). After the explosion, seismic source (Main phase2) migrates back toward west almost to the same location of the Preceding phase (Figure 15e). This migration is most probably caused by a resonance process which induces the oscillation of the upper part of the magma column while is recovering the equilibrium after the explosive magma/gas released (Ripepe, Delle Donne, et al., 2021). Figure 15. Schematic illustration of explosion process at Stromboli volcano as inferred from the seismic source. (a) About 150-200 s before the explosion onset, ground inflation starts due to the pressure increase caused by the accumulation of gas-rich magma below a crystal-rich and dense magma mush and the gas pushing up the magma column. (b) About 10-20 s before, the acceleration of magma and/or gas bubble motions beneath the crater generates small seismic waves (Preceding phase). (c) From about 5 s before the explosion onset, the pressure source moves eastward to the crater along the geometry of the shallow feeding system, exciting VLP1 and VLP2 seismic waves (Main phase1) and LP seismic waves. (d) Explosion occurs and generates large seismic waves and acoustic waves. (e) The pressure source migrates westward and then moves back to almost the same location of Preceding phase, exciting large and resonant VLP1 seismic waves (Main phase2). See detail in text.

Conclusions
We have examined seismic signals generated by small explosive eruptions at Stromboli volcano recorded at only 100-300 m away from the active craters. Large amplitude of VLP and long-period (LP) signals mainly occur during and after the explosive onset and are thus reflecting by the explosive dynamics. We show for the first time that small amplitude seismic signal in the 0.05-0.2 Hz (VLP1) frequency band is detected 10-20 s prior (Preceding phase) the explosive onset. This is the seismic response of the sharp ground acceleration occurring at the end of the 150-200 s long inflation cycle associated with each explosion at Stromboli volcano (Genco & Ripepe, 2010;Ripepe, Delle Donne, et al., 2021).
Moment tensor inversion reveals that the source of the large VLP1 (0.05-0.2 Hz) amplitude is located around the edge of the southwest (SW) crater at an elevation of 600 m a.s.l. This location does not coincide with the explosive vent, but is 150-200 m southwest of NE crater where explosion are located by acoustic waves. In agreement with previous observation (Ripepe, Delle Donne, et al., 2021), we infer that the large amplitude seismic VLP1 frequency band is excited in the very shallow portion of the magma conduit and is probably caused by a reaction force related to the withdrawal of magma/gas during the explosive eruption. This interpretation finds an evidence in the source mechanisms which is dominated by the vertical dipole component of the moment tensor.
At higher frequency band, seismic sources associated with VLP2 (0.2-0.5 Hz) and LP (0.5-1.0 Hz) spectral component are located almost beneath the NE crater. Their location coincide then with the exploding vent and are excited at the explosion onset. Semblance analysis of the small Preceding phase and of the Main phase reveals that the pressure source moves back and forth from the western side of the crater area to the east below the exploding crater. This migration of the seismic source is in agreement also with the source locations estimated by moment tensor solution, and seems tracking the movement of the gas-rich magma batch during the last tens of second before the explosion. East-to-west movement of the relative location of VLPs has been at Stromboli also observed during lateral effusive eruption and have been explained as following the lowering of the magma level during the drainage of the magma from the very shallow portion of the magma reservoir (Ripepe et al., 2015;Valade et al., 2016). We suggest that migration of VLPs is controlled by the geometry of the shallow feeding conduit located beneath the crater terrace which is bending from the southwest toward the northeast following the main structural feature of the volcano edifice.

Appendix A: Moment Tensor Inversion Including Tilt Motions: Method and Computation System
Observed seismogram generated by a point source moment-tensor can be represented in the frequency domain as where ω is the angular frequency,  The Green's functions are convolved with a cosine function to stabilize the inversion: where , r E t is the rising time of the source time function. To obtain the moment tensor solution independent of the assumed source time function, this cosine function is also convolved into the data vector that consists of the observed seismograms.
Our computational domain consists of a grid with 401 × 401 × 201 nodes in the north-south (NS), east-west (EW) and up-down (UD) directions equispaced by 10 m for a total length of 4 × 4 × 2 km. The center of the domain is set at the summit. Positive E x and E y directions are set to be north and east, respectively, and E z direction is vertical one with downward positive. This Cartesian coordinate is based on the computational one in OpenSWPC. The Perfectly Matched Layer (PML) boundary condition (Chew & Liu, 1996) is adopted to minimize artificial reflections. Note that the output waveform is changed to be upward positive in this software. First, the Green's functions for translational motion G np q , ( ) trans  are yielded at the eight broadband seismic stations by convolving the cosine function in Equation A5 with  0.5 r E t s. We calculate the Green's functions in velocity response and use the observed velocity seismograms for the data vector E d . This is because the Green's functions include static offset in displacements due to near-field term, which may lead to phase reversal in the time function of s E m due to the cyclic convolution. Subsequently, we calculate the Green's functions for tilt motion G np q , ( ) tilt  by computing the vertical displacement at the grids around the stations with 10 m interval, taking the difference of the vertical displacements in the NS and EW directions, and converting to the velocity responses. M , and zx E M ): the resolutions of the diagonal components range from 0.54 to 0.88, and those of the deviatoric components from 0.93 to 0.99. The covariances of the diagonal components range from −0.39 to −0.20 and those of the deviatoric components from −0.02 to 0.02. These results indicate that there are large trade-offs between one diagonal moment tensor and the other two diagonal ones, although trade-offs between the deviatoric tensors and the diagonal ones are small. That is, the diagonal components of the moment tensors are not very well constrained even when using our verynear-field observation data. On the other hand, the well-constrained deviatoric components indicate the dominance of the diagonal components in the source mechanism. The model resolution matrix for the deeper source ( Figure A2c) shows that the resolution of the diagonal components of the moment tensors is estimated to be 0.33, 0.90, and 0.79, respectively. Although the contribution of the smallest eigenvalue to the mechanism solution and the resolution of zz E M for the deeper source are higher than those for Location A, we may not be able to discuss the differences of source mechanism such as crack, cylindrical and spherical pressure source.
We calculate the model resolution matrices using the Green's functions computed at the pseudo seismic network ( Figure S5 in Supporting Informarion S1) and compare them with the results of the very-near-field network. Distribution of eigenvalues obtained by the singular value decomposition of the Green's functions is shown in Figure  S8a in Supporting Informarion S1. For Location A, the smallest eigenvalue is 1.49% of overall contribution. For the source 100 m below Location A, the smallest eigenvalue is 2.47% of overall contribution. These results show that the smallest eigenvalues contribute to the mechanism solutions compared with the case of the very-near-field network. When we place the smallest eigenvalues with zero to consider the pseudoinverse problem, the model resolution matrices show the trade-offs between one diagonal moment tensor and the other two diagonal ones (Figures S8b and S8c in Supporting Informarion S1): the resolution of the diagonal components of the moment tensors is estimated to be 0.76, 0.87, and 0.45, respectively in the case of Location A, and 0.62, 0.80, and 0.64 in the case of the deeper source. Although the contribution of the smallest eigenvalue to the mechanism solution and the resolution of xx E M are improved by the pseudo seismic network, the results suggest that further ingenuity in network configuration will be necessary to improve the inversion solutions for understanding the volumetric changes in detail.