Seismogenic Structure Orientation and Stress Field of the Gargano Promontory (Southern Italy) From Microseismicity Analysis

Historical seismic catalogs report that the Gargano Promontory (southern Italy) was affected in the past by earthquakes with medium to high estimated magnitude. From the instrumental seismicity, it can be identified that the most energetic Apulian sequence occurred in 1995 with a main shock of MW = 5.2 followed by about 200 aftershocks with a maximum magnitude of 3.7. The most energetic earthquakes of the past are attributed to right-lateral strike-slip faults, while there is evidence that the present-day seismicity occur on thrust or thrust-strike faults. In this article, we show a detailed study on focal mechanisms and stress field obtained by micro-seismicity recorded from April 2013 until the present time in the Gargano Promontory and surrounding regions. Seismic waveforms are collected from the OTRIONS Seismic Network (OSN), from the Italian National Seismic Network (RSN), and integrated with data from the Italian National Accelerometric Network (RAN) in order to provide a robust dataset of earthquake localizations and focal mechanisms. The effect of uncertainties of the velocity model on fault plane solutions (FPS) has been also evaluated indicating the robustness of the results. The computed stress field indicates a deep compressive faulting with maximum horizontal compressive stress, SHmax, trending NW-SE. The seismicity pattern analysis indicates that the whole crust is seismically involved up to a depth of 40 km and indicates the presence of a low-angle seismogenic surface trending SW-NE and dipping SE-NW, similar to the Gargano–Dubrovnik lineament. Shallower events, along the eastern sector of the Mattinata Fault (MF), are W-E dextral strike-slip fault. Therefore, we hypothesized that the seismicity is locally facilitated by preexisting multidirectional fractures, confirmed by the heterogeneity of focal mechanisms, and explained by the different reactivation processes in opposite directions over the time, involving the Mattinata shear zone.

Historical seismic catalogs report that the Gargano Promontory (southern Italy) was affected in the past by earthquakes with medium to high estimated magnitude. From the instrumental seismicity, it can be identified that the most energetic Apulian sequence occurred in 1995 with a main shock of M W = 5.2 followed by about 200 aftershocks with a maximum magnitude of 3.7. The most energetic earthquakes of the past are attributed to right-lateral strike-slip faults, while there is evidence that the present-day seismicity occur on thrust or thrust-strike faults. In this article, we show a detailed study on focal mechanisms and stress field obtained by micro-seismicity recorded from April 2013 until the present time in the Gargano Promontory and surrounding regions. Seismic waveforms are collected from the OTRIONS Seismic Network (OSN), from the Italian National Seismic Network (RSN), and integrated with data from the Italian National Accelerometric Network (RAN) in order to provide a robust dataset of earthquake localizations and focal mechanisms. The effect of uncertainties of the velocity model on fault plane solutions (FPS) has been also evaluated indicating the robustness of the results. The computed stress field indicates a deep compressive faulting with maximum horizontal compressive stress, S Hmax , trending NW-SE. The seismicity pattern analysis indicates that the whole crust is seismically involved up to a depth of 40 km and indicates the presence of a low-angle seismogenic surface trending SW-NE and dipping SE-NW, similar to the Gargano-Dubrovnik lineament. Shallower events, along the eastern sector of the Mattinata Fault (MF), are W-E dextral strike-slip fault. Therefore, we hypothesized that the seismicity is locally facilitated by preexisting multidirectional fractures, confirmed by the heterogeneity of focal mechanisms, and explained by the different reactivation processes in opposite directions over the time, involving the Mattinata shear zone.

INTRODUCTION
The high seismicity rate of the Italian territory periodically causes serious damage and human losses. The investigation of active tectonic structures is of fundamental importance for the mitigation of the seismic risk. An improvement in the knowledge is then required to gain further insights on the seismogenic structures of the area. In particular, the seismogenic structures of the Gargano Promontory (hereafter GP) are still under discussion; the peculiar geographical position of the GP, the long recurrence time for large earthquakes, and the poor instrumental coverage by the RSN are the main causes for the delay in the knowledge of this seismic sector. Geodynamically, the GP and the Apulia region are considered to belong to the Adriatic microplate or Adria plate (Anderson and Jackson, 1987; Figure 1A). This microplate seems to be originated as a result of the split-up of the old collision boundaries (McKenzie, 1972;Favali et al., 1990) and played an important role in this geodynamical context, being considered as the foreland of the Alps, Apennine, and Dinarides (Di Bucci and Angeloni, 2013).
The GP represents the highest sector of the Apulian foreland (around 1,000 m above the sea level) that extends to the East toward the Adriatic Sea ( Figure 1B) and is mainly composed by slightly deformed carbonatic successions (Del Gaudio et al., 2007) as shown in the geological map in Figure 1B, where the main stratigraphical and geological features are reported. Based on previous geological maps (Bosellini and Morsilli, 2001;Cotecchia, 2014), we reported the main stratigraphical units of the area according to the typical depositional profile of the carbonate environment. Concerning the tectonic structures, many faults are recognized along the entire GP. However, the characteristics of these structures are poorly known. This is the reason why we included in the tectonic map only the main structures for which some information can be found: the Mattinata Fault (MF) (Chilovi et al., 2000), the Apricena Fault (AF) (Patacca and Scandone, 2004), the Candelaro Fault (CF) (Mongelli and Ricchetti, 1970), and the Sannicandro Fault (SF) (Salvi et al., 1999) as red lines.
The GP tectonic evolution is rather complex and is a topic of debate. Although many authors agree that the MF (Figure 1B) is the main structure involved in the development of the GP sector, its kinematics is controversial. Many authors reported that MF had a left-lateral strike-slip kinematics (Funiciello et al., 1988;Salvini et al., 1999), but several researchers considered the hypothesis of a reactivation with right-lateral motion during the late Pliocene until today (Chilovi et al., 2000;Tondi et al., 2005). In contrast, many other authors have always considered this structure as a right-lateral strike-slip structure (Piccardi, 1998). According to these authors, this structure was active since its first activation, recognized in lower Cretaceous (Piccardi, 1998), as confirmed also by the 2002 Molise earthquake (Valensise et al., 2004;Piccardi, 2005). Other studies ascribed the evolution of this area to a system of parallel E-W faults that border the promontory to the North and South (with the MF system) generating compressive stress (Brankman and Aydin, 2004) or that split this area promoting a general uplift linked to the high subduction angle of the Adria microplate (Doglioni et al., 1994). Recent results from the stress field inversion show active tectonics currently characterized by a compressive regime, with a NW-SE dominant shortening according to the NE-SW minimum horizontal stress axis orientation (Montone and Mariucci, 2016). Evidence of compression was also reported offshore, as inferred from the analysis of seismic profiles (Argnani et al., 1993(Argnani et al., , 2009Scisciani and Calamita, 2009) and onshore with the study of well logs and geological data (Bertotti et al., 1999).
The GP and surrounding regions have been hit by several strong earthquakes, as the 1627 (Mw = 6.7 ± 0.1), the May 31, 1646 (Mw = 6.7 ± 0.3), and the March 20, 1731 (Mw = 6.3 ± 0.1) earthquakes (Rovida et al., 2019). All these historical earthquakes have an estimated magnitude larger than six. In the period 1985-2013, the RSN recorded a lot of microearthquakes in the GP (de Lorenzo et al., 2017), despite the small number of seismometers located in the area. Since 2013, cooperation between the University of Bari, the INGV (National Institute of Geophysics and Volcanology), and the INFN (National Institute of Nuclear Physics) in the frame of the OTRIONS project (Tallarico, 2013) allowed a further improvement in the seismic coverage of the area, with the installation of 12 new seismic stations (Filippucci et al., 2021a). Since then, a large number of microearthquakes have been recorded. This allowed us to gain several important geophysical results for this area.
In particular, de Lorenzo et al. (2017) collected the first earthquake database by using the recordings of the network stations in the period 2013-2014 and developed a 1D velocity model for GP. Furthermore, by using the same earthquake database, the rheological behavior of GP crust (Filippucci et al., 2019a), the level of energy absorption (Filippucci et al., 2019b), and the stress field (Filippucci et al., 2020) were analyzed. In 2015, the network geometry was updated: some stations were switched off, and other stations were added in the northern part of GP, allowing a better seismic coverage and a less azimuthal gap of the recorded earthquakes. In this paper, we added to this dataset the earthquakes recorded in the period from 2015 to 2020. The earthquake seismic recordings of the period 2013-2018 were recently released on the Mendeley data repository (Filippucci et al., 2021b); the earthquakes recorded in the period 2019-2020 can be found on the EIDA platform. In this paper, using a more robust earthquake database, we extended the region of interest moving to the northern sector of the GP and tried to find the relationships between the known fault structures and the presentday microseismicity of this area. Finally, we retrieve the stress field using first-motion polarities focal mechanisms.

Seismic Database
The OTRIONS Seismic Network (OSN) ( Figure 1C) is operating in the GP and Capitanata area since April 2013. The OSN stations are integrated, in real time, with those of the RSN (Figure 1C). During the first 6 years of OSN operation, a high number of microearthquakes, below the magnitude of minimum detection of the RSN, were recorded even if the event detection was affected by several problems of data transmission and archiving (Filippucci et al., 2021a). During the first 15 months of operation, a manual work of event detection was done (de Lorenzo et al., 2017). Since July 2014 to the end of 2018, the event detection was demanded to the SeisComP3 software (Helmholtz Centre Potsdam GFZ German Research Centre for Geosciences and gempa GmbH, 2008) with manual picking revision. Waveform data from 2013 to 2018, with picking times marked, were recently released (Filippucci et al., 2021b). From the beginning FIGURE 1 | (A) Adria microplate (brown area) as reported by Mantovani et al. (2006). The red square indicates the area investigated in this paper. (B) Simplified geological map. The main stratigraphical and geological features are reported (Salvi et al., 1999;Chilovi et al., 2000;Bosellini and Morsilli, 2001;Patacca and Scandone, 2004); colors of the stratigraphic units refer to the ages of geological formations: Jurassic facies, light blue; Cretaceous facies, green; Eocene facies, orange. We reported in the map the Mattinata Fault (MF), the Apricena Fault (AF), the Candelaro Fault (CF), and the Sannicandro Fault (SF) as red lines, both continuous (that is observable on the ground surface) and discontinuous (that is hypothesized under the ground surface) and their dipping direction. (C) Enlargement on the red area of GP. The triangles refer to different seismic networks according to colors: green refers to RSN; light blue refers to OSN; yellow refers to National Accelerometric Network (RAN) managed by the Civil Protection Department (DPC).
of 2019 to the present day, data are archived and distributed by INGV-EIDA. The waveforms were downloaded from the EIDA archive, and the picking was manually revised. Furthermore, for the period 2019-2020, seismic data were integrated with data from the National Accelerometric Network (RAN) ( Figure 1C). Since these instruments operate with a triggering system, the recordings are not available for all the events, but, where possible, the azimuthal gap was significantly reduced, and the number of available first motion polarities was increased.
The collected events constitute a database of 1,894 microearthquakes. All the collected events were relocated using Hypo71 (Lee and Lahr, 1972) and a 1D velocity model computed for GP (de Lorenzo et al., 2017). Station corrections have been included in the analysis as computed by de Lorenzo et al. (2017). Hypo71 was run with different starting depths, in order to better explore the parameter space. The best location for each event was selected minimizing the following dimensionless parameter S (Michele, 2016): where RMS is the computed root mean square, RMS 0 = 1 s is a reference residual time, ERH is the horizontal error of the estimated hypocenter, ERZ is the vertical error, d min is the minimum source to the receiver distance, and Z is the foci depth. The network sensitivity allowed us to record thousand events. Earthquakes belonging to neighboring zones were discarded due to the high azimuthal gap, θ. We also recorded very small events, clustered around some active quarries and probably due to quarry explosions, which we removed from our database. Because of the GP geographical position, it is difficult to obtain earthquake location with θ less than 180 • , especially for the events located along the coastline and offshore. For this reason, we considered acceptable all the solutions with θ < 240 • . We also discarded all the solutions with RMS greater than 0.5 s. This selection process led us to reduce the database to 635 events, mapped in Figure 2 and listed in Supplementary Table 1. The quality of the database is described in Figure 3, where the frequency histograms of ERZ and ERH and the histogram of the event RMS are shown. The 99% of hypocenter locations have errors within 5 km; in particular, 93.5% of the horizontal errors and 92.4% of the vertical errors are less than 2 km, with an average horizontal error of 0.55 km and an average vertical error of 0.60 km; 73% of the residuals are less than 0.2 s with an average of 0.2 s.  These histograms indicate the high quality of the localizations. Therefore, an update of the velocity model is not necessary. Since December 2018, we estimated the vertical local magnitude (Mlv); from January 2019, the horizontal local magnitude (Ml) was considered as computed by INGV. The two major events recorded in the analyzed period are the April 23, 2017 with Mlv = 4.0 and the April 18, 2020 with Ml = 3.7 (respectively, ID = 333 and ID = 602, Supplementary Table 1). In particular, Mlv ranges between 0.1 and 2.8 and Ml ranges between 0.5 and 3.6 with the 90% of both Ml and Mlv < 2.0 (Supplementary Table 1).

Focal Mechanism Determination
We used the FocMec code (Snoke et al., 1984) to retrieve focal mechanisms from the inversion of P-wave first-motion polarities. Since we deal with earthquakes of magnitude from very low to low, the number of available polarities is very low too. Starting from a database of 635 earthquakes, we selected a total of 200 earthquakes having a minimum number of six P-wave polarities data to compute focal mechanisms.
We carried out several tests with different weights attributed to the polarity data. The FocMec code allows us to assign an error for each discrepant polarity observation according to the weight conditions imposed by the user (Snoke, 2017). Considering the small number of P-wave polarity data for each event of our database, the best choice consists of assigning a unity-weighting equal to zero; this choice is the most restrictive, since only fault plane solutions (FPS) without discrepant polarities can be accepted by the code. For each event, the code returns all the FPS compatible with the polarity data. In some cases, the acceptable FPS are so different that we cannot uniquely identify the best one, and the event must be discarded. In Figure 4, we show an example of accepted and discarded FPS, and all the 36 accepted focal mechanisms.
Following the approach described by Hardebeck and Shearer (2002), for each event, the average of the inferred fault planes is considered as the preferred one, and the uncertainty is represented by the standard deviation φ, δ, and λ of strike φ, dip δ, and rake λ. FPS is considered constrained only if the set of planes is tightly clustered around the average fault plane.  Table 1).
At the end of this procedure, we accepted only solutions with φ, δ, and λ ≤ 40 • simultaneously. By using the abovedescribed criteria, 36 well-constrained focal mechanisms were considered as acceptable. In Table 1, we reported the average value and the standard deviation for φ, δ, and λ of the 36 FPS couples (FPS 1 and FPS 2), together with the localization parameters. We also indicated the faulting type (FT) according to the classification of Zoback (1992) in the version modified by Frepoli et al. (2017) ( Table 1). The accepted solutions are shown on the GP geographical map both as beachballs ( Figure 5A),   colored according to FT (Heidbach et al., 2016) and as P-T axes ( Figure 5B) sketched according to FT too (Heidbach et al., 2016). Although the S-polarities and amplitude ratios may improve the focal mechanism results, in our case, this improvement is not possible because there are several problems affecting the S-wave data. First of all, the very small magnitude of events reflects into uncertainty that affects the estimated S-wave arrival time: S-wave onset pickings are often associated to a SAC weight from 1 to 3 that corresponds to an uncertainty ranging between 50 and 300 ms. This error bar generally impedes to pick S-polarity. Moreover S-waves are affected by the splitting of the horizontal components. This effect is due to the medium anisotropy. Unfortunately, an anisotropic model of the GP crust is not yet available, and therefore, the S-waves cannot be used to compute the focal mechanisms. In order to use the S-wave information (polarities and amplitudes), it will be therefore necessary to perform a specific study on the anisotropy of GP and Northern Apulia, which is beyond the scope of this paper.

Error Propagation and Validation
The use of a 1D velocity model can bias the earthquake locations, especially in the complex tectonic GP area. In order to account   Table 2). The distances ED and VD among hypocenters, due to the error bar of the velocity model, are comparable with the localization errors ERH and ERZ (Figure 3).
We also recomputed focal mechanisms by using the three hypocenters for each event (the MIN, the MAX, and the AVG) and plotted the results (Supplementary Figures 1-6). It can be observed that the three focal mechanisms for each event do not show significant variations, indicating stable results. It is worth to note that, since we set FocMec (Snoke et al., 1984) to return only solutions without discrepant polarities, errors in the velocity model can move the polarity to be discrepant and, in some solutions, can be not accepted (Ns = 15,16,17,30,32,and 34 in Supplementary Figures 3, 5, 6).

The Stress Tensor Inversion
We used the FMSI package (Gephart, 1990) to perform different tests on our focal mechanism database. By performing this inversion, we estimated four of the six components that describe the stress tensor; the three eigenvectors that represents the principal stress direction (σ 1 , σ 2 , and σ 3 , respectively, maximum, intermediate, and minimum compressive stress direction) and the stress ratio, R = (σ 2 − σ 1 )/(σ 3 − σ 1 ), a dimensionless scalar quantity indicating the shape of the stress tensor.
The data set consists of 36 high-quality focal mechanisms: 7 TF and 15 TS, 4 NF, 1 NS, 5 SS, and 4 U ( Figure 5A and Table 1). To constrain the results obtained by the stress tensor inversion, we assigned a weight for each nodal plane. We used the criteria described in Filippucci et al. (2020) and adapted the quality factors of FPFIT (Reasenberg and Oppenheimer, 1985) to the FocMec output. The fitting quality Q f quantifies how the solution fits polarity data. In this case, all the 36 returned focal mechanisms are perfectly compatible with the polarity observations, so we assigned a default Q f = 1 to all the events. The nodal plane quality factor Q p quantifies the quality of the retrieved nodal plane: We used the corresponding standard deviation of these three parameters to describe the range of perturbations of φ, δ, λ. We assigned Q p1 = Q p for FPS 1 and Q p2 = Q p for FPS 2, where FPS are the nodal planes identified by FocMec ( Table 2). Q p12 is the total weight assigned to the two FPS: Q p12 = 3 for the AA class; Q p12 = 2 for AB, BA, BB, AC classes; Q p12 = 1 for BC, CB, CA classes was assigned to each event. We also considered the number of polarity data, N pol , to assign the weight Q pol to the focal mechanism. Since the number of available observations is very low, we assigned Q pol = 0 to focal mechanisms computed with N pol = 6 and N pol = 7. Then, by scaling for the maximum number of N pol per event, we assigned the following weights: Q pol = 1 for 8 ≤ N pol ≤ 10; Q pol = 2 for 11 ≤ N pol ≤ 15; Q pol = 3 for 16 ≤ N pol ≤ 20, Q pol = 4 for 21 ≤ N pol ≤ 25; Ns, identification number of the focal mechanism; Q p1 and Q p2 , the quality factors relative to the FPS 1 and FPS 2, respectively; N pol , the number of polarities; Q p12 , the weight assigned to Q p1 and Q p2 ; Q pol , is the weight assigned to N pol ; W TOT , the total weight used in the FMSI inversion.
Q pol = 5 for 26 ≤ N pol ≤ 30; Q pol = 6 for 31 ≤ N pol ≤ 35. W TOT = Q f + Q p12 + Q pol is the total weight used in the FMSI inversion for each focal mechanism. The frequency histogram of foci depths Z ( Figure 2B) shows one main seismogenic layer at 14 < Z < 30 km. Starting from this observation, we grouped focal mechanisms according to depth. A first inversion was performed with 23 events (listed in Table 2), but the misfit angle (Misfit = 8.56 • in Figure 5C) was too high to accept this solution, according to the homogeneous criteria proposed by Lu et al. (1997). We then performed a second run by excluding three events (Ns = 7, 11, 26 in Table 2) located borderline of GP. We obtained Misfit = 6.74 • (Figure 5D), which indicates a homogeneous compressive stress field with a subvertical σ 3 and a maximum horizontal compressive stress, S Hmax , trending NW-SE according to the Italian Stress Map for GP (Carafa and Barba, 2013;Carafa et al., 2015b;Montone and Mariucci, 2016). The misfit angle (Gephart and Forsyth, 1984) is computed through an angular rotation given by the angular difference between the observed slip direction on a fault plane and the shear stress on that fault plane derived from a given stress tensor; the stress tensor orientation that provides the average minimum misfit is assumed to be the best stress tensor for a given population of focal mechanisms (Frepoli et al., 2011).
Finally, in order to observe the distribution of the stress field within the GP area, we calculated the orientation of S Hmax (Figure 5E) by interpolating all the 36 focal mechanisms; to this end, we used the method described in Carafa and Barba (2013) and the SHINE code (Carafa et al., 2015a).

RESULTS AND DISCUSSION
The seismicity collected in the period from April 2013 to June 2020 indicates a dense and sparse seismic activity along the whole GP. This seismicity is not related to any seismic sequence. The recorded seismicity (Figure 2A) is distributed throughout the GP following roughly a NE-SW direction, with deepest events located in the NE sector. From the frequency histogram of the hypocenter's distribution (Figure 2B), a seismicity peak at 18 < Z < 27 km can be observed. At these depths, earthquakes are related to a diffuse seismicity without a clear time and/or spatial relation, with relative distances also greater than tens of km. We relocated with HypoDD (Waldhauser and Ellsworth, 2000) this deeper seismicity, but no significative change was observed, and we discarded this result.
We investigated the possible bias on locations and focal mechanisms of the 36 final events, due to the error bar of the velocity model of de Lorenzo et al. (2017). These tests led us to assess that the velocity model computed by de Lorenzo et al. (2017) gives rise to reliable results (Supplementary Table 2 and Supplementary Figures 1-6).
By dividing events into four depth classes, a substantial absence of shallow seismicity is visible in the NE of GP (Figure 6). For 0 < Z < 20 km, earthquakes are mainly located S of the MF zone and NW of the investigated area; for 20 < Z < 30 km, the seismicity is distributed along the whole GP; for Z > 30 km, the seismicity is located only NE in the GP. In Figure 6, it can be seen how the seismicity is distributed both S and N of the MF system. Different kinematics can be observed at different depths ( Figure 5). Normal to normal-strike events are detected at S in the GP, at shallow depths with different T-axis orientations. Twentytwo of the 36 focal mechanisms present thrust and thrust-strike kinematics, and many of them show a P-axis trending NW-SE (Figures 5B, 6). In Figure 5A, five strike-slip events, located W of the GP along the MF at Z = 20 km (except for the shallower event Ns = 18), show a right-lateral mechanism. The P-axes of these five strike-slip solutions (Figure 5B), except for the event Ns = 1, which is slightly rotated, are trending NW-SE similar to the majority of the thrust and thrust-strike P-axes. According to the results, we agree with the interpretation FIGURE 6 | Geographical maps of epicentral distribution of all 635 earthquakes and 36 focal mechanisms, grouped by four depth ranges, as superimposed on each map. Epicenters are colored according to depth. described the different depth ranges. The black line is the MF System as reported by Zunino et al. (2012). Arrows represent the P-T axes with colors referring to kinematics, as shown in Figure 5. Hypocenters are plotted as black circles in the vertical cross-section maps. Above each cross-section, the topographic profile is reported. In each cross-section map, the focal mechanisms are projected in cross-section. Beachball colors are described in Figure 5 caption. Focal mechanisms dimensions are scaled according to magnitude ranges (small for Ml < 3; big for Ml ≥ 3). Green line represents the MF trace reported in depth. The surrounding area highlighted in light green represents the rock volumes involved in the MF shear zone. Red lines are for the Apricena AF, the Candelaro CF, and the Sannicandro SF Faults, dashed lines indicate that no information is available for the fault dip and/or for the fault depth.
of the MF as a right-lateral strike-slip structure. Moreover, the inferred maximum horizontal compressive stress, S Hmax , is SE-NW oriented (Figures 5C,D), and it is compatible with the recognized dextral strike-slip tectonics in the foreland sector of the Italian Apennines. As known from experiments, a strikeslip fault can produce a compression or an extension of the FIGURE 8 | The same as Figure 7 for the N-S cross-sections. All cross-sections have 50 km length and 8 km semi-width per side.
Frontiers in Earth Science | www.frontiersin.org surrounding rock volumes according to its relative sense of shear; regarding the MF case, it is observed that a W-E right-lateral strike-slip structure can cause a NW-SE compression that occurs along the NE-SW-oriented strikes (Piccardi, 1998(Piccardi, , 2005Chilovi et al., 2000;Valensise et al., 2004;Tondi et al., 2005).
As shown in Figure 5A, nodal planes trend approximately NE-SW, with low to moderate dip angles, underlining a general compression toward NNW-SSE with NW-SE trending P-axes ( Figure 5B) consistent with the data reported in the Italian Stress Map (Montone and Mariucci, 2016). Moreover, the stress inversion (Figures 5C,D) indicates a deep homogeneous compressive stress field with S Hmax oriented NW-SE. This trend is also observed by computing the stress field in any point of the GP area. As shown in Figure 5E, the resulting S Hmax orientation seems to be constant along the whole promontory, highlighting this NW-SE regional compression from the SW of GP to the Adriatic Sea, also according to the Italian Stress Map data (Montone and Mariucci, 2016) interpolated and shown in Carafa and Barba (2013).
Different kinematics are observed W of GP in contrast with previous studies (Carafa and Barba, 2013;Montone and Mariucci, 2016). In fact, from focal mechanism interpolation (Figure 5E), a change in the S Hmax orientation is observed W of GP. According to Zoback (1992), this spatially concentrated variation of the maximum horizontal compressive axes may be associated to a "second-order" stress pattern due to a local perturbation of the stress field whose causes are still unknown.
The observed relation between seismicity distribution and faulting kinematics seems to suggest three hypothesis that we will investigate in the next subsections by performing vertical crosssections: (1) the presence of a seismogenic layer trending toward NE (SW-NE cross-sections); (2) the role of the MF structure in the observed seismicity pattern (S-N cross-sections); (3) the presence of a seismogenic layer dipping toward NW (SE-NW cross-sections).

SW-NE Cross-Sections
We considered five SW-NE cross-section lines (Figure 7), orthogonal to the geomorphological evidence of compression (Bertotti et al., 1999) and parallel to the direction of surface heat flow, which decreases toward NE (Filippucci et al., 2019a) to investigate the existence of a seismogenic layer trending toward NE. Each cross-section line has a horizontal box-section with a length of about 52 and 16 km of width. Each earthquake and each focal mechanism falling into the box-section is plotted in the vertical cross-section map.
From cross-sections 3 and 4 (Figure 7), it is evident that a great number of events, located at Z > 15 km, occur along a slightly tilted surface that seems to dip toward the NE. It is worth to note that this deepening of seismicity probably shows the trend of the Adriatic Moho, whose depth is reported to range between 35 and 40 km (Chiarabba et al., 2005;Piana Agostinetti and Amato, 2009). Doglioni et al. (1994) identified the high angle of the subduction plane of the Adria plate below the Apennine chain as the cause of a partial rise of the Moho and of the general uplift of the GP sector. However, the majority of focal mechanisms projected in the vertical cross-sections show nodal planes NNW-SSE oriented orthogonally to the SW-NE observed seismicity trend; therefore, these plane solutions do not seem to highlight any seismic structure trending SW-NE. SW in the GP, normal focal mechanisms and the related shallow seismicity seem to highlight a SW dipping structure, which, however, does not seem to be related to the CF (red dashed line in Figure 7), although these events seem to indicate that some stress release involves this fractured volume of the upper crust.
It can be also observed, in the vertical cross-sections 3 and 4 (Figure 7), that in the NE of GP, seismicity is not present above 20 km and that earthquakes seem to concentrate around the area of the MF system, whose depth is estimated to range between 15 and 20 km (Piccardi, 1998;Valensise et al., 2004;DISS Working Group, 2018). For the MF shear zone, we assumed a depth of 20 km and a width of 4 km (Figure 7, green line and rectangle in vertical cross-sections) smaller than that assumed by Salvini et al. (1999) in order to reduce over-interpretation errors.

S-N Cross-Sections
We considered five SN cross-section lines oriented orthogonal to this structure (Figure 8) to investigate the role played by the MF zone in seismicity pattern and stress field. Excluding crosssections 1 and 5 that border the MF zone, from W toward E, there is a general increase in seismicity distributed along two FIGURE 10 | The same as Figure 7 for SE-NW cross-sections. All cross-sections have 69 km length and 8 km semi-width per side.
Frontiers in Earth Science | www.frontiersin.org main probable alignments. In these three cross-sections, a lack of shallow seismicity N to the MF zone is visible (Figure 8, green line and green rectangle in the vertical cross-sections). Deep seismicity is shown in all cross-sections, but it is dominant in cross-sections 3 and 4 (Figure 8). In fact, all earthquakes appear to be located along a submerging N surface with a slightly more inclined angle moving toward E.
The lack of shallow seismicity in the NE part of GP starting from the MF zone corroborates the hypothesis of a "natural barrier" role of the MF zone in the seismicity pattern. At the S of MF, earthquakes with normal to normal-strike faulting are observed. In cross-sections 2 and 3 (Figure 8), focal mechanisms do not allow to assert whether this seismicity belongs to CF (red dashed line in Figure 8). No events are localized along AF and SF.
Along the eastern segments of the MF at 15 < Z < 20 km (Figure 8, cross-section 3), it is possible to observe a large concentration of seismicity. This could be a hint of the presence of a splay of this deeper surface or could be caused by an area of extreme weakness that can facilitate the stress release. This second hypothesis can be confirmed by the heterogeneity of the focal mechanisms. The cross-section overlapping in Figure 9 permits to better identify the seismic horizons, the kinematic heterogeneity in the proximity of the MF structure, and the good correlation between seismicity distribution and focal mechanisms. In Figure 9A, the general seismicity trend is shown, revealing the presence of a deep surface. In fact, a great number of events appear to be aligned in depth along a slightly tilted surface. In Figure 9B, we plotted only focal mechanisms falling in proximity of the MF structure in order to highlight the kinematic heterogeneity in correspondence of the MF shear zone (Figure 9, green rectangle). The observed heterogeneity validates the idea of the existence of fractured rocks that could break with faulting kinematics, which are different from the strike-slip faulting ascribed to MF. In this frame, a focal mechanism of events located in and around the MF rectangle seems to show right-strike-slip kinematics and hints of partial reuse of the MF system in the compressive tectonic domain.

SE-NW Cross-Sections
We considered four cross-section lines SE-NW oriented (Figure 10), perpendicular to those of Figure 7 in order to investigate the presence of a seismogenic layer dipping toward NW. Deep events are few and quite widespread along crosssection 1 (Figure 10) so any interpretation is not possible. In cross-section 2 (Figure 10), deep seismicity seems to follow the trend previously observed in Figures 8, 9, highlighting a slope slightly dipping toward NW. Along the MF zone (Figure 10 green line and green rectangle in cross-sections), an abrupt interruption of the shallower seismicity and a greater concentration of events at Z > 15 km is observed. These features are also shown in cross-section 3 (Figure 10) where the events, close to the MF, are more concentrated at Z = 20 km. Moving from cross-sections 2-4 (Figure 10), the maximum foci depth increases, highlighting the deepening of the Moho toward the Adriatic Sea, also shown in cross-sections 3 and 4 of Figure 7. The seismicity pattern seems to confirm the presence of a low-angle surface slightly dipping N-NW that hosts this seismicity. This hypothesis seems to find an additional constraint in the plotted thrust and thrust-strike focal mechanisms (Figure 10, cross-sections 2, 3, and 4) that show one nodal plane in agreement with this hypothesized structure.
In Figure 11, the vertical cross-sections 2, 3, and 4 of Figure 10 are overlapped in order to highlight the presence of a highangle shallower seismogenic layer and a low-angle deeper one, which trends SW-NE and slightly dips NNW. In Figure 11A, all focal mechanisms are plotted, excluding normal and normalstrike events. The seismic/aseismic transition seems to occur  along a curved surface that is sketched in Figure 11B. Two interpretations can be advanced: the first one considers a lowangle deep surface (red solid line in Figure 11B), in agreement with kinematic data (red boxes, Figure 11B), with the seismicity distribution and with the general compressive/transpressive kinematic with S Hmax NW-SE oriented; the second one considers as a seismic/aseismic transition the surface (orange dashed line in Figure 11B sketched by joining the focal mechanisms grouped in the orange boxes of, Figure 11B) around the MF shear zone (green rectangle in Figure 11B). This second interpretation agrees with a change in the dip angle of this transition surface and also agrees with the hypothesis of a local release of microearthquakes, rupturing on planes with different kinematics. The heterogeneity of the obtained focal mechanisms around the MF can be justified by the sequence of different reactivation processes over time. This observation is also confirmed by the seismicity distribution, which clusters around the MF (Figures 9, 11).

CONCLUSION
New insights of the GP seismicity have been presented in this paper. We used the micro-earthquakes recorded from April 2013 to June 2020, with the aim of understanding a possible correlation between the low-magnitude seismicity and the known seismogenic structures. From the seismicity pattern analysis, we found out that the whole crust is seismically involved up to about Z = 40 km with an earthquake concentration at 14 < Z < 37 km. We interpret the GP seismicity as generated by a deep compressive stress regime with S Hmax oriented NW-SE with a deep faulting NE-SW oriented. This interpretation is verified by the focal mechanisms, the stress field inversion, and stress field interpolation. Our results agree with the Italian Stress Map (Montone and Mariucci, 2016) and do not agree with the extensional present-day interpretation of Billi et al. (2007) and with the surface fault lines, which the contractional tectonic evolution of the GP is attributed to Bertotti et al. (1999).
We hypothesize the presence of a low-angle seismogenic surface trending SW-NE and dipping SE-NW. This surface could also extend in the Adriatic Sea following the same trend, according to the focal mechanisms from Pondrelli et al. (2006), to the Italian Stress Map (Montone and Mariucci, 2016) and to the interpolated stress field computed by Carafa and Barba (2013). Moreover, the seismicity localized along the Adriatic sector (Giardini et al., 2014) seems to connect the GP with the Dinaric front, spreading with a SW-NE trend along the Gargano-Dubrovnik Fault (Battaglia et al., 2004; Figure 12B). Based on GPS data (Oldow et al., 2002), Adria is supposed to be divided into two microplates, Northern and Southern Adria plates (NAd and Sad, respectively), and connected by the G Fault (Figure 12).
Regarding the role of the MF system, it cannot be stated that the whole MF is an active fault, but that some sectors of this structure are releasing the present stress field. Regarding the lack of shallow seismicity in the North-Northeast area of the GP, it can be hypothesized that deep seismicity (between 15 and 30 km) may be favored by the presence of weak preexisting fractures, lubricated by the circulation of deep fluids, as it happens in other Italian seismogenic areas (Miller et al., 2004). The presence of fluids in the upper crust was supported by geochemical analysis (Salvi et al., 1999) and resistivity imaging (Tripaldi, 2020) in proximity to CF and could explain the seismicity observed near the CF system. To disentangle this issue, which has important implications in seismic hazard, deep geophysical imaging and thermal-rheological modeling, which includes fluid convection, should be conducted.
Obviously, we cannot exclude that the lack of shallow seismicity at N will produce an earthquake energetic enough to affect the entire crustal volume, but at the same time, we cannot assert it. Some of the strongest historical earthquakes that affected this area, as the 1414 earthquake (Mw 5.8 ± 0.46), the 1646 earthquake (Mw = 6.72 ± 0.25) and also the 1889 and the 1893 events (Mw = 5.47 ± 0.12 and 5.39 ± 0.19, respectively) (data from Rovida et al., 2019), were located in this sector of GP without any information about foci depths. The September 30, 1995 earthquake (Mw = 5.15 ± 0.10) occurred at a depth of about 27 km (Del Gaudio et al., 2007), with a controversial kinematic, whether it is reverse faulting (Pondrelli et al., 2006) or E-W dextral strike-slip (Del Gaudio et al., 2007), both compatible with the stress field computed in this work.
Possible future developments of this research could be the integration of geodetic information, InSAR data (Atzori et al., 2007), and GPS data (Carafa et al., 2020), with seismic data to estimate the strain rates of the area. Moreover, new seismic data analysis and a more extended area under study around the GP could be useful to investigate how the GP is inserted in the Italian and Adriatic geodynamic context.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. These data can be found here: Seismic waveforms for the period 2013-2018 can be found at: https://doi.org/10.17632/7b5mmdjpt3.3, while seismic waveforms for the period 2019-2020 can be found at https://www.orfeus-eu.org/data/eida/nodes/INGV/.

AUTHOR CONTRIBUTIONS
SM organized the database and produced figures and tables. SL, MF, AF, and SM contributed to earthquake localizations. SL, MF, SM, and PP contributed to focal mechanism solutions. SM and PP performed the stress inversion. AT supervised and financed the work. MF and SM wrote the original draft. All authors contributed to conception and design of the study, manuscript revision, read, and approved the final version.

FUNDING
The computational work has been executed on the IT resources of the ReCaS-Bari data center, which has been made available by two projects financed by the MIUR (Italian Ministry for Education, University and Re-search) in the "PON Ricerca e Competitività 2007-2013" Program: ReCaS (Azione I-Interventi di rafforzamento strutturale, PONa3_00052, Avviso 254/Ric) and PRISMA (Asse II -Sostegno all'innovazione, PON04a2_A). This work was partially supported by Project PRIN n. 201743P29 FLUIDS (Detection and tracking of crustal fluid by multiparametric methodologies and technologies).

ACKNOWLEDGMENTS
We would like to thank Federica Ferrarini and two reviewers for their suggestions, which greatly improved the manuscript. We would also like to thank Michele C. M. Carafa for stimulating discussion and for assistance using the software SHINE. Figures were obtained by employing the GMT freeware package by Wessel and Smith (1998) and by employing the QGIS free software, v 3.8 (www.qgis.org).