## 3 Materials and Methods
### GPS Data  
Geodetic monitoring on Etna started early in the 1970s with pioneristic trilateration networks; since then, a continuous evolution of methods and technologies has led to the current configuration and complexity designed for several applications and studies at different space and time scales [10; 11; 12]. Now, with its 39 permanent stations, the Etnean GPS network is one of the largest deployed on an active volcano (Etn@net, [13]). For surveillance purposes, most of the stations are acquired every second and processed in real-time at high-rate (1 Hz) by using the epoch-by-epoch algorithm of Geodetics \({}^{\text{\textregistered 0}}\) RTD software package [14]. High-rate GPS time series contain both systematic errors (e.g., multipath) and unmodeled stochastic noise, which can be in the same order of magnitude of the displacement signal. In this study, we applied an ad-hoc strategy in order to minimize the noise level of the time series. To this end, the high-rate data from continuous GPS stations belonging to the INGV Etnean GPS network were processed in kinematic mode (which means assuming that the station is continuously moving) by using three different algorithms: (i) precise epoch-by-epoch solutions [14], by using the Geodetics RTD software which provides independent estimates of positions for each epoch of dual-frequency; (ii) Precise Point Positioning (PPP) solutions [15] by using the GIPSY software; and (iii) differential phase kinematic [16] by using the TRACK software [17], which uses the ionosphere-free LC observable double-differences with respect to nearby fixed stations. All the available raw GPS data were decimated at 30 s and the solutions were sub-sampled with medians in windows of 10 min. As a reference fixed station we chose ENIC (see Figure 3). To highlight the displacement signals in the noisy time series of GPS components, we calculated the baseline variations; this allowed us to significantly reduce common sources of noise (e.g., atmosphere, constellation geometry, etc.) which were predominant in the high-rate GPS data. In order to enhance the useful variations against the background noise and  
Figure 5: Observed tilt vectors cumulated during the first lava fountain indicating an edifice deflation. MMT, MDZ, MSP and EC10 stations were out of order. to further reduce the geodetic correlated noise, we combined the obtained time series from the different processing software by a weighted mean optimized to minimize the resulting associated uncertainties [18]. This can be obtained by considering the inverses of estimated variances as weights (covariance matrix is also provided in the Supplementary Materials). Figure 6 shows the time series calculated as the mean of all the considered baseline relative variations, expressed in parts-per-million (ppm). To estimate the displacements and their associated uncertainties, a piecewise linear regression was carried out on the baseline time series. The choice of using a piecewise line to infer the slopes followed a physical time continuity. The break points for the segmented line were set to the beginning and ends of the paroxysms. A robust regression algorithm [19] was applied, by which the piecewise regression line was made to fit the data set as closely as possible while minimizing the sum of squares of the differences. To verify the stability of the regression results, we tried to linearly fit independently the episodes (e.g., yellow line in Figure 6C). It is worth noting that the variations are quite similar and comparable with the ones obtained by piecewise regression. The only discrepancy with our approach is the period in between the second and third paroxysms when the independent linear fit does not show the same positive step as in the piecewise line. However, this does not affect the coherent deflective behavior during the third paroxysm. To infer the sources, we considered the total variations of the regression lines in the eruptive periods.