## 2 Methodology
### Elastic Dislocation Modelling Using InSAR Observations  
To provide an accurate estimate of the earthquake source parameters based on a single rectangular fault plane and to determine where in-depth slip distributions are located, the following analysis of InSAR data and modelling was performed. After processing the four interferograms of Table 2, we started the modeling of displacement data to get the source parameters. In order to reduce the amount of data to handle in the inversion, particularly needed in this area where the radar coherence is maximized, we sampled the InSAR raster maps adopting a regular grid sampling of 2000 m in the far field and 500 m in the near field. With this sampling, we obtained a subset of 3500 points for the ascending and descending datasets. Then, co-seismic displacement modeling was conducted using a consolidated two-step inversion method: a non-linear inversion to solve for the geometry, location and uniform shear dislocation on a rectangular source, and a linear inversion to solve for the slip distribution over the constrained fault plane [30]. In both steps, the underlying model is the rectangular shear dislocation in a homogeneous and elastic half-space [31]. Further, we adopt a strategy to account for the effects of topography [32], ranging in this area from zero to around 600 m and playing a small but not negligible role. For every event, we jointly inverted the ascending and descending displacement maps. The non-linear inversion exploits the Levenberg-Marquardt minimization algorithm to find the global minimum of the cost-function based on the least-squares of the residual between the observed and predicted data [33]; it consists of a Gauss-Newton method, implemented with multiple restarts to guarantee catching of the global minimum. Details of the implementation can be found in [30].  
Figure 3: Ascending (**a,c**) and descending (**b,d**) pattern of Interferometric fringes and displacements including the May 21, M 5.0 aftershock. The biggest dark star indicates the location of the mainshock and the small dark star shows the location of May aftershock.  
After the definition of the source geometry with a uniform slip and rake, we calculated the slip distribution along the fault plane, by using the algorithm described in [34]. It consists of a least squares optimization based on linear inversion of the matrix:  
\[\left[\begin{array}{c}d_{a}\\ d_{d}\\ 0\end{array}\right]=\left[\begin{array}{c}G\\ k\cdot\
abla^{2}\end{array}\right]\cdot m\]  
where \(d_{a}\) and \(d_{d}\) are the observed ascending and descending data vectors, respectively, \(m\) is the vector of the unknown slip values, \(G\) is the Green's function matrix, \(\
abla^{2}\) is a Laplacian operator used to regularize the slip distribution, weighted by the empiric coefficient \(k\)[35; 36]. To prevent the change of the slip direction (back-slip) along the fault, we used additional constraints of non-negativity of the parameters.  
We remark that all the inversions (non-linear and linear, for the mainshock and the aftershock) are carried out with a joint contribution of ascending and descending InSAR datasets, automatically weighted according to the approach described in [37].Furthermore, for every inversion the parameters of a possible ramp that affects the InSAR data were simultaneously assessed together with the source parameters.