On the use of dense seismo-acoustic network to provide timely early warning of volcanic eruptions

. The Stromboli volcano is well known for its persistent explosive activity. On July 3 rd and August 28 th 2019, two paroxysmal explosions occurred, generating an eruptive column that quickly rose up to 5 km above sea level. For the first eruption, the Toulouse Volcanic Ash Advisory Center (VAAC) issued a volcanic ash advisory to the civil aviation users with a two-hour delay. The various processes of these events were monitored in the near and far fields by infrasonic arrays up to distance of 3700 km and by the Italian national seismic network at a range of hundreds of kilometres. Using state-of-the-art propagation modelling, we identify the various seismic and infrasound phases for precise timing of the eruptions. We highlight the advantage of a dense seismo-acoustic network to enhance the monitoring capability of a global network at a regional scale for providing both a reliable source characterisation and a timely early warning to VAACs. localization of both natural and man-made events with great precision. In this study, we identify the various infrasound wavefields of the 2019 eruptions recorded at near and far-field infrasound arrays as well as by seismic stations part across Italy. We then characterize the volcanic source for hazard assessment by using far-field measurements capable of providing timely early warning to VAACs.


Introduction
Located in the Aeolian Islands, in Southern Italy, the Stromboli volcano (38.789N 15.213E, 920 m) is known worldwide for its persistent explosive activity of mild intensity from its open vent summit craters. Its ordinary explosive activity, repeating regularly at a rate of ~10 explosions/hour [1] , emits scoria and ash up to a height of ~100-200 m above the craters, with ejecta fallout typically confined to the crater terrace. This activity has been extensively studied with multiple geophysical observations, spanning from ground deformation and seismicity [2,3,4] , infrasound acoustics [5,6] , infrared thermometry [7,8] and videogrammetry [9,10] . The ordinary activity at Stromboli is occasionally punctuated by major explosions and paroxysms [11,12] . Paroxysms represent the largest-scale explosive events being able to generate convection plumes up to a height of 10 km and eject meter-sized blocks up to distances of 2 km from the vents, thus possibly affecting settled areas of Stromboli island. At least, 25 paroxysms were recorded between 1558 and 1993 [11] , while 4 events have occurred during the last 20 years.
The main risk at Stromboli is due to tsunamis, which are able to affect not only the island but the whole coastal regions of southern Italy [12] . Caused by landslide or sector collapses, as well as pyroclastic flows entering the sea, at least 8 tsunamis related to the Stromboli volcano happened between 1879 and 2003 [11,13] . On December 2002, a subaerial and submarine landslide at Stromboli volcano [14] triggered a tsunami wave that reached locally a height of 10 meters, and impacted the other Aeolian islands and on coastal regions in Southern Italy [13,15] .
Based on the petrology of erupted magma, paroxysms are driven by deep magma batches rising mostly in closed system conditions [16] that eventually fragments in the shallow system producing the observed eruptive columns [17,18] . The recent 2019 eruptions of July 3 rd and August 28 th occurred without any clear precursor. The 2019 episodes consisted in large volcanic explosions, few seconds apart, from different summit craters, releasing an eruptive column that reached ~5 km height and produced a fall-out of lithic blocks, decimeter-sized scorias and ash affecting the summit areas as well as the inhabited settlements of the island. The collapse of these eruptive columns produced pyroclastic flows along the Sciara del Fuoco that entered the sea and trigger tsunamis that locally reached wave heights of 1 m [19] .
Following the July 3 rd explosions, the eruptive plume rose up to 4 km above the summit. Because of the unusual suddenness of the eruption, no eruption notification with a VONA (Volcano Observatory Notice for Aviation) was delivered by Italian observatories. The Toulouse Volcanic Ash Advisory Center (VAAC), which contributes to the International Airways Volcano Watch developed by the International Civil Aviation Organization (ICAO) issued a volcanic ash advisory (VAA) to the ICAO with a two-hour delay. On August 28 th , a VONA was sent by the Istituto Nazionale di Geofisica e Vulcanologia (INGV) about 30 minutes after the first eruptions, alerting for the presence of an ash plume rising up to 4 km above sea level. A VAA was then issued by Toulouse VAAC notifying aviation users about the presence of a drifting ash cloud in the vicinity of the volcano.
The 2019 paroxysmal explosions were recorded locally by the multi-technology monitoring network deployed on the island in the range of 2 km from the crater (acoustic and seismic sensors and video camera operated by the University of Florence), at hundred of kilometres across the permanent Italian seismological network operated by INGV, and in the far-field, up to ~3700 km, by infrasound stations part of the International Monitoring System (IMS) completed by national arrays.
Combining observations from seismo-acoustic arrays allows improving operation monitoring methods to discriminate between natural and anthropogenic phenomena [20] . The new era of massive datasets offers an opportunity to examine the propagation of infrasound wavefield in more detail than previously possible. Analyzing records of dense seismic networks such as the transportable USArray [21] or the European AlpArray [22] provides unprecedented spatial detail of the infrasound ground footprint, allowing the detection and localization of both natural and man-made events with great precision.
In this study, we identify the various infrasound wavefields of the 2019 eruptions recorded at near and far-field infrasound arrays as well as by seismic stations part across Italy. We then characterize the volcanic source for hazard assessment by using far-field measurements capable of providing timely early warning to VAACs.

Results
Infrasound observations. The increased number of operating IMS stations and the establishment of regional infrasonic arrays have evidenced an unprecedented potential of such enhanced network in terms of detection capability, in particular for remote volcano monitoring [23,24] . In addition, national seismoacoustic monitoring systems have been developed in central Europe over several decades filling a gap in the global IMS network. Nowadays, the detection and location capability of such a combined network is significantly enhanced, offering a unique opportunity to investigate improved methods for discriminating between natural and artificial acoustic sources, as well as to better understand seismoacoustic coupling mechanism at the earth's interface [25,26] . Figure 1 presents the considered infrasound network surrounding Stromboli volcano over the Euromediteranean region. At the time of the July 3 rd and August 28 th eruptions, easterly stratospheric wind flow prevailing at 30-60 km altitude favored westward stratospheric propagation. Infrasound records at four IMS stations (IS26, IS42, IS37, IS48) and three national infrasound arrays (AMT, OHP, CEA) are analysed. Shown in Table S1 are the detection results of the July and August eruptions in short to medium distances (AMT: 543 km; IS48: 618 km; OHP: 977 km; IS26: 1124 km; CEA: 1509 km), as well as long distances (IS37: 3379 km; IS42: 3711 km), covering western (IS48, IS42) to northern (IS26, IS37) directions. Table S1 summarizes the onset time, duration, back-azimuth, apparent velocity, frequency range and peak-to-peak amplitude of the detections. All arrivals are stratospheric, with average celerities ranging from 284 to 312 m/s. Only the detection at IS37 in August might be thermospheric with a celerity of 267 m/s, about 20 minutes later than expected for purely stratospheric returns. Back-azimuth deviations from the true direction towards Stromboli volcano are all within 2°, only the case of IS37 in August has a 6° deviation, likely related to higher crosswinds on southward propagation path.
Seismic observations. Infrasound waves from the July 3 rd and August 28 th paroxysmal explosions were also detected by the permanent Italian seismological network up to a range of ~700 km. Figure 2 shows the hodochrons of the waveforms recorded at 198 and 154 seismic stations for the July and August eruptions, respectively. Records exhibiting abnormal noise level likely related to sensor malfunction have been discarded. Time is relative to the origin of the main explosion (14:45:42 and 10:17:16 for the July and August eruptions, respectively) [27] . Aligning the envelops of band-pass filtered data between 0.5-2 Hz using a constant celerity of 300 m/s highlights the signature of stratospherically ducted air-to-ground phases up to a propagation range of ~600 km [21] . Figure S1 presents the seismic records, band-pass filtered between 2 and 8 Hz, at stations located up to a distance of ~170 km from Stromboli. Out of the 38 available stations, impulsive signals are observed at four stations (MILZ, IFIL, MSRU, USI) with arrival times consistent with T-waves [28] . These signals are associated with a seismic precursor which occurred at 14:44:44 UTC, about 1 minute before the first explosion. The observed signals are not compressional P-waves usually generated by explosive events but hydroacoustic waves travelling at a velocity of approximately 1500 m/s, converted to seismic waves on the flanks of an island or on the continental coastal shelf.

Long-range propagation modelling. For the July and August eruptions, 2D parabolic equation (PE)
simulations clearly highlight stratospheric ducting conditions to the west (Figure 1), well in agreement with the summer conditions of the predominant westward stratospheric wind. The first five to ten ground hits of stratospheric returns are visible by concentric rings of reduced attenuation to the west of Stromboli, also reaching arrays AMT, IS48 and OHP. These stations are all within 1000 km and in the cone where stratospheric returns are modeled and observed, resulting in 2D PE attenuation values between 50 and 65 dB. At larger distances, stations CEA, IS26 and IS42 capture stratospheric phases with higher attenuation values between 65 to 95 dB, while IS37 located at 3379 km on the fringe of stratospheric ducting suffers from high attenuation predicted between 100 and 110 dB. The station markers in Figure 1 are color coded with the observed peakto-peak amplitude, converted to attenuation in dB respective to a reference amplitude at 1 km from the source.
This reference is estimated by applying a least square fit between modeled and observed attenuation values (Table S2). Reference source amplitudes ranging from 1 Pa to 10 kPa with 1 Pa steps are applied to define the respective observed attenuations at the different stations. The minimized quadratic difference occurs for a source amplitude of 1510 Pa on July 3 rd and 466 Pa on August 28 th . The converted station attenuations deduced from the latter results are in the range of 50 to 65 dB for AMT, IS48 and OHP, between 65 to 85 dB for CEA, IS26 and IS42 (excluding the low amplitude station observation of IS42 in August) and 90 to 115 dB for IS37. This method of deriving source signal intensity independently from near-source observations yields results consistent with near-field observations [29] . Figure 2 presents the hodochrons of the seismic waveforms recorded at INGV seismic stations up to a distance of 700 km from Stromboli for the July and August eruptions. Superimposed are theoretical arrival times of direct (up to ~80 km, celerity of 340 m/s) and stratospheric phases (first branch ranging from ~200 km to ~400 km and second branch ranging from ~400 km to ~650 km with celerities between 270 and 305 m/s). At almost all stations, the maximum of the signal amplitude coincides with the first and second stratospheric branches with ray turning heights ranging between 35 to 50 km. The summer atmospheric situations of July 3 rd and August 28 th provide comparable propagation conditions. Using a constant celerity of 300 m/s, mean time residuals reduce from 20±4 s to 5±2 s when accounting for station dependent celerity models. Figure 3 compares the recorded and modeled waveforms obtained from normal mode simulations of the wave equation for the July and August eruptions. For each station, range-dependent effective sound speed profiles are considered. The source time function is a ricker wavelet with a dominant frequency of frequency of 1 Hz. Overall, a good agreement between the observed and predicted arrival times and signal duration is observed. Considering the onset times of dominant arrivals (one per station), the mean arrival time residual is 4.8±0.5 s. Beyond 500 km, dominant arrivals switch from the first to the second stratospheric branch for both events. Due to the stronger stratospheric jet in July, the first stratospheric branch extends beyond 500 km.
Location accuracy. Figure 4 compares the time residuals (SSE) for the July and August eruptions considering a constant celerity of 300 m/s and travel times derived from 3D ray tracing simulations. The considered source region extends from 36N to 44N and 10E to 20E with grid resolution of 0.5 km. The travel times are calculated within a time window of 2 minutes centered on the origin time of the paroxysmal eruptions ( Figure  S1) with a time step of 1 s. The minimum of SSE is calculated over the studied region (blue curves), from which the distance from the corresponding source point and Stromboli is calculated (red curves). Using a constant celerity of 300 m/s, the minimum of SSE is reached at 14:45:48 on July 3 rd and 10:17:46 on August 28 th , with corresponding errors of 6 s and 30 s compared with the true origin time. When considering modeled celerity values using atmospheric models from the European Centre for Medium-Range Weather Forecasts (ECMWF-IFS), the difference with the true origin time reduces to ~2 s (14:45:40) in July and ~3 s (10:17:13) in August, and the corresponding location error decreases from 1.6 km to 1.1 km in July, and 3.7 km to 1.3 km in August. Figure 5 presents the geographical distribution of the time residuals for the July and August eruptions. Time residuals are calculated using the origin time given by the minimum of SSE (Figure 4, solid blue curve). Due to the geographical distribution of the stations, the location capability of the seismic network is worse west and east of Stromboli. For comparison, the location results obtained using back-azimuths at seven infrasound arrays only are superimposed. The red ellipses show the 95% confidence location using wind (ECMWF-IFS) corrected back-azimuths (Table S1) and a measurement uncertainty of 1° for arrays of kilometric aperture [30] . When considering observations at infrasound arrays only, the source is located 31 km and 16 km north of Stromboli with median error ellipse area of 4780 and 8705 km 2 for the July and August eruptions, respectively.
Notification of volcanic activity. During the last decade, the deployment of infrasound stations at local, regional and global scales has increased significantly the potential of infrasound measurements to mitigate volcanic hazards [31] . In particular, it was demonstrated that infrasound can play a key role in supporting early warning systems by notifying in near real-time the onset of explosive eruptions, thus reducing the risk of aircrafts encountering volcanic ash in poorly monitored regions. Applying the procedure described by Marchetti et al. [32] , the amplitude of infrasound detections associated with the July and August eruptions are corrected for propagation effect along each source-to-receiver paths in order to infer the source amplitude at a reference distance of 1 km from the crater [33] . The source amplitude is then multiplied by the normalized number of detections per minute, which reflects the time persistency of the signals recorded at the array. The derived infrasound parameter (IP) is used to define the occurrence of significant explosive volcanic activity when a fixed threshold of 60 is exceeded [32] . Figure 6 shows the IP variations at stations AMT, IS48, OHP, and IS26, color-coded according to the alert level. By applying such an approach, notifications could have been delivered up to a distance of 1124 km (IS26) and 618 km (IS48) for the July and August events, respectively. For the July 3 rd eruption, a notification based on infrasound observations at the closest IMS array in Tunisia would have preceded by 105 minutes the VAA which was issued at 17:00 UT. On August 28 th , a VONA was sent by INGV at 10:48 UTC, 31 minutes after the event occurrence, allowing a VAA to be issued by Toulouse VAAC at 10:56 UTC. In this case, the first infrasound based notification using records at IS48, would have been possible at 11:10 UTC only. A notification based on the arrival times of air-to-ground coupled waves recorded by the INGV network would have been delivered ~15 minutes before the VAA, thus highlighting the added-value of regional seismic network to contribute to an operational alert system.

Discussion
The 2019 paroxysmal eruptions of Stromboli were clearly detected by IMS infrasound stations up to a distance of 3711 km. In addition to the IMS network, signals were also recorded by the national seismological network across Italy. The infrasound signature of the eruptions as measured by infrasound and seismic sensors has been modeled at local and regional distances. The amount and variety of observed infrasound phases represents a unique dataset for statistically evaluating atmospheric models, numerical propagation modeling and localization methods.
In a range of ~170 km from Stromboli, transient signals were observed at four seismic stations (MILZ, IFIL, MSRU, USI) with arrival times consistent with T-waves. Such observations suggest the occurrence of an earthquake inside the volcanic edifice and below sea level prior the July 3 rd eruption. The lack of T-wave observations at other stations is explained by masking effect ( Figure S2). For example, along the propagation paths, Panarea and Filicudi Islands may have blocked or diffracted energy such as T-waves cannot be observed at stations IVPL, IVG and IACL (see inset on Figure S2). Furthermore, the presence of shoals between Stromboli and southern Italian coastlines may have reduced the coupling efficiently into P-waves. Modeling the generation and propagation of T-waves is challenging since several important aspects such as the radiation pattern of the source, and the wave coupling at the seafloor interface should be resolved [34] . Therefore, propagation modeling using high-resolution bathymetric data should be realised to confirm this hypothesis and quantify the T-wave attenuation along the propagation paths.
An accurate picture of the ground return footprint of stratospheric arrivals is recovered using 2D PE and ray-tracing simulations, as well as normal mode simulations. When applying the least square estimation of the reference source amplitude, a model-predicted and array-observed infrasound attenuation ranges correspond within values below ±10 dB for July and August (Table S2). The average deviation in dB between model and observation are 3.4 dB for the July and 8.7 dB for the August eruption, respectively. The higher accuracy of the modeling on July 3 rd is explained by a stronger eruptive source producing signals of larger signal-to-noise ratio, thus lowering uncertainty in the amplitude estimation. Another reason is the more stable stratospheric ducting conditions occurring in July, yielding better model prediction with less spreading of the modeled amplitude compared to August. In August, the uncertainty between model prediction and observation is especially high at IS26, as the station is located outside of the modeled stratospheric cone. By minimizing the quadratic difference between the observed and predicted attenuation, a source amplitude of 1510 Pa on July 3 rd and 466 Pa on August 28 th is inferred. This independently derived source amplitude is consistent with near field observations (near 1300 and 650 Pa for the July and August eruptions, respectively) [29] .
Using only IMS infrasound detections, the cross-bearing location shows errors of 31 km and 16 km for the July and August eruptions, respectively. Using the arrival times of air-to-ground coupled waves and realistic celerity models reduce the location errors down to ~1 km with an error in the origin time less than 3 s. This study underlines the benefit of integrating local and regional infrasound networks as well as seismic networks into the global IMS network to significantly improve source location and origin time estimates. Combining the IMS network with dense national seismo-acoustic networks can provide a crucial surplus for improving location and discrimination methods which are used as effective tools for the CTBT verification regime.
These results also highlight the added value of exploiting the synergy between complementary networks to develop at a low-cost efficient monitoring systems for disaster prevention or mitigation. For example, in the field of seismology, rapid earthquake information is essential for both the public, authorities, and policy makers. It has been demonstrated that coupling different information technology detections to seismological arrivals accelerates the detection and the localization of seismic events at a global scale [35] . Pursuing such evaluation would contribute to improve methods for natural hazards monitoring in the perspective of building a timely early warning system to VAACs.

Methods
Data processing. The progressive multi-channel correlation method (PMCC) was applied to the raw pressure sampled at 20 Hz to detect coherent infrasound signals within the background noise and estimate the wave front parameters [36] . The time and frequency window parameters for PMCC processing were set using 1/3-octave band scheme in the 0.1-5.0 Hz band to facilitate detection of broadband signals. Window time durations were set to vary linearly with the period, with a 90% overlap. Such log-scale configuration with variable window length improves discrimination between interfering signals [37] . At a sampling rate of 20 Hz, the resolution of the azimuth is about 1° and 5 m/s for the horizontal trace velocity [30] .
To quantitatively determine the optimum source location using the arrival time of air-to-ground coupled waves observed across the Italian seismic network, a grid-search approach for a single-source model developed by Walker et al. [38] is implemented. The optimum location is estimated by calculating the normalized sum squared errors (SSE) of the misfit between the observed and predicted arrival times at seismic stations.
Propagation modeling. To predict the pressure wave attenuation, infrasound propagation is modeled using a two-dimensional parabolic equation approach (2D PE) up to distances of 4500 km from the Stromboli volcano [39] . Simulations are computed covering the full azimuth range with a step witdth of 0.1°. Atmospheric variability is taken into account using time-varying atmospheric models from the European Centre for Medium-Range Weather Forecasts (ECMWF, http://www.ecmwf.int) (HRES IFS cycle 38r2). ECMWF IFS analysis products contain gravity waves (GWs) originating from different sources such as orography, convection, wind shears, jets and fronts. IFS captures well the synoptic scale meteorology below ~30 km altitude. However, GW activity and GW interaction with the mean flow are not predicted accurately due to imperfect parameterisations [40] , or the lack of fine scale orography leading to underestimated or missing GW activity [41] , for instance. The lack of assimilated observations at stratospheric altitude and a sponge layer starting already in the stratosphere are also responsible for significant temperature and wind biases in the middle atmosphere [42] . In order to account for unresolved wind perturbations, small scale atmospheric variations from the spectral gravity wave model developed by Gardner et al. [43] are added to the ECMWF atmospheric fields. Such empirical corrections significantly affect infrasound propagation in the stratopause region where ray-turning heights are quite sensitive to small-scale atmospheric perturbations [44] . This approach was already successfully employed to explain explosion detections through diffraction and scattering of acoustic energy into the geometric shadow zone in the infrasound range across Europe [45] .
For identitying the propagation paths of the recorded signals, long-range infrasound propagation is simulated using the windy atmospheric sonic propagation ray theory-based method (WASP-3D) which accounts for the longitudinal variation of the atmospheric model along the propagation paths [46] . The ray canonical variables (slowness vector, position, and propagation time) are numerically solved by linearized hydrodynamic equations in spherical coordinates. Following a shooting procedure, the azimuthal deviations of eigenrays sought between Stromboli and infrasound arrays are calculated. Rays are classified and labeled depending on their turning heights and number of ground reflections before reaching the station. Extracted celerity models and azimuthal deviations are median values among the selected rays. This procedure, preferred to costly eigenray techniques, yield reliable station-dependent travel times and azimuthal corrections needed for accurate source localization [26] .
Synthetics waveforms are also calculated at seismic stations where the signal-to-noise of the observed air-to-ground coupled signals allows unambiguous arrival time phase picking. The full-wave propagation method used is a low order range-dependent propagation model obtained from a normal mode waveletbased decomposition of atmospheric perturbations [45] . Selecting the components leading to the most sensitive eigenvalues yields accurate simulations with significantly reduced computational cost.

Data availability
Seismic records were obtained from the national seismic network operated by INGV (http://esm.mi.ingv.it, last accessed August 2020). The vertical component data of broadband seismic stations were downloaded from the INGV International Federation of Digital Seismograph Networks (FDSN) web services (http://terremoti.ingv.it) [48] . ECMWF products, including the atmospheric model analysis are available via https://www.ecmwf.int/en/forecasts/dataset (last accessed August 2020) under CC-BY 4.0 License. Bathymetry and topographic data are from the Shuttle Radar Topography Mission (SRTM30) digital elevation model (DEM) (http://eros.usgs.gov). The DTK-GPMCC software included in the NDC-in-a-Box package designed for infrasound array processing is available upon request. Infrasound records at IMS and national arrays are freely available in the Open Science Framework repository (https://osf.io/ehy3z/). Figure 1 -Attenuation maps derived from range dependent PE simulations at 1 Hz for the July 3 rd (left) and August 28 th (right) eruptions. Color-coded backgrounds represent the surface attenuation predicted by the PE model. The colored station markers represent attenuation of the array-observed amplitudes compared to a reference source amplitude at 1 km of 1510 Pa (July) and 466 Pa (August). For the sake of visibility, the northern most station IS37 is not included in the graphs. Table S2 provides details of the observed and modeled attenuation at all infrasound stations. The stratospheric wind field at 50 km altitude is represented by black arrows (50 m/s reference provided in the lower left corner).   (Figures 4 and 5), performed on the recorded waveforms bandpass filtered between 0.5-2 Hz. To ease the comparison between data and simulation, picks are reported on the modeled waveforms (right). A good agreement beween the predicted and observed arrival times, but also shapes and durations is observed.