Aseismic transient during the 2010–2014 seismic swarm: evidence for longer recurrence of M ≥ 6.5 earthquakes in the Pollino gap (Southern Italy)?

In actively deforming regions, crustal deformation is accommodated by earthquakes and through a variety of transient aseismic phenomena. Here, we study the 2010–2014 Pollino (Southern Italy) swarm sequence (main shock M W 5.1) located within the Pollino seismic gap, by analysing the surface deformation derived from Global Positioning System and Synthetic Aperture Radar data. Inversions of geodetic time series show that a transient slip, with the same mechanism of the main shock, started about 3–4 months before the main shock and lasted almost one year, evolving through time with acceleration phases that correlate with the rate of seismicity. The moment released by the transient slip is equivalent to M W 5.5, significantly larger than the seismic moment release revealing therefore that a significant fraction of the overall deformation is released aseismically. Our findings suggest that crustal deformation in the Pollino gap is accommodated by infrequent “large” earthquakes (M W ≥ 6.5) and by aseismic episodes releasing a significant fraction of the accrued strain. Lower strain rates, relative to the adjacent Southern Apennines, and a mixed seismic/aseismic strain release are in favour of a longer recurrence for large magnitude earthquakes in the Pollino gap.

MIDAS provides uncertainties based on the scaled median of absolute deviations of the residual dispersion. The uncertainties have been shown to be realistic and do not require further scaling 9 .

The high-rate GPS analysis strategy
In addition to the standard GPS processing strategy, we performed also a high-rate analysis of the GPS data at the two closest stations, MMNO and VIGG (Fig. S2). The high-rate GPS analysis strategy was performed by using the GIPSY-OASIS II software (ver. 6.3) released by Jet Propulsion Laboratory (JPL, http://gipsy-oasis.jpl.nasa.gov) and the JPLs final fiducial GPS orbits and high-rate (30 sec) clocks products 6 . The analysis strategy performed in this study was mainly based on three steps: first, using the Global Pressure and Temperature 2 (GPT2) troposphere refractivity function, a static solution with a 5-min sampling rate was performed to estimate the tropospheric dry and wet zenith delay and the horizontal gradients as stochastic random walk parameters; second, the dry troposphere contribution was extracted from the previous step and used, as well as the GPT2 mapping function, as input in the second static solution to improve the estimation of the wet troposphere zenith delay contribution; third, both the dry and wet troposphere contributions were used as input to estimate the station position epoch-by-epoch during the kinematic solution. For all the steps, the precise point positioning method was applied to ionosphere-free carrier phase and pseudorange data 2 and a 2nd-order ionospheric correction was estimated by using the IRI ionospheric model 10,11 . Ocean loading was also computed from the FES2004 tidal model coefficients delivered by the Ocean Tide Loading Provider at Chalmers University 5 (http://holt.oso.chalmers.se/loading) to model station motion. Finally, absolute antenna calibration was used to improve the model of the GPS antennas and the ambiguity resolution was performed by using the wide-lane and phase bias (WLPB) method 6 .

InSAR processing strategy
Interpretation of single interferograms encompassing the main shock generated during the seismic crisis was impaired by the presence of large atmospheric patterns that can be only handled by multi-temporal analysis carried out with Advanced Differential SAR Interferometry (A-DInSAR) approaches. Data have been processed with a two-scale (A-DInSAR) approach whereby a sequence of processing steps were carried out at low resolution (small scale) and high resolution (large scale), respectively. The low resolution processing involved a Small Baseline Subset (SBAS) DInSAR processing algorithm as described in Fornaro et al. 12 , with a spatial resolution of about 50mx50m. This algorithm retrieved an estimation of the atmospheric propagation delay patterns as well as the deformation time series on a small spatial scale on almost the whole frame of 40Kmx40Km. In all the CSK interferograms, topographic fringes were removed using the 1 arcsec resolution Shuttle Radar Topography Mission (SRTM) Digital Elevation Model. Moreover, a specific processing has been applied for the estimation and compensation on slowly temporally variable atmospheric components associated with seasonal stratified atmospheric components. The method, similar to that discussed in Fornaro et al. 13 , applies a filter that removes components that would not be captured by the standard spatially low-pass and temporally high-pass filtering of atmospheric patterns. The signal component retrieved at low resolution were used to calibrate the data for the subsequent higher resolution analysis. The latter was carried out via a tomographic based approach which has been demonstrated 14 to achieve improvements in the detection and estimation of topographic and motion parameters of persistent scatterers with respect to classical Persistent Scatterers Interferometry (PSI) methods.
The higher resolution processing benefited from a filtering method, named "Component extrAction and selEction SAR (CAESAR)" which has been recently proposed and patented. CAESAR achieves a favorable tradeoff between the final resolution and the spatial density of monitored pixels. It exploits a covariance based analysis to extract dominant components and achieve a dramatic increase in the detection of persistent scatterers by slightly sacrificing the spatial resolution of the final image. Further details on the CAESAR method are given in Fornaro et al. 15,16 : the final image has a resolution of 9mx9m, achieved averaging on a 3x3 box the full resolution CSK data with a very high density of monitored pixels, see Supplementary Fig. S4. Before modeling, the DInSAR interferograms were downsampled using a resolution-based resampling algorithm 17,18 . The resampling technique permits us to reduce the number of data points from several millions to a set of about some hundreds of points, with the highest density of data points close to the source of deformations.

Time dependent inversion
Geodetic time series contain several signals, related both to the steady tectonic motions and to the transient motions. Time dependent inversions were performed with TDEFNODE 19 . In the TDEFNODE approach the parameters describing the steady motions are estimated simultaneously with the transients. The long term constant (or steady state) geodetic surface velocities are estimated by finding the best fit slope for each time series individually. As regard, the distribution of transient slip on the fault, in TDEFNODE it can be represented in a number of ways. The approach does allow considering both for uniform slip faults and distributed slip on discrete nodes, as well as simple functions describing the spatial distributions of slip on the fault. The slip rate, s, on the fault during an event is described by: where ( ) • ( ) is the spatial dependence, and x is along-strike position on the fault, w is the down-dip position, t is time and A is the amplitude. The time dependence, S(t), of the slow slip event can be set to an impulse, a Gaussian function, a box-car function or a series of overlapping triangles as is done for earthquake time functions 20 . The surface displacement history is found by integrating S(t) over time and by applying the appropriate Green's functions. In TDEFNODE the inversions were done with a combination of grid search and simulated annealing. The quantity that is minimized is the sum of the reduced chi-square of the weighted misfit residuals plus any penalties (i.e., penalties that are assessed for exceeding parameter bounds). Parameter uncertainties are estimated by a linearization at the best fit parameter values, and for this reason are likely underestimated. Here, we used a total of 87772 observations (63972 observations are the ENU component position of the GPS time series, while 23.800 observations are the LOS displacements) and a time history comprising overlapping triangles. In this case, the free parameters for the time history are T 0 , the origin time, and the triangle amplitudes A i (where i is the number of triangles in the time function and A i is given in mm/yr). The aseismic slow slip transient event was represented by an origin time, 24 triangular time function elements, one slip amplitude and fault dimensions (length and width), fault strike, dip and rake and position. The resulting chi-square of the misfit residuals is 2.04, suggesting that the observation uncertainties are slightly underestimated or that a uniform slip fault model is a too simple model to explain the observed deformation.

Slip distribution inversion
The input data set for the inversion are the horizontal cumulative displacements observed at 12 GPS sites (Table T4) and the cumulative CSK interferogram (Fig. S7g). To test variable slip on the fault plane, the best fit fault was discretized into smaller patches of 1 x 1 km. We extended the fault from the surface to a depth of about 15 km in order to capture the entire area affected by the swarm sequence and the concentric displacement pattern revealed in the DInSAR interferograms, and to allow for down dip and along strike slip. We computed the Green's functions, which relate unit slip on individual fault patches to surface displacements at individual points using rectangular dislocations in an elastic, homogeneous, and isotropic half space 21 . The inversion was carried out using a bounded-values, weighted least squares inversion algorithm 22 , to impose a nonnegative constraint on the estimated slip allowing unbounded values for InSAR orbital parameters (that is, offset and ramp), and using a Laplacian smoothing operator to avoid implausible slip distribution.

Stress transfer
It is well established 23-26 that stress transfer after a slip episode (both coseismic and aseismic slip) may load or discharge adjacent structures. It is therefore important to identify which other possible faults in the Pollino seismic gap area that may have been brought closer to failure by the stress changes following the mixed coseismic/aseismic transient slip episode occurred during the 2010-2014 Pollino swarm sequence. For this reason, we calculate the stress changes due to our best-fit model of the aseismic slip on the northern and southern normal faults in the region (ME, CPST and PF faults) using the Coulomb 3.1 code 26 . We assume a dip of 50° and a rake of -90° for each fault and subdivide them into 2X2 km patches (fault dimension about 20X16 km for each fault, friction coefficient 0.4). Obviously, the largest stress increases are on the patches of the causative fault surrounding the area that is modelled as being affected by the aseismic slip during the swarm sequence (Fig. S8). On the southern PF fault, we find little increased stresses (about 0/0.2 bar) only in the north-western tip of the fault plane, whereas on the northern MF and CPST faults do we find decreased stress (about -4/0.8 bar) on the two south-eastern halves of both the planes, and increased stresses (up to 1.5 bar) on their north-western parts. If we consider in the calculation the small patch of estimated aseismic slip from our variable slip model (located in the northern part of the extended fault plane), we find that stress is further decreased on the central parts of the northern faults. Since this feature might not be well resolved by our data, we suggest that the stress variation on the faults of the Pollino area is likely due only to the main patch of our estimated aseismic slip distribution. Tables   Supplementary Table T1. Coordinates and velocities of continuous GPS station in southern Italy. Coordinates, velocities (in Eurasia reference frame) and errors for each GPS station in Southern Italy processed and used for reference frame realization. Site name, latitude and longitude (in degrees), east and north component of the velocity and relative errors (in mm/yr), time of acquisition (start/end in decimal year), GPS network.