Shifts in the eruptive styles at Stromboli in 2010–2014 revealed by ground-based InSAR data

Ground-Based Interferometric Synthetic Aperture Radar (GBInSAR) is an efficient technique for capturing short, subtle episodes of conduit pressurization in open vent volcanoes like Stromboli (Italy), because it can detect very shallow magma storage, which is difficult to identify using other methods. This technique allows the user to choose the optimal radar location for measuring the most significant deformation signal, provides an exceptional geometrical resolution, and allows for continuous monitoring of the deformation. Here, we present and model ground displacements collected at Stromboli by GBInSAR from January 2010 to August 2014. During this period, the volcano experienced several episodes of intense volcanic activity, culminated in the effusive flank eruption of August 2014. Modelling of the deformation allowed us to estimate a source depth of 482 ± 46 m a.s.l. The cumulative volume change was 4.7 ± 2.6 × 105 m3. The strain energy of the source was evaluated 3–5 times higher than the surface energy needed to open the 6–7 August eruptive fissure. The analysis proposed here can help forecast shifts in the eruptive style and especially the onset of flank eruptions at Stromboli and at similar volcanic systems (e.g. Etna, Piton de La Fournaise, Kilauea).


GBInSAR Data
The NE portion of the summit area of Stromboli has been continuously monitored since January 2003 6 by a GBInSAR system located on a stable section of the flank, about 1.5 km away from the crater terrace ( Fig. 1). GBInSAR Line Of Sight (LOS) is mostly sensitive to the N-S horizontal component of displacement (average azimuth angle = 15°). Negative and positive values of displacement indicate, respectively, a movement toward and away from the sensor. Given the location of the system (Fig. 1), flank instability could be a possible interpretation for movements towards the sensor but not for movements away from the sensor. Since we have observed both movements away and forward the sensor and the modelling indicates that the location and depth of the source are stable, we believe that the changes in the LOS correspond to either inflation (negative values -downslope movement toward the sensor) or deflation (positive values -upslope movement away from the sensor) of the summit area 1-6 .
The GBInSAR system consists of a transmitting and a receiving antenna moving along a rail (3 m long in the configuration deployed at Stromboli) 1 . GBInSAR measures ground displacement along the LOS by computing, via cross correlation, the phase differences between the backscattered signals associated with two consecutive synthetic aperture radar (SAR) images. The ability of GBInSAR to measure . The white triangle identifies the location of the instrument (in the insert, the location of Stromboli volcano is shown). The cumulated displacement maps is produced by the LiSA (Linear SAR) system produced by Ellegi LLC using proprietary GBInSAR technology by LiSALab LLC, a European Commission Joint Research Centre spin-off, and installed at Stromboli by the Dipartimento di Scienze della Terra-Università di Firenze (owner of the system), in the framework of the research agreements (SAR.net, SAR.net2, InGrID and InGrID2015 projects) with the "Presidenza del Consiglio dei Ministri-Dipartimento della Protezione Civile" (Presidency of the Council of Ministers-Department of Civil Protection). Map was generated using ESRI ArcGIS 8.2 platform. volcano deformation depends on the persistence of phase coherence over time. The loss in coherence is primarily due to ground movements, e.g., lava flows or rock avalanches 1,2 . A coherence threshold equal to or above 0.8 is required to recognize deformation areas from a GBInSAR interferogram 2 . Due to the short time (11 min) between two subsequent measurements, interferometric displacements are usually smaller than half wavelength, and phase unwrapping procedures 21 are not necessary. Both the range and cross-range resolutions are on average 2 m × 2 m, with a precision in displacement measurements of less than 1 mm 3 . Displacement rates are computed by differencing the displacements obtained from two consecutive images and dividing by the time spanned. Displacement time series (Fig. 2) are acquired using an algorithm to sum, pixel by pixel, the displacements for every consecutive pair of images and then average that rate over an 8-hour interval 3,22 . Displacement time series of selected points (averaged over 10 pixels) are obtained from cumulative displacement maps with a precision in the displacement measurement of 0.5 mm. Daily displacement can be calculated from the time series.  Fig. 2), corresponding to periods of high-intensity eruptive activity, characterized by frequent and strong explosions and overflows 6 . Time series in the areas characterized by the largest displacement (Fig. 2b) (Fig. 2b,c). The discrete Fourier transform of the GBInSAR time series of the daily displacement of the external rim of the crater terrace (Fig. 2b), shows that the most energetic peaks correspond to a period of 256 days (Supplementary Information). We observed that the GBInSAR displacement time series of areas affected by continuous debris deposition/erosion below the NE crater (NEC) have the same trend of the upper part of the crater terrace (Fig. 2c). Cross-correlation analysis between the two daily-averaged time series reveals that the displacement fluctuation occurred simultaneously (highest absolute-value correlation 0.69 at lag time t = 0). This may mean that the either the cone of the NEC is more affected by the deformation of the conduit, regardless of the erosional/depositional processes or that the swell induced by magma in the conduit is also accompanied by depositional processes on the NEC flank.

Modelling the Deformation
To determine the main parameters of the deformation source, we inverted the GBInSAR displacement employing the software dMODELS 23 . A number of source geometries (spherical source 24 , prolate spheroid 25 , horizontal penny-shaped source 26 and tensile dislocation 27 ), all in a flat, elastic, homogeneous, isotropic half-space are available in dMODELS 23 for several geodetic techniques: leveling, tilt, GPS and InSAR. In the case of InSAR measurements, the software models the changes in range along the radar LOS 28 . Actual volcanic sources are not embedded cavities of simple shape but we assume that these models may reproduce the strain field created by actual storage areas and transport pathways. Given the location of the GBInSARsystem, the angle of incidence of the radar (LOS direction almost perpendicular . The data were acquired using the Leica ADS80 sensor which instrumental vertical and horizontal accuracy is ± 10/20 cm and ± 25 cm, respectively. Map was generated using ESRI 8.2 tm platform. to the volcano flank) and the fact that the system monitors only one flank of the volcano (Fig. 1), we can neglect (in first approximation) the volcano topography and model the deformation as a LOS displacement over a flat half space. The dMODELS software employs a nonlinear inversion algorithm to determine the best-fit parameters for the deformation source by searching the minimum of following the cost function 23 : where, N is the number of data points, P the number of model parameters, d k are the experimental data, m k the modelling results, and σ k the data uncertainties. The non-linear inversion algorithm is a combination of local optimization (interior-point method 29 ) and random search. This approach is more efficient for hyper-parameter optimization than trials on a grid 30 . We use the empirical variogram, a measure of spatial correlation 31 to determine which one the proposed source geometries best fit the deformation (see also Table 1). When two sources would fit the data with a similar precision (e.g, the deformation episode from 08 August 2012 to 08 March 2013 can be explained by either inverse faulting or deflation of a spherical source), we choose the source with the least number of parameters 32 .
Examples of the inversion of InSAR measurements to determine the location of the deformation sources are shown in Fig. 3. We inverted the run-up phases for each period of high-intensity eruptive activity, choosing cumulative maps in the existing dataset (two maps per day; see Table 1). To minimize the influence from sources not related to magma accumulation, we discarded maps affected by atmospheric disturbance. The models reveal a substantial stability of the deformation source over the considered time interval. The best fitting source geometry for the 130 January 2010-August 2014 inflation is a sphere 149 ± 41 m beneath the volcano flank (Table 1; Fig. 4). The source depths can be transformed where H asl is the height of the summit crater terrace above the sea level and L is the distance between the source surface location and the summit crater terrace. Equation (2) allowed us to estimate a source depth d asl = 482 ± 46 m a.s.l. (Fig. 5).

Discussion and Conclusive Remarks
Modelling of the deformation confirmed the presence of a very shallow reservoir consistent with the persistency of magma within Stromboli's conduits, whose existence has been proposed before from the analysis of geochemical data 33,34 , that broadly corresponds with the source location of syn-explosive deformation (350-600 m a.s.l.) 35 . The presence of a very shallow reservoir has been suggested by seismic and deformation data 19 and by petrological studies of lithic ejecta, and in particular from the evidence of pyrometamorphism in tephra accumulated within the crater terrace during persistent activity 36 . At Stromboli, zones of magma accumulation at different depth have been identified. The comparison of the source location recognized by this study with geophysical [8][9][10]19,38 , gas chemistry 39 and petrological 33,34,39,40 data argues that Stromboli's magma plumbing configuration is a multiple-zone storage system (Fig. 6) composed by: i) a deep storage area that feed the most primitive magmas towards the surface (residence time> 55 years) 40 ; ii) an intermediate storage, mainly activated during energetic explosions and flank effusions (residence time = 2-10 years) 40 ; iii) a shallow storage that is involved in all the surface and near-surface phenomena, including explosive activity, central and flank effusions, and non-eruptive  dike injection (residence time = 10-213 days) 33,34 . The comparable period of the GBInSAR displacement cycles and residence time in the shallow storage system 33,34 suggests that the ground displacement in the crater terrace area is controlled by the accumulation of magma in the shallow storage system. The modelling of GBInSAR data allows us to estimate the volume of the shallow magma accumulation and to evaluate its energy budget. The January 2010 -August 2014 unrest was characterized by the accumulation of 2.76 × 10 4 m 3 of magma. This value can be corrected to take into account both the effect of the magma compressibility (β m ) and host rock (β c ) stiffness 41,42 . Compressibility of basaltic magma falls in the range 0.4-2 × 10 −10 Pa −1 ref. 41, while the host rock stiffness has been evaluated following the relationships: where μ is the shear modulus, E the Young's modulus and υ the Poisson's ratio. A range of values for E and υ can be obtained from ref. 43, considering both intact and damaged basalts (β c in the range 0.5-1 × 10 −10 Pa −1 ). The corrected cumulative magma volume during January 2010-August 2014 is 4.7 ± 2.6 × 10 5 m 3 (Table 1), in agreement with previous estimates (7 ± 2 10 5 m 3 ) based on geochemical data 33,34 . These values are also of the same order of magnitude of the volume of magma drained in the starting phase of the 2002-03 and 2007 flank eruptions 11,17,44 , suggesting that the ephemeral vents first depleted the very shallow source. The accumulation rate fluctuated over time (Table 1)  molten and open at a very shallow depth, need heat to its margins at a rate of < 1 MW, much lower than heat release at the crater terrace (423 ± 226 MW) evaluated by ref. 46. The 2014 Stromboli flank eruption was characterized by the development of an eruptive fissure that moved the eruptive centre from the summit crater terrace at 750 m a.s.l. to an ephemeral vent at 650 m a.s.l. The initiation and propagation of a fracture depend primarily on a first order on the potential energy stored in a volcanic edifice when it is loaded 51,52 . Here the loading is primarily related to inflation of the shallow magma storage area. We calculate the strain energy of the shallow intrusions using the strain-nucleus U n where V re is the volume of magma in a single deformation episode, V c is the total magma volume, assumed here as the volume accumulated during the previous time intervals, a is the magma body radius and γ is the melt fraction (on average 50%) 39,53 . The strain energy in each analyzed interval is presented in Table 1. The total strain energy stored during the period January 2010-August 2014 was U n = 3.3 ± 1.8 × 10 14 J ( Table 1). The surface energy needed to open the 6-7 August eruptive fissure, considering a 170 m long (strike dimension, planar distance between the NE crater and the ephemeral vent) and 100 m tall (dip dimension, high difference between the NE crater and the ephemeral vent) fracture and an energy density of about 1.3-43 × 10 6 J m −2 (see ref. 51), is in the range of 10 9 -10 11 J, much lower than the total energy stored at Stromboli during the January 2010-August 2014 inflation. The maximum rate and strain energy accumulated at Stromboli coincide with the inflation preceding the flank eruption of August 2014 (Fig. 7). Ground displacement in the crater terrace area is controlled by the accumulation of magma in the shallow storage system. Fracture opening and propagation is controlled by an increase in the magma accumulation rate that allows building up significant potential energy in a short amount of time. Variations in GBInSAR time series reflect variation in strain energy stored in the shallow source and could be used to forecast the shift between summit and flank eruptions.  Table 1). The maximum rate and strain energy accumulated at Stromboli coincide with the flank eruption of August 2014.