Fault Kinematic Modeling Along a Widely Deformed Plate Boundary in Southern Italy

Convergent plate boundaries are often characterized by widely deformed zones, where coexisting tectonic processes and variable fault kinematics can occur. Here, we quantify this variability along the Africa‐Eurasia deformed boundary in southern Italy, based on the evaluation of geodetic strain rate by recent space geodesy observations and plate motions, which are integrated by main geometric properties of detected faults in the area. We propose a compilation of 160 known faults. We use numerical methods to predict fault kinematics and net slip rate, due to the geodetic deformation field with the inclusion of fault strain accommodation. The obtained tectonic setting is compared with the observable, showing a fault rake agreement of the 73%, which allows us to consider this approach potentially favorable to improve the knowledge of fault kinematics along diffuse plate boundaries, when fault properties are not directly available.

Space geodesy helps to characterize the areas of tectonic compression and extension (Cuffaro et al., 2011;Kreemer et al., 2014;Palano et al., 2015).In those cases, geodetic strain rate computed by Global Navigation Satellite System (GNSS) velocities does not consider the location and geometry of faults in numerical solutions.
In this paper, we compute a new geodetic strain rate field, fault slip rakes and net slip rates in southern Italy, integrating GNSS velocities, global plate motions and available geologic data on fault geometry.We tested this approach by adopting numerical methods with the NeoKinema code (Bird & Liu, 2007), which computes long-term-average anelastic strain rate, derived fault kinematics and fault offset rates.To this scope, we compiled 160 fault lineaments as well as a recent geodetic velocity field.Our main aim is (a) to increase the grid resolution (5 km) of the numerical solutions along the short segment of the AF-EU boundary and (b) to predict tectonic setting in southern Italy with the integration of geodetic and geologic record, comparing results with observable.

Data and Methods
Our approach consists of evaluating regional-to-local geodetic crustal deformation with the inclusion of fault properties within our study area.It is based on the analysis and processing of the available data following three main steps: (a) compilation of the traces of active faults with any available rakes and/or offset rates; (b) collection (i) In southern Italy, we compiled 160 digitalized fault traces located in 12 sectors which are the Southern Tyrrhenian Faults (STF), the Southern Tyrrhenian Transcurrency (STT), the Sicily Onland Transcurrency (SOT), the Ionian Fault System (IFS), the Western Lobe Faults (WLF), the Alfeo-Etna Fault System (AFS), the Calabria Offshore Faults (COF), the Eastern Lobe Compression (ELC), the Western Lobe Compression (WLC), the Calabria-Sicily Faults (CSF), the Southern Tyrrhenian Extension (STE) and the Gioia Basin Faults (GBF), plus two compression fronts (i.e., the Drepano and the Apennines fronts), which characterize the active deformation in the study area (Figures 1-3 and in Supporting Information S1).We merged data sets mainly from Billi et al. (2007), Carminati et al. (2010), Polonia et al. (2016), andDoglioni et al. (2012) with auxiliary information by existing literature (Figures S1-S5 in Supporting Information S1).An additional sector is the Database of Individual Seismogenic Sources (DISS) (DISS Working Group, 2021), but we did not include the sources in southern Tyrrhenian and Ionian Sea to avoid overlaps of crustal deformation results with faults used in Billi et al. (2007) and Polonia et al. (2016).Fault offset rate components are described by the primary and secondary slip, that is, the dominant and the minor sense of the offset for a specific fault.They are not largely available and were constrained as initial condition only for the DISS database (Table S1 in Supporting Information S1).All compiled fault properties are reported in Table S1 in Supporting Information S1.
(ii) A recent GNSS velocity field in the ITRF14 reference frame (Altamimi et al., 2017) was obtained by Billi et al. (2023).We integrated that solution computing additional GNSS rates spanning the 1995.00-2023.00time period, and extending the field westward to Sardinia and northward to the 42°N parallel (Figure 1), to include more continuous GNSS stations.To neglect local velocity outliers (e.g., Hammond et al., 2016;Jiang et al., 2022;Pan et al., 2021), we removed from the final solution all sites biased by large velocity uncertainties and/or showing suspicious movements with respect to nearby sites (∼0.2% of the data set here analyzed).For example, stations close to Mt. Etna were not considered to avoid the short-term volcanic deformation signature (Palano et al., 2022).The velocity field can be considered as an interseismic solution, since no large earthquakes (M > 6) occurred in the study area during the measurement time-window.A final database of 392 GNSS velocities is provided (Figure 1 and Table S2 in Supporting Information S1), using processing methodologies adopted in Billi et al. (2023).
(iii) We modeled data obtained in (i) and (ii) using the NeoKinema code v.5.4 (Bird & Liu, 2007;Shen & Bird, 2022), a kinematic finite element code that estimates long-term average horizontal velocities at the Earth's surface by weighted least squares fitting of geodetic, geological and geophysical data.It solves an inverse problem, based on a spaced gridded interpolation, with the distance weighted approach (e.g., Cardozo & Allmendinger, 2009), where geological heave rates (i.e., the horizontal component of the net slip) and geodetic velocities are merged to estimate regional tectonic deformation field (i.e., strain rate).Stress directions (Heidbach et al., 2018) are useful to constrain the solution (Bird, 2009;Bird & Carafa, 2016;Bird & Liu, 2007;Carafa & Bird, 2016), but, here, we did not utilize them to obtain our deformation field (see in Supporting Information S1 for details).Other adopted constraints are the use of a refined 5 km mesh grid and the AF-EU Euler vector by Altamimi et al. (2017) (Figure 1).Additional methodological details can be found in the supplementary information and in Bird and Liu (2007), Bird (2009), and Carafa and Bird (2016).

Results
The obtained geodetic deformation field with the integration of fault geometries describes the various tectonic settings of southern Italy (Figure 2).We show the strain rate first invariant (2D-dilatation) and second invariant (2nd INV), as defined in Riguzzi et al. (2012) and Kreemer et al. (2014).Three profiles (Figures 2c-2e) are presented to evaluate simulation results from Ionian Sea to South Tyrrhenian Sea, passing through Calabria (section AAʹ, Figure 2d), Messina Strait (section BBʹ, Figure 2e), and Sicily (section CCʹ, Figure 2f).Along those sections, 2D-dilatation and 2nd INV provided by global strain rate model GSRM v.2.1 (Kreemer et al., 2014) are also reported for comparison.Grid solution computed with a 5-km spaced resolution, with no smoothing factor, highlights tectonic deformation as obtained in previous studies (Cuffaro et al., 2011;Palano et al., 2012), for example, compression in the southern Tyrrhenian Sea, extension in central Sicily and Calabria, but it also provides a general compression in the Ionian Sea, due to plate motion constraints.
Higher values with respect to those previously computed by Kreemer et al. (2014) are obtained, sometimes exceeding 300 nstrain/a for both 2D-dilatation and 2nd INV estimations (Figure 2).The deformation field in the Ionian Sea is like the one by GSRM model v.2.1, but higher and oscillating values are observed onshore Sicily and Calabria (Figure 2).
Predicted fault rakes and net slip rates are obtained and a comparison between observed (OFK) and predicted (PFK) fault kinematics is also provided (Figure 3 and Table S1 in Supporting Information S1).Observed fault kinematics are inferred from literature and the kinematics is assigned to faults according to common rake values for normal (N = 270 ± 45), thrust (T = 90 ± 45), right-lateral (R = 180 ± 45), and left-lateral (L = 0 ± 45) faults (see in Supporting Information S1).Input fault properties are here also used as benchmark to validate the proposed approach and model results, and quantitatively to evaluate model goodness, which is the ratio between the number of PFK concordant with the OFK and the total number of OFK.The computed primary and secondary offset rates for each fault (Table S1 in Supporting Information S1) allowed us to derive net slip (Figure 3).
Overall, 117 predicted fault rakes align with the observed ones over the entire database (Tables S1-S13 in Supporting Information S1), with an obtained ratio of 117/160.The fit is best in the DISS, ELC and WLC sectors, worst in the CSF, COF and SOT ones.Obtained net slip rates are comparable with ones used in input for the DISS sector and, overall, are mostly confined under 0.50 mm/a, except for some other values over 0.80 mm/a, such as 1.37 mm/a for the fault n. 3 (e.g., the Aspromonte-Peloritani seismogenic source, Figures 3a and S2a in Supporting Information S1), or 0.83 and 3.05 mm/a for faults n. 37 and n. 149, respectively (Figures S2b and S5n in Supporting Information S1).Moreover, net slip rate is also evaluated along the two compression fronts of the Drepano and the Apennines systems (Table S14 in Supporting Information S1), providing low values of 0.03 and 0.05 mm/a, respectively.

Discussion and Conclusions
This geodetic/geological integrated approach constitutes a procedure which allowed us to compute the interseismic deformation field in and around southern Italy and obtain the long-term tectonic setting.
We suggest that the higher strain rate values in our model relative to those previously obtained (Billi et al., 2023;Devoti et al., 2011;Kreemer et al., 2014;Palano et al., 2012) are due to the dense network of GNSS stations, the dimension of the computation mesh (from 5 km) and the lack of smoothing function application when the solution grid images are generated (Figure 2).This is observed for closer GNSS stations with variable velocities.For example, it can be verified at the Aeolian Islands or in the region south of the Messina Strait (Figure 2), where the huge number of stations provides a 2D-dilatation and/or second invariant over 300 nstrain/a and a compression south of the Messina Strait, respectively (Figure 2 and Data set S6 in Supporting Information S1).In the former region, the inferred strain is due to the superposition of regional tectonic strain and local strains related to the long-term contraction of a magmatic body at Aeolian Islands (e.g., Cintorrino et al., 2019).In the latter region, compression is also powered by the horizontal component of the imposed initial offset rate for the Aspromonte-Peloritani seismogenic source (fault n. 3, Table S1 in Supporting Information S1).
Since our main purpose is to predict long-term fault kinematics by the integration of space geodesy and geological record, discrepancies with respect to previous strain rate solutions are not discussed here in details.Although strain rate field is not here a primary result, but it contributes to evaluate fault tectonics, it is worth noticing that higher values for geodetic strain rate solutions were previously computed in the Aeolian Island region (Bortoluzzi et al., 2010) and recently proposed by Serpelloni et al. (2022) in northwestern Sicily.A similar strain rate order of magnitude is obtained also using our GNSS solution on a 5 km meshgrid and the SSPX software by Cardozo and Allmendinger (2009) (Figure S10 in Supporting Information S1).Hence, we suggest that the regional signature of our GNSS solution is robust and reliable, providing a detailed snapshot of the ongoing crustal deformation over the study area.The use of the NeoKinema code also contributes to predict fault kinematics offshore in the South Tyrrhenian Sea, and in the Ionian Sea, because imposed plate motions are used as velocity constraints, when GNSS rates are not available.The estimated uncertainties of primary and secondary slip rates (i.e., the components of the offset rate, Table S1 in Supporting Information S1) are in order of the 10%, and the overall fault kinematic solution describes the behavior of fault displacements in southern Italy (Figure 4), denoting mostly compression in the Ionian Sea, a strike slip displacement with a normal component in the Sicily Channel and a complex kinematics in southern Tyrrhenian Sea, Sicily and Calabria.
Model goodness in the 13 sectors is measured considering the percentage of alignment between predicted and observed fault kinematics (Figure 4).We find an overall agreement with the existing literature in three sectors at 100% (DISS, ELC and WLC), in one over the 80% (STT), and in two over the 70% (AFS and STE).Lower performances are provided in two other sectors over the 60% (IFS and GBF), in two around 50% (STF, and WLF) and in three less than 40% (SOT, COF, and CSF).Overall, the global comparison over the 13 sectors gives an obtained model goodness of 73% (Figure 4).
We suggest that lower performances from onshore sectors can be attributed to the dense number of GNSS stations and by the refined computation meshgrid.Indeed, the SOT and CSF sectors are within the GNSS network (Figure 4), so that we guess that the analyzed faults locally respond to the near deformation field provided by variable GNSS displacements in Sicily and Calabria, lowering the model goodness with respect to the observed fault kinematics derived for instance by Billi et al. (2010), Polonia et al. (2016) and Barreca et al. (2018).The same reasons, but with different percentages, can be proposed for the STF, GBF and STE sectors (Figure 4), where GNSS velocities at Ustica and Aeolian Islands and variable fault directions strongly constrain predicted fault kinematics (Figure 1), or for the case of Capo Peloro Fault (Fault n. 160. Table S13 in Supporting Information S1), where discrepancy between predicted and observed fault kinematics can be attributed to the high number of different GNSS velocity directions in Sicily and Calabria.Predicted offshore fault kinematics (Figure 4) are mostly constrained by fault location and direction and by far-field kinematic constraints (i.e., GNSS) and mostly by imposed plate motions.This could be the case of the STT, WLF and AFS sectors, where misfit between observed and predicted data can be explained with fault directions with respect to the differences in GNSS velocities between Sicily and Aeolian Islands or Ustica, or relative to the general plate motion directions.Same conclusions can be derived for the IFS sector, where a general right-lateral movement is observed (Polonia et al., 2016), also if in its westernmost part a left-lateral displacement is detected (Gutscher et al., 2017).We suggest that the obtained left-lateral kinematics in the western IFS sector (Figure 4) can be derived by the combination of plate velocity directions and the fault orientation in that domain.
In the COF sector (Figure 4), instead, the extensional kinematics derived by Polonia et al. (2016) and Capozzi et al. (2012) is not aligned with the predicted compression that we obtained, because we suppose that the major constraint due to the imposed convergent AF-EU plate motion at domain boundaries poorly combines with the fault locations and directions in that domain, resulting in overall compression (Figure 4).On the contrary, the distribution of thrusts in the ELC and WLC sectors fully align with the directions of the AF-EU plate motion, resulting in the 100% of model goodness in both the sectors.Finally, the DISS sector (Figure 4) gives a full model goodness also because slip rates are used there as initial constraints, so that the NeoKinema code considers those fault lineaments as active features, when computing the general deformation field coupled with GNSS velocities and imposed plate motions.
Computed net slip rates are in the order of 1 mm/a in average, except for a few higher values (Figure 3) in the DISS sector or in those domains where highly variable GNSS velocity directions are detected (e.g., STF, SOT, CSF, STE); in all the other sectors they are very low, suggesting that the imposed AF-EU plate motion is not capable to provide remarkable computed displacements along fault planes, or emphasizing the major effect of a temporary locking along the AF-EU subduction interface (Carafa et al., 2018).
We remark that the use of the distance weighted approach, the primary and secondary slip uncertainty computation, and the integration of plate kinematics in the resulting deformation field make the NeoKinema code a robust tool to predict tectonic setting and fault kinematics at least in southern Italy.On the contrary, faults are treated as single surfaces, so that with our geometric parameters we cannot investigate parallel structures that converge at depth to one large structure, collecting a larger amount of strain.For this reason, we do not provide consider ations on seismic hazard, also if our rake prediction is qualitatively comparable with the focal mechanisms from the European-Mediterranean Regional Centroid-Moment Tensor catalog (Pondrelli, 2002), or from Orecchio et al. (2014), Polonia et al. (2016) and Sgroi et al. (2021), in southern Tyrrhenian and Ionian Seas.
In conclusion, the approach of combining geodetic observation with geological record, using numerical methods, results in computation of a high-resolution deformation field and long-term tectonic setting in a study area, integrating space geodesy, plate motions and fault geometries.The possibility to obtain predicted fault displacements and net slip rates allows us to consider this procedure as a favorable method to improve the knowledge of fault kinematics along diffuse plate boundaries, when faults are known but their kinematic attributes elude us.The procedure level of accuracy can be improved onshore with the increase of GNSS data in a study area, and overall, with a refinement of the angular vectors describing plate motions.Furthermore, the precision increase of fault location onshore and at seafloor is fundamental to avoid artifacts in fault kinematic modeling, and to obtain robust computed crustal deformation fields with the integration of geological constraints.

Figure 1 .
Figure 1.Tectonic settings in southern Italy and 2-D computation grid.Brown lines are the compiled fault lineaments derived from the Database of Individual Seismogenic Sources (DISS) database (thick lines, with associated slip rate, DISS Working Group, 2021), and additional literature (thin lines, with no slip rate, e.g., Billi et al., 2007; Polonia et al., 2016).Solid and dashed black lines are the compression fronts and the Eurasia (EU)-Africa (AF) plate boundary from DeMets et al. (2010).Shaded areas are the Ionian Fault System (IFS) and the Alfeo-Etna Fault System (AFS).Gray and black vectors are the Global navigation satellite system velocities after Billi et al. (2023) and the computed rates relative to fixed EU, respectively.Blue and red arrows denote compression and extension areas.In the lower left inset, the second invariant of the GSRM v.2.1 geodetic strain rate model by Kreemer et al. (2014) is reported, and the solid black lines are plate boundaries from DeMets et al. (2010).EL-Eastern Lobe, WL-Western Lobe.

Figure 2 .
Figure 2. Geodetic strain rate model (this study) for the Southern Tyrrhenian and Ionian Seas.(a) 2D-dilatation and (b) second Invariant with location of the AAʹ, BBʹ, and CCʹ profiles.Gray areas are those where computed values exceed 300 nstrain/a.Brown lines are the compiled fault lineaments and the dashed black lines are the EU-AF plate boundary from DeMets et al. (2010).(c-e) sections of the AA', BBʹ, and CCʹ profiles with 2D-dilatation and second Invariant strain rate values and topographic elevation from the Emodnet (https://emodnet.ec.europa.eu/en/bathymetry)and the ETOPO 2022 Global Relief Model (www.ncei.noaa.gov/products/etopo-global-relief-model)databases.Along the profiles, a comparison with the GSRM v.2.1 model (Kreemer et al., 2014) is proposed.m.s.l.-mean sea level.

Figure 3 .
Figure 3. Observed versus predicted fault kinematics for the 160 selected lineaments at different sectors in southern Italy (see text for acronyms).Predicted net slip rate is reported in every upper inset.T-Thrust; N-Normal; R-Right Lateral; L-Left Lateral.As an example, the term N-L represents primary (N) and secondary (L) slip respectively.

Figure 4 .
Figure 4. Predicted fault rakes in southern Italy, with primary (solid color) and secondary (dashed color) slip type.Lower right panel shows model goodness, indicating percentage agreement with existing literature (gray bars).See text for acronyms.T-Thrust; N-Normal; R-Right Lateral; L-Left Lateral.