## 4 Modeling  
We used high-rate measurements from both GPS baselines and tilt variations to estimate the causal sources of each paroxysmal and, then, the integrated 3D SISTEM ground deformation map to search for a cumulative source feeding the entire sequence. In this way we exploited the availability and potential of HR data showing the continuous and rapid deformation, even with a non-optimal signal-to-noise ratio, and the powerful integrated map, in terms of both spatial resolution and quality of the signal. This allowed us to identify the single-shot sources feeding each paroxysmal and the entire feeding system involved during the whole sequence. Among all the possible HR baselines, we chose the baselines that showed clear patterns, avoiding too-noisy baselines and redundancy (using a single GNSS site for at most two baselines). To calculate the displacements and their associated uncertainties, we fitted the baseline time series with a piecewise linear regression line. The variations were calculated from the segments at the eruptive periods. Similarly, the tilt variations in the same periods were considered in the joint inversion. The analytical model adopted to fit the geodetic data was the finite spherical source model [29] and we used the GAME software [30] to invert the data, applying a hybrid approach of genetic algorithm and pattern search [31]. The objective function was the sum of the squared model residues weighted with their associated uncertainties. The model parameter uncertainties were estimated by using a Jackknife approach [30,32]. The resulting sources show similar depths and coherent horizontal positions (within the area of summit craters) as shown in Table 1. The location of the source was very stable for the three paroxysms; taking into account the differences among the three solutions, the results clearly define a volume located at around 1.8 km depth beneath the summit craters area. The model fitting with the measured deformation data degrades as the energy of the eruption decreases; it is good for the first event, becoming just sufficient for the third event (see Figure 9). In fact, the first two, stronger, paroxysms are those with the strongest signal-to-noise ratio and this allows a better match between predicted and measured deformation. The third one shows a less-fitting modeled deformation, especially on baselines 6-4 and 1-7 (Figure 9), due to a lower signal-to-noise ratio in the deformation data. The last fountain did not show enough deformation to be coherently measured by a sufficient number of considered GPS baselines. Thus, we avoided modeling the source. The model's fit with tilt data is shown in Table 2. \begin{table}
\begin{tabular}{c c c c} \hline \hline  & **Paroxysm 1** & **Paroxysm 2** & **Paroxysm 3** \\ \hline UTM Easting (m) & 498,921 \(\pm\) 192 & 498,694 \(\pm\) 95 & 499,017 \(\pm\) 150 \\ UTM Northing (m) & 4,177,360 \(\pm\) 166 & 4,177,110 \(\pm\) 105 & 4,179,300 \(\pm\) 172 \\ Depth b.s.l. (m) & 1876 \(\pm\) 450 & 1840 \(\pm\) 632 & 1763 \(\pm\) 1032 \\ Volume Change (M m\({}^{3}\)) & \(-1.4\)\(\pm\) 0.2 & \(-1.0\)\(\pm\) 0.2 & \(-0.9\)\(\pm\) 0.2 \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Parameters of the single-shot sources for the first three paroxysms, from HR GNSS and tilt data, together with the associated uncertainties at 2-sigma. Figure 9: Model from the HR data relevant to the first three paroxysms. (**A**) Red circle indicates the position of the source, gray arrows indicate the expected displacements and blue lines indicate the inverted baselines; (**B**) model fitting: green points indicate the measured baseline variations, red points indicate the model predictions. The baselines are sorted in descending order of the error (misfit). In order to model the source feeding the entire sequence, that should represent the average location of the origin of the magma drained during the four paroxysms, we analyzed the longer-term (12 days) but lower-noise data coming from GNSS and DInSAR.