Repeating earthquakes and ground deformation reveal the structure and triggering mechanisms of the Pernicana fault, Mt. Etna

Structure and dynamics of fault systems can be investigated using repeating earthquakes as repeatable seismic sources, alongside ground deformation measurements. Here we utilise a dataset of repeating earthquakes which occurred between 2000 and 2019 along the transtensive Pernicana fault system on the northeast flank of Mount Etna, Italy, to investigate the fault structure, as well as the triggering mechanisms of the seismicity. By grouping the repeating earthquakes into families and integrating the seismic data with GPS measurements of ground deformation, we identify four distinct portions of the fault. Each portion shows a different behaviour in terms of seismicity, repeating earthquakes and ground deformation, which we attribute to structural differences including a segmentation of the fault plane at depth. The recurrence intervals of repeating earthquake families display a low degree of regularity which suggests an episodic triggering mechanism, such as magma intrusion, rather than displacement under a constant stress. Repeating earthquakes display variable behaviour across distinct portions of the Pernicana fault system, Mt. Etna, indicating structural differences and episodic triggering, according to a 20 year earthquake dataset, and GPS ground deformation measurements.

R epeating earthquakes, or repeaters, also called multiplets and earthquake families 1,2 , are earthquakes characterized by the same location and fault geometry, but with different occurrence times 3 . Hence, repeating earthquakes affect the same fault area, exhibit the same slip, and share the same waveforms, provided that the medium interposed between source and seismometer has not changed 3 . Repeating earthquakes can be triggered by repeated ruptures of a fault patch driven to failure by aseismic creep on the surrounding fault plane 4,5 . They can be identified by constraining the source areas and/or analyzing the waveform characteristics 6 . Repeating earthquakes have been identified in different settings, such as plate boundary zones 4,7 , active volcanoes 8,9 , and glaciers 10,11 .
Repeating earthquakes have many seismological applications, such as detecting temporal changes in the Earth's structure due to earthquakes 2,12 and volcanic activity 8,13 by coda wave interferometry, reconstructing the fault geometries by accurate hypocenter determinations [14][15][16][17] , measuring the cumulative slip associated with the repeaters activity 4,18 and characterizing the seismic source through the empirical Green's Function approach 19 .
At Mt. Etna volcano, repeating waveforms have often been observed both for long period/very long period (LP and VLP) events and volcano-tectonic (VT) earthquakes 20,21 . Concerning the VTs, most of them are shallow (focal depth <7 km b.s.l.), located in the eastern flank [22][23][24] , and have a magnitude <2.0 23 . The VT spatial distribution is closely linked to the complex structural framework of the volcano. In fact, Mt. Etna is located on the structural domain of the Gela-Catania Foredeep, at the intersection between the front of the Apennine-Maghrebian Chain and the NNW-SSE Hyblean-Maltese Escarpment fault belt 25,26 (Fig. 1a). The eastern and southern flanks of Mt. Etna are dominated by instability, causing sliding of these volcano sectors to the east and south, respectively [27][28][29] . Numerous hypotheses have been proposed to explain the causes of this instability, such as increases in magma pressure in the plumbing system/magma intrusions and/or gravitational spreading and reorganization 27,[30][31][32] . In addition, it has been suggested that flank sliding can facilitate the passive rise of magma 33 .
To the north, the~E-W-trending 18-km-long transtensive Pernicana Fault System (hereafter referred to as PF), mainly characterized by left lateral-normal motion, confines the instability of the Mt. Etna eastern flank 29 (Fig. 1a). PF extends from the NE Rift, affected by magma intrusion and formation of eruptive fissures (such as during the 2001 and 2002-2003 eruptions 34,35 ), downslope to the sea 29 . PF accommodates the seaward sliding of the Etna's eastern flank with an average velocity of 1-2 cm/y 36 . PF has generated intense seismicity instrumentally recorded at least since the 1980s 37 , and the earthquakes associated with its activity are mainly located in the western and central portions, while they are nearly absent in the eastern portion 38 . Depths of the foci are very shallow, mainly between the surface and about 3 km b.s.l. 38 . From 1999 to 2009, a large number of earthquake multiplets with very long repeatability over time affected the PF 21 . It is worth noting that observations of repeating earthquakes generated by PF date back to 1990 39 .
While ground deformation studies have investigated the dynamics of PF 28,36,40 , earthquake data have never been used to characterize the behavior of PF in terms of cumulative slip, as well as to provide insights into the segmentation of the different PF portions. Moreover, it is not clear whether the fault activity is influenced by magma pressurization in the volcano plumbing system and/or by gravitational forces. Hence, we analyzed the repeating earthquakes taking place on the north-eastern sector of Mt. Etna during 2000-2019 to shed light on these aspects. To reconcile seismic and ground deformation information, the dynamics constrained by repeating earthquake analysis has been compared with the results from the analysis of GPS data.  29,61 ). The white triangles indicate the seismic stations, whose signals were used to identify the repeating earthquake families (EMNR, EZPO, and ECBD). The inset in the bottom left corner of a shows a map of Sicily with some structural features (AMC: Apennine-Maghrebian Chain; HME: Hyblean-Maltese Escarpment; CC: Calabride Chain; HP: Hyblean Plateau 75 ). The red star indicates the location of Voragine (one of the summit craters). b Digital elevation model of Mt. Etna 74 with the VT earthquake locations (colored dots) and the Pernicana Fault (black thick line 29,61 ). The size of the dots depends on the VT earthquake magnitude (see black dots in the lower right corner) and the color on the focal depth (see the colorbar). The dashed rectangle in b indicates the area investigated in this work.

Results
between the first and the last earthquakes in a family, ranges from less than 1 h to the duration of the entire analyzed catalog, that is almost 20 years (Figs. 3a, c and 4a, b). On the basis of the definitions given by Igarashi et al. 41 , most of the PF families (47 out of 61) are "nonburst-type", with a lifetime longer than 3 years.
In addition, we calculated the coefficient of variation (hereafter referred to as COV), defined as the standard deviation divided by the mean and showing the probabilistic behavior of a random variable, to quantify the variability within the repeating earthquake families 4 . The COV computed on magnitude shows low values, ranging between 0.1 and 0.4, suggesting relatively constant magnitude values within each family (Fig. 4c). On the other hand, the COV computed on inter-event times (defined as the time span between successive events) shows higher values of~0.6-1.3 (Fig. 4d), as well as a wider range suggestive of a broad variability in the occurrence behavior of the earthquakes belonging to the families. COV values close to 0 indicate periodicity or quasi-periodicity, COV equal to 1 Poissonian recurrence, and COV greater than 1 temporal clustering (e.g., ref. 42 ). It is worth noting that we did not calculate COV on inter-event times in case of doublets. We ruled out the possibility of "artificially high" COV values for small inter-event times by examining the relationship between COV computed on inter-event times and the mean interevent times values ( Supplementary Fig. 4).
Repeating earthquake analysis allowed identifying four PF portions with different features (Fig. 3), called eastern portion (longitude comprised between 512 and 520 km UTM; hereafter referred to as EP), central portion (longitude 508-512 km UTM; CP), western portion (longitude 503-508 km UTM; WP1) and westernmost end of the fault (longitude 501-503 km UTM; WP2). All the repeating earthquake families are located along CP, WP1, and WP2 (Fig. 3), while EP exhibits no repeating earthquakes at all. This was expected, as the area where the eastern portion of the PF is located appears to be almost unaffected by shallow (focal depth <5 km b.s.l.) seismicity, which characterizes the CP, WP1 and WP2 (Fig. 1b). Regarding the temporal distribution of the PF repeating earthquakes, there are three main periods characterized by intense activity (Supplementary We divided the time series of cumulative slips of each PFrelated family (calculated by Eq. 2) into groups according to the longitude of the family centroid ( Supplementary Fig. 6). Successively, per each 1-km-long longitude range, we selected the family with the highest value of total cumulative slip during the analyzed period ( Fig. 5b and Supplementary Fig. 6). The part of the PF  with the maximum cumulative slip (~130-180 cm at longitude of 509-511 km UTM) is the CP, while the WP2 shows much lower cumulative slip, ranging from~15 to~20 cm (Fig. 5). Moreover, CP is also characterized by the most populous families and the highest magnitude values (Fig. 3c, d). On the other hand, it must be noted that WP1 (at longitude of 503-506 km UTM) shows a much higher number of earthquake families (12-19) compared to CP (6) ( Supplementary Fig. 6).
Changes over longitude of the cumulative strain release, the maximum cumulative slip by repeating earthquakes and the GPS ground deformation data were plotted in Fig. 6 to show the spatial variability of the PF. In particular, we estimated the first parameter by considering all the earthquakes of the catalog (not just the repeating earthquakes) falling within 2-km distance from the PF and grouped into 1-km-long longitude intervals, as well as only the PF repeating earthquakes. We computed the second one as the maximum  cumulative slip calculated based on the repeating earthquake families again grouped into 1-km-long longitude intervals ( Fig. 5b and the thick black line in Supplementary Fig. 6). Concerning the GPS results, the stations Crisimo and NS01, belonging to the GPS periodic network, provided East-West displacement data for the whole analyzed period 2000-2019 (Fig. 2c). In addition, the GPS permanent network supplied information only for the most recent time span (2012-2019) but with higher spatial resolution, that allowed constraining both slip and opening variations along PF, modeled as composed of four linear segments (Fig. 7a, b).
Regarding the time-related families, we found 47 temporal links among the families (see examples in Fig. 8). It is worth noting that, although the maximum temporal difference among the events was fixed to 5 days, the median lag value calculated among the detected pairs of close events in time is 1.3 days. Moreover, the median distance among the centroids of the timerelated families is 1.7 km. Hence, the events belonging to such time-related families tend to occur closely in space and time, suggesting how the mechanism temporally linking the different families is quick and effective especially at very close distances.

Discussion
We found numerous repeating earthquakes at PF, which have allowed us to investigate its dynamics. It is possible to divide the whole length of PF into four portions, characterized by different behaviors in terms of VT seismicity and more specifically of repeating earthquakes, which reflect important structural differences.
Moving from East to West, EP not only lacks any family of repeating earthquakes (Fig. 3), but it is almost completely devoid of any shallow seismicity (Figs. 1b and 6a; highlighted also by previous studies, e.g., refs. 23,38 ) and is characterized by both minimum slip and opening as measured by the GPS permanent network (Fig. 7a, b). Such a PF portion appears as a complex shear system made up of en échelon left-lateral faults 29 , whose location is suggested by surface evidence and damage caused to edifices by creep phenomena 29,34,43 .
Conversely, CP is characterized by intense seismicity (Fig. 1b) and a few families with the highest number of events and magnitude, as well as with the highest corresponding cumulative slip (up to 174 cm at longitude 510-511 km UTM; Figs. 3, 5, 6b). This repeating earthquake-derived slip roughly matches the slip measured by the E-W component of Crisimo GPS station during the entire period analyzed (~200 cm; Fig. 2c). In addition, CP is also characterized by higher slip values as constrained by permanent GPS data during 2012-2019 than EP (Fig. 7a). It is also worth noting that for one segment of such a fault portion (longitude 510-511 km UTM) the cumulative strain release, computed by only repeating earthquakes, is similar to the cumulative strain release by the whole seismicity in the area (Fig. 6a), suggesting that most of the earthquakes taking place in this fault part are multiplets. From the structural point of view, CP shows a sinistral slip with an individual N110°trending steep fault escarpment 29 . All these features suggest that PF is here characterized in depth by a single fault plane with asperities working as creep gauges, allowing to measure the fault slip.
Moving westward (WP1), we find the portion of PF with the highest strain release (Fig. 6a) and the highest number of repeating earthquake families (up to 19 families at longitude of 503-504 km UTM; Supplementary Fig. 6). Each of these families is generally composed of fewer events and smaller magnitude compared to the families in CP (Fig. 3). Hence, unlike the highest slip values constrained by GPS data (Fig. 7a), the cumulative slip computed by repeating earthquakes in WP1 is lower (~30-90 cm) than the one computed for CP (Fig. 5), as well as than the one measured by NS01 GPS station (~220 cm; Fig. 2c). On the surface, WP1 appears as an ENE-WSW trending fault defined by a set of extensional fractures, directly, structurally, and kinematically connected to the NE Rift 29,44,45 . All these features are likely to be suggestive of a segmentation of the fault plane in depth. Thus, each fault segment, exhibiting asperities and then giving rise to repeating earthquake families, probably accommodates only part of the total fault slip. It is worth noting how most of the strongest earthquakes occurring on PF (with magnitude ranging from 3.6 to 4.3) take place along WP1 and CP at longitude ranging from 503 to 510 km UTM, characterized by most of the detected repeating earthquake families (50 out of 61; Supplementary Fig. 7).
Finally, the westernmost end of the fault (WP2) is characterized by a decrease in the VT seismicity (Figs. 1b and 6a) and by a few families (1-5) with a few events (2-3) composing each family, as well as by the lowest cumulative slip (Figs. 3, 5 and Supplementary Fig. 6). WP2 is connected to the NE Rift, and has been affected in the past by magma intrusions and formation of eruptive fissures 34,35 . This is also confirmed by the highest opening values, reached in this PF sector, constrained by data from the permanent GPS network (Fig. 7b). For this reason, this PF sector is probably not capable of effectively generating many repeating earthquakes, as well as more generally VT seismicity. It is also worth noting that the earthquakes taking place in WP2 show low magnitude (mostly lower than 2; Fig. 3d). Hence, we cannot exclude that some events generated by this fault sector have been neglected, especially during the first part of the analyzed time interval, when one-component short-period seismometers were used to monitor Mt. Etna seismicity.
Important differences among the portions of PF also concern the shape of the slip histories, as shown in Fig. 5 and Supplementary Fig. 6. Indeed, while WP1 and WP2 show slips temporally concentrated in particular time spans, CP is characterized by slips spread over almost the entire period analyzed.
Moving from a spatial to temporal analysis, three periods showed intense repeating earthquake activity along PF (Supplementary  34,46 , while the third one (April 2010) is related to the tectonic loading due to the sliding of the eastern flank 47 . Hence, also the repeating earthquakes suggest how different phenomena can affect the PF activity, first of all magma intrusion/pressurization of the volcano plumbing system and gravitational forces 27,[30][31][32] . It is also worth pointing out that 11 out of the 13 strong earthquakes with M > 3.5 took place during these three time periods, suggesting that the link between intense seismicity and repeating earthquakes is not only spatial but also temporal.
In this respect, the recurrence behavior of the events belonging to the PF earthquake families can provide further evidence. Indeed, COV values, ranging from~0.6 to 1.3, indicate variable inter-event times and lack of periodicity, and hence a low degree of regularity (Fig. 4d). Thus, it is possible to infer that such events do not derive from a constant stressing rate acting on the asperities, as obtained in some tectonic areas along plate boundaries, where inter-event time COV values are much lower than 1 42,48 . The observed variability in the inter-event times can be due to the interactions of PF with other faults by static stress transfer or even by dynamic stresses, in case of remote triggering from distant earthquakes [49][50][51] . In addition, the occurrence times of the events belonging to the families are influenced by the role played by PF in accommodating the seaward sliding of the Etna's eastern flank, which as aforementioned is mainly driven by intrusion/pressurization of the volcano plumbing system and gravitational forces. Focusing on the westernmost part of the PF connected to the NE Rift (WP2), intrusions and/or pressurization phenomena of the central plumbing system could directly affect the occurrence times of the events. We also calculated the median value of the inter-event time COV per each portion of the fault, except for EP which does not show any repeating earthquake family. WP2 has only one family with more than 2 events and then we calculated only one COV value equal to 1.41, while WP1 and CP showed 0.93 and 0.75, respectively. Hence, no evident differences among the PF portions were found.
We also highlighted temporal links among families (Fig. 8). In particular, the mechanism linking the different families turned out to be quick and effective especially at very close distances, as testified by median lag and distance, among the detected linked pairs of events, equal to 1.3 days and 1.7 km, respectively. As suggested by Chen et al. 49 who analyzed the repeating earthquake sequences in the Parkfield region (California), such temporal relationships among families can be due to: i) common triggering by a local slow-slip transient spanning the clusters, or ii) shortterm triggering between very close-by events. GPS data show a clear acceleration in the cumulative slip during the 2002-2003 eruption (Fig. 2c), and hence the detected pairs of close events in time belonging to the different families, taking place during this period, are likely to result from the acceleration in the eastern flank sliding triggered by the magma intrusion (hence a phenomenon similar to "i"). On the other hand, other event pairs did not occur during periods with evident acceleration in cumulative slip as measured by GPS data, suggesting that the phenomenon (ii) can be also relevant. We identified four portions of PF, characterized by different features of the repeating earthquake families. The eastern portion of PF (EP), whose location is only suggested by surface evidence and damage caused to edifices by creep phenomena, does not show any family of repeating earthquakes. The central portion (CP) shows few very populous families, whose corresponding cumulative slip roughly matches the slip measured by GPS stations, suggesting how here the repeating earthquakes work as effective creep gauges allowing to measure the fault slip. On the other hand, the western portion of PF (WP1) exhibits numerous families with a low number of events and then a much smaller cumulative slip compared to the slip measured by GPS stations. This likely results from a segmentation of the fault plane in depth; each fault segment, exhibiting asperities and then giving rise to certain repeating earthquake families, accommodates only part of the total fault slip. Finally, the westernmost tip of PF (WP2), connected to the NE Rift and affected by magma intrusions, is largely aseismic and characterized only by a few families with low repeatability.
The events belonging to the PF earthquake families lack periodicity and show a low degree of regularity. This suggests that the earthquakes composing the families do not derive from a constant stressing rate acting on the asperities, but rather from episodic triggering phenomena, linked to the role played by PF in accommodating the seaward sliding of the Etna's eastern flank. Such a sliding can be triggered by both magma intrusion/pressurization of the volcano plumbing system and gravitational forces. Furthermore, temporal links among repeating earthquake families have been identified. The mechanism behind such links turned out to be quick and effective especially at very close distances, and hence it consists of (i) common triggering by the acceleration phenomenon in the eastern flank sliding, and (ii) short-term triggering between very close-by events.
Finally, our study shows how the integration of repeating earthquakes and ground deformation data can help investigate the dynamics of a fault in detail, even in a complex volcanic system as Mt. Etna, dominated by a complicated interplay of eruptive and tectonic phenomena. The same approach could be extended to other faults characterized by repeating earthquake activity, to provide insights into their behavior and its changes over time.

Methods
Data. The dataset comprises 1863 VT earthquakes (hereafter simply referred to as earthquakes) with magnitude from 0.5 to 4.3 (average and median magnitude values are equal to 1.6 and 1.5, respectively), located in the north-eastern sector of Mt. Etna and recorded from 1 January 2000 to 31 May 2019 by the permanent seismic network managed by the Istituto Nazionale di Geofisica e Vulcanologia-Osservatorio Etneo (INGV-OE) ( Fig. 1b and Supplementary Fig. 1 52-54 ). The earthquakes were located using the Hypoellipse algorithm 55 and a 1D crustal velocity model, proposed for Mt. Etna by Hirn et al. 56 and subsequently modified by Patanè et al. 57 . In particular, the top of the crustal velocity model is set at 1600 m a.s.l and the altitude of the seismic stations was taken into account.
The temporal distribution of these earthquakes, as well as their cumulative seismic strain release, are plotted in Fig. 2a, b. The strain release of each earthquake was calculated as the root square of the energy, calculated by using the equation: 58 log E ð Þ ¼ 9:9 þ 1: where E is the energy in erg and M the local magnitude. To perform the cross-correlation analysis, we used three seismic stations located on the north-eastern sector of Mt. Etna: EMNR, EZPO, ECBD (Fig. 1a). During the 20 years analyzed, the quality of the sensors equipping these stations improved over time, from analog, one-component, short period (1 s) to digital, three-component, broadband (40 s) seismometers. We selected these three stations because of their locations, very close to the PF, the long-time recording period, and the good signal to noise ratio.
Repeating earthquake detection. We band-pass filtered the signals between 1 and 20 Hz by a Butterworth 2-pole filter. The high-pass filter at 1 Hz was applied to minimize the possible waveform differences at the lower frequency end due to the distinct sensors (short-period and broad-band) installed during the considered 20 years. The low-pass filter at 20 Hz was implemented to reduce the anthropogenic high frequency noise. Once the signals were filtered, 5-second-long signal windows, starting 0.5 s before the P-wave arrival time, were extracted from the vertical component of the seismic signals recorded by the three stations. Due to the short distance between seismic sources and stations, 5-second-long windows comprise both P-and S-phases. Successively, we computed a cross-correlation matrix for each station, populated by the cross-correlation coefficients obtained by comparing the waveforms of earthquake pairs recorded at the same station.
We applied the following method to extract the repeating earthquake families from each cross-correlation matrix 59 (see Supplementary Fig. 2 for a schematic example of the repeating earthquake detection method). The earthquake with the highest number of correlation values above the threshold (0.9, as also reported below) was identified, and all the events well-correlated with it were stacked to find an average family waveform. This stacked waveform was again cross correlated with the original seismic record. All the earthquakes with a correlation greater than the threshold were grouped into a family and removed from the matrix. The same process was applied on the remainder matrix, until the entire matrix was organized into distinct families. Unlike the "bridging technique" 60 , the algorithm described guarantees that all the earthquakes belonging to a family are "similar" to a single waveform, namely the stacked one. A cross-correlation threshold equal to 0.9 was chosen to extract the families. Finally, we merged the results obtained by the three stations: 21 if the earthquake "a" belongs to the family "1" at the station "STA1" and to the family "2" at the station "STA2", the families "1" and "2" are unified into a single family. Examples of waveforms of the repeating earthquakes belonging to a family (named #3) and recorded by all the three used stations are shown in Supplementary Fig. 3.
Successively, to extract the families related to the PF activity, firstly the centroid location of each family was computed as the point having the average coordinates among the earthquakes belonging to the family. Then, the families whose centroid is located at a maximum horizontal distance of 2 km from the PF line (as indicated by the literature; e.g., ref. 29,61 ) shallower than 5 km b.s.l. were considered as related to the PF activity (Fig. 3). The horizontal distance threshold of 2 km was chosen based on the location horizontal errors of the earthquakes, particularly high (up to ̴3 km) in the first part of the catalog (2000 -2003). In addition, such a threshold should account for possible deviations from a perfect vertical dip of the fault in depth. On the other hand, the 5 km depth threshold is typical for the PF seismicity, as suggested by previous studies 38,62 .
Once the repeating earthquakes related to PF were identified, Supplementary  Fig. 5 Families time-related to each other. Following the work by Chen et al. 49 , we sought evidence of possible temporal "links" among the different PF families to understand whether the earthquake families interact with each other in a way that is observable in the relative occurrence times. In particular, we investigated the temporal relationships among the events belonging to the PF families: (i) a given family was taken into account as the reference family; (ii) the temporal differences between the occurrence times of all the events belonging to the reference family and the occurrence times of all the other events not belonging to the reference family were calculated; (iii) all the families showing at least three events with a time difference shorter than 5 days with the events of the reference family were considered "time-related" to the reference family. These steps were repeated as many times as the number of PF families, each time considering a different family as reference.
Ground deformation data. To study the ground deformation taking place in the PF area, we used GPS stations belonging to both GPS permanent 66 and GPS periodic 67 networks, managed by INGV-OE. The permanent network performs a continuous monitoring of the ground deformation, while the periodic network consists of geodetic benchmarks measured by GPS in survey mode at least yearly, to increase the spatial details needed for imaging the complex ground deformation pattern of Mt. Etna volcano. We used Crisimo and NS01 stations belonging to the NE sector of the GPS periodic network (Fig. 1a). These stations, used to monitor the deformation related to the PF especially during seismic and volcanic crises 68,69 , provided ground deformation data for the whole time interval (2000-2019). In particular, we considered the displacement along the East-West direction, reflecting the slip along the PF (Fig. 2c).
In addition, we used eight GPS stations, spatially distributed on both sides of the PF, belonging to the GPS permanent network to estimate the slip distribution along the fault from ground deformation data (see blue diamond markers in Fig. 1a). We processed the raw GPS observations using the GAMIT/GLOBK software 70 to produce time series of daily station positions. To ensure consistency in the results, we processed the longest available period in which all stations were regularly functioning (2012-2019).
To restrict the analysis to the differential deformation due only to the fault kinematics, thus avoiding the definition of an outer reference system, we considered the relative displacements measured along baselines. The baseline variations were calculated from the position time series for all the pairs of stations whose link crosses the fault. We used cumulated baseline horizontal variations in the considered period to estimate the kinematic parameters along the PF. To this end, the PF was divided into 4 segments and the predicted deformation field at the surface was calculated as the superposition of 4 rectangular dislocation models 71 associated with fault segments. The dislocations are vertical, with top at 50 m and bottom at 3000 m of depth. Since the PF mostly shows strike-slip and opening components, we considered only these components of the source mechanisms of the 4 dislocation models. For fixed model geometries, after computing the Green's functions for unit slip and opening, and for each model/segment, the relationship between predicted deformation and segment kinematics (i.e., strike-slips and openings) becomes linear. Thus, the strike-slip and opening parameters of the 4 segments can be simultaneously estimated from measurements by superimposing their contributions and inverting a linear system of equations. Since the system was overdetermined, we solved it in a least-squares sense. In order to suppress unrealistic scattered solutions and guarantee a consistency of contiguity, a spatial Laplacian constraint was introduced in the inverse problem formulation 72 . The Lagrange multiplier was obtained by the U-curve method 73 .

Data availability
The seismic and GPS data are available from the INGV but restrictions apply to the availability of these data, which were used under license for the current study, and so are not publicly available. Data are however available from the authors upon reasonable request and with permission of the INGV (Salvatore Alparone, salvatore.alparone@ingv. it; Andrea Ursino, andrea.ursino@ingv.it). We provided an xls file as Supplementary Data. 1 with the features of all the detected repeating earthquake families of PF (https:// doi.org/10.5281/zenodo.4727613). This file contains: a first sheet with the main features of all the families; a sheet per each family with information about all the earthquakes composing the family, as well as figures with waveforms and cross-correlation matrices.

Code availability
Codes used in this study are available from the authors upon request. Matlab version 2020a was used to generate all the figures.