Crustal Structure of the Seismogenic Volume of the 2010–2014 Pollino (Italy) Seismic Sequence From 3D P- and S-Wave Tomographic Images

A tomographic analysis of Mt. Pollino area (Italy) has been performed using earthquakes recorded in the area during an intense seismic sequence that occurred between 2010 and 2014. 870 local earthquakes with magnitude ranging from 1.8 to 5.0 were selected considering the number of recording stations, the signal quality, and the hypocenter distribution. P- and S-wave arrival times were manually picked and used to compute 3D velocity models through tomographic seismic inversion. The resulting 3D distributions of VP and VS are characterized by high resolution in the central part of the investigated area and from surface to about 10 km below sea level. The aim of the work is to obtain high-quality tomographic images to correlate with the main lithological units that characterize the study area. The results will be important to enhance the seismic hazard assessment of this complex tectonic region. These images show the ductile Apennine platform (VP = 5.3 km/s) overlaying the brittle Apulian platform (VP = 6.0 km/s) at depth of around 5 km. The central sector of the area shows a clear fold and thrust interface. Along this structure, most of the seismicity occurred, including the strongest event of the sequence (MW 5.0). High VP (>6.8 km/s) and high VP/VS (>1.9) patterns, intersecting the southern edge of this western seismogenic volume, have been interpreted as water saturated rocks, in agreement with similar geological context in the Apennines. These fluids could have played a role in nucleation and development of the seismic sequence. A recent study revealed the occurrence of clusters of earthquakes with similar waveforms along the same seismogenic volume. The hypocenters of these cluster events have been compared with the events re-located in this work. Jointly, they depict a 10 km × 4 km fault plane, NW-SE oriented, deepening towards SW with a dip angle of 40–45°. Instead, the volume of seismicity responsible for the ML 4.3 earthquake developed as a mainshock-aftershock sequence, occurring entirely within the average-to-low VP/VS Apennine platform. Our results agree with other independent geophysical analyses carried out in this area, and they could significantly improve the actual knowledge of the main lithologic units of this complex tectonic area.


INTRODUCTION
Mt. Pollino area (Italy) is located at the transition between the southern edge of the Apennines chain and the Calabrian subduction forearc ( Figure 1A). It is one of the most structurally complex tectonic regions in the intra Apennine seismogenic belt (Michetti et al., 1997;Cinti et al., 2002). The area is considered a seismic gap in the Apennine chain because no earthquakes of magnitude M > 6 occurred in historical time, in striking contrast with the nearby north and south sectors where the occurrence of many strong earthquakes is well documented ( Figure 1B). During the last centuries, the seismicity of Pollino area has been of moderate magnitude (M L < 6.0, CPTI15 database; Rovida et al., 2016), but paleo-seismological studies have recognized active faults south-east of the study area that generated at least two earthquakes of magnitude larger than M L 6.5 in the last 10,000 years (Cinti et al., 2002;Cheloni et al., 2017). The geology of the area is characterized by the overlay of a low strain-rate extensional deformation occurred during Early Pleistocene (D'Agostino et al., 2011) over pre-existing strikeslip Late Pliocene-Early Pleistocene structures (Ghisetti and Vezzani, 1982;Schiattarella, 1998). The current NE-SW extensional stress regime (Totaro et al., 2016;Napolitano et al., 2021) is driven by the simultaneous retrograde S-Eward motion of the subduction zone and by the opening of the Tyrrhenian back-arc basin (D'Agostino et al., 2011). The result of such lithospheric evolution is a complex structure characterized by W-and E-dipping sub-parallel normal faults in the upper crust (above 15 km depth) and by strike-slip faults below 15 km depth (Totaro et al., 2015;Brozzetti et al., 2017;Ferranti et al., 2017;Napolitano et al., 2021).
The recent seismicity of the study area has been characterized by low and moderate magnitude earthquakes. A seismic sequence followed the M W 5.6 earthquake that occurred on September 9, 1998, in the Mercure basin (Guerra et al., 2005;Brozzetti et al., 2009). Between 2010 and 2014 a seismic sequence of more than 10,000, low-to-moderate earthquakes occurred in the study area, on the western side of the Pollino Massif, with a time-varying rate of seismicity (Passarelli et al., 2015;Totaro et al., 2015). The sequence developed within two seismogenic volumes and showed a combination of both swarm-like and classic mainshockaftershock sequence (Passarelli et al., 2015;Napolitano et al., 2021). Two main events occurred in 2012: a M L 4.3 in the eastern volume on 28 May, and a M L 5.0 (the largest magnitude earthquake of the sequence) in the western volume on 25 October (yellow stars in Figure 1A). These events were characterized by normal fault mechanism (black beach balls in Figure 1A) in perfect agreement with the local extensional stress field of the area ( Figure 1C, modified after Napolitano et al., 2021). In addition, a slow earthquake, located at the western seismogenic volume of the sequence, occurred from mid-2012 to mid-2013, releasing aseismically an amount of strain energy equivalent to M W 5.5 (Cheloni et al., 2017).
Many studies have been carried out to assess the seismic hazard of the Pollino area seismic gap. Earthquakes from the 2010-2014 sequence have been analyzed by many authors to FIGURE 1 | (A) Map of the study area, located between the Southern Apennines and the Calabrian subduction fore arc (CFA). Seismic events analyzed are shown with red circles. Yellow stars represent the two largest events of the 2010-2014 Pollino sequence. Beach balls represent their focal mechanism. Permanent and temporary seismic stations operated by INGV (blue triangle) and UNICAL (cyan circle) and shared between INGV and UNICAL (blue triangle in cyan circle) and by GFZ (green diamond) are shown. (B) The historical earthquakes of M L >6 (CPTI15, Rovida et al., 2016) are shown with blue squares. The target area, clearly a seismic gap, is enclosed in a magenta square and the orange circle is the location of the swarm area. (C) The local stress field, modified after Napolitano et al. (2021). provide new insights about the structure of the seismogenic faults (Totaro et al., 2015;Brozzetti et al., 2017) and to study site effects (Napolitano et al., 2018). Part of them have been used by Napolitano et al. (2021) to relocate with high precision clusters of events with similar waveforms to compute focal mechanisms of small earthquakes, to obtain a very local stress field, and to investigate the swarm evolution in space and time. Geological-structural, seismic, and geophysical data have been combined by Ferranti et al. (2017) to characterize the seismotectonic frame of the study area from surface to deep crust. The presence of fluid-filled connected faults has been inferred from 2D imaging of scattering and absorption contributions to the total attenuation of high frequency coda waves Sketsiou et al., 2020). The codanormalization method has been applied to perform a 3D attenuation image of fluid storage and tectonic interactions across the Pollino fault network (Sketsiou et al., 2021).
A detailed study of the elastic properties of the crust is useful to assess the seismic hazard of complex tectonic areas. When many well recorded earthquakes are available, this goal can be achieved by passive seismic tomography, based on the use of P-and S-wave travel times. The technique is particularly effective in areas characterized by appropriate distributions of hypocenters and seismic stations. Until now, the Pollino area has been only marginally included in large scale tomographic studies, whose main target area was the Southern Apennines-Calabrian arc region (Barberi et al., 2004;Orecchio et al., 2011). These papers analyzed earthquakes occurred before the 2010-2014 Pollino seismic sequence. Barberi et al. (2004) used earthquakes recorded from 1978 to 2001 and 84 artificial sources, while Orecchio et al. (2011Orecchio et al. ( ) used data recorded between 1981Orecchio et al. ( and 2008. In both papers, the authors performed the inversion using a horizontal and vertical grid spacing of 40 and 10 km, respectively. More recently, Totaro et al. (2014) used earthquakes occurred from 1981 to October 2012, partially including events of the 2010-2014 seismic sequence, to perform a seismic tomography at regional scale. Also in this case, the authors focused on the likely transition from the thinner Tyrrhenian crust to the thicker crust beneath Calabria. Finally, Pastori et al. (2021) provided a new local 1D velocity model with the goal of performing refined 1D and double-difference locations using the events from the last Pollino sequence.
Our work aims to provide a highly detailed 3D tomographic image of the shallow crust in the Pollino area, taking advantage of the 2010-2014 seismic sequence. We mean to enhance our knowledge about the elastic properties of the seismogenic volume responsible for the long-lasting sequence and to investigate the possible role of fluids in triggering the sequence, as suggested by many authors (Passarelli et al., 2015;Sketsiou et al., 2020;Sketsiou et al., 2021). To reach this goal, we preliminarily checked the quality of the available waveforms recorded at the permanent and temporary stations installed in the area between 2010 and 2014. For each of these waveforms, we manually picked P-and S-wave arrivals and checked their quality. A preliminary location of these events has been obtained using one of the large scale models available for Southern Italy. Using this reliable dataset, through a statistical procedure, we obtained a 1D model, used as starting model of the 3D tomography which simultaneously inverts the arrival time of body waves for both velocity model and earthquake location parameters. Synthetic tests were carried out to validate the results. Then, the results were compared with other independent geological and geophysical information available for this area. A reliable interpretation of the properties of the shallow crust beneath the Pollino area might be an essential starting point for future more detailed analysis about seismogenesis and evolution of the 2010-2014 sequence.

DATA
We selected 870 earthquakes occurred in the Pollino area, within a volume of 100 × 120 × 25 km 3 ( Figure 1A). These events, characterized by magnitude between 1.8 and 5.0, were recorded between January 30, 2010, and November 25, 2014. We manually picked P-and S-wave arrival times of each waveform to guarantee the robustness and reliability of the analyses to be carried out. We collected 9981 P-and 6862 S-wave arrival times recorded by 29 seismic stations ( Figure 1A) operated in the area by three institutes: Università della Calabria (Unical), Istituto Nazionale di Geofisica e Vulcanologia (INGV), and Deutsches GeoForschungsZentrum of Potsdam (GFZ) (Passarelli et al., 2012;Margheriti et al., 2013). Weights from 0 to 4 were assigned to each picking based on its accuracy ( Table 1).
We applied the modified Wadati diagram method to evaluate the picking accuracy of the P-and S-wave arrival times (Chatelain, 1978). Furthermore, the Wadati diagram provides an estimate of the average V P /V S value for the study area. The average V P /V S , estimated considering all points within one standard deviation from the best fit line, results to be V P V s 1.786 ± 0.004 (slope of the red line in Figure 2A). The V P /V S obtained from this analysis is very similar to those of previous studies characterized by different scales and datasets (Barberi et al., 2004;Piana Agostinetti and Amato, 2009;Ferranti et al., 2017;Pastori et al., 2021).
As a first step, we performed a preliminary hypocenter location of the initial dataset using the probabilistic approach implemented in the software NonLinLoc (Lomax et al., 2000). This preliminary location takes advantage of the 1D P-wave velocity model by Barberi et al. (2004) and of the V P /V S estimated from the modified Wadati diagram ( Figure 2A). We established a lower threshold of 6 direct wave arrival times per event including at least 2 S-wave picks and a root-mean-square (RMS) smaller than 0.6 s to retain location results. We combined seven location parameters ( Figure 2C) in order to evaluate the earthquake location quality factor following the formula proposed by Michele et al. (2019) (q f , Figure 2B). We obtained 779 A-class and 75 B-class locations, which confirm the high quality of the dataset and the reliability of signal pre- The seven location parameters used to obtain q f : Err (km) is the error on horizontal (red) and vertical (green) hypocenter locations, azimuthal gap in degrees (Gap), root mean square error in seconds (RMS), number of phases (Nph), distance between the probability density function (pdf) expected value and its maximum likelihood (LocDist), and radius of the volume described by the scatter points radius of pdf (PDFRad).
Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 735340 4 processing. The quality assessment allowed us to define a highquality dataset comprising 854 well located earthquakes, for which 9643 P-and 6649 S-wave arrival times are available, useful for the following analysis.

TOMOGRAPHIC METHOD AND INVERSION STRATEGY
In this section, we describe the procedure that was applied to obtain the 3D tomographic inversion and explain how we tuned the selected parameters. We first describe the selection of the most reliable 1D model used as starting model for the tomographic inversion, and then we focus on the inversion strategy.
A reliable final 3D velocity model must necessarily start from a trustworthy 1D model (Kissling et al., 1994). We chose the P-wave velocity model by Barberi et al. (2004) as the starting reference 1D velocity model for the preliminary hypocenter location and data pre-processing. Then, we applied a statistical procedure that allows to find the most reliable local reference 1D velocity model to be used as an input in the 3D tomography of the study area (Vanorio et al., 2005;De Matteis et al., 2010). Following this procedure, we generated 120 1D P-wave models characterized by velocity values randomly chosen in the range ±500 m/s around the value of the reference model for each layer ( Figure 3A). For each of these 120 models, we performed a preliminary relocation of the dataset using NonLinLoc (Lomax et al., 2000) and then we performed a 3D tomographic inversion using the same velocity model (Latorre et al., 2004). We then compared the RMS time residuals reduction between the first and 12th iteration (red and light blue points in Figure 3B). Even though the starting RMS values span from 0.17 s up to 0.44 s, all the final solutions reach RMS time residuals smaller than 0.1 s. We selected the 47 velocity models whose starting RMS time residuals were smaller than 0.25 s. The velocity at each layer was averaged and a new local reference 1D P-wave velocity model was obtained (solid red line in Figure 3A). This new velocity model turns out to be faster (0.23 km/s on average) than the starting model by Barberi et al. (2004) for almost any layers and agrees with the most recent 1D model obtained by Pastori et al. (2021). The S-wave velocity model has been obtained dividing the P-wave velocity by V P /V S 1.786, estimated from the modified Wadati diagram (Figure 2A).
We used a linearized, iterative travel time tomography which simultaneously inverts the arrival time of body waves for both velocity model and earthquake locations parameters (Latorre et al., 2004). This method has been used in many applications in tectonic regions and volcanic and geothermal areas using data from passive (e.g., Vanorio et al., 2005;Amoroso et al., 2014;Amoroso et al., 2018) and active surveys (De Landro et al., 2017;De Landro et al., 2020). The hypocenter and station distributions FIGURE 3 | (A) 120 randomly generated initial 1D P-wave velocity models (grey lines). Blue line represents the starting P-wave velocity model by Barberi et al. (2004), black lines are the upper and lower velocity bounds to perform the random selection. Red line is the new local P-wave velocity model achieved in this work and used as starting velocity model for the 3D velocity tomography. Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 735340 allow us to set the mesh spacing of the inversion grid to 5 × 5 × 1 km 3 . A finer grid of 0.5 × 0.5 × 0.5 km 3 was set to compute the theoretical travel times at each station through a finite-difference solution of the eikonal equation (Podvin and Lecomte, 1991). Travel times were recalculated by numerical integration of the slowness along the ray path. The model parameters (velocities, hypocenter coordinates, and origin time) were inverted using the least squares root (LSQR) method of Paige and Saunders (1982). The damping parameter was set to 0.3 through the empirical approach of the trade-off curve between the model parameters (P-and S-wave model, spatial hypocenter location, and origin time) and the data variances L-curve (Supplementary Figure S1). The velocity distribution in a continuous medium is described by a trilinear interpolating function based on a grid of regularly spaced nodes.
The derivative weight sum (DWS), that provides a measure of the density of the rays within each grid cell, and the checkerboard test were performed to mark the well-resolved areas of the map. Supplementary Figures S2, S3 in the Supplementary Material show the checkerboard tests performed on both V P (maximum amplitude 400 m/s) and V S (maximum amplitude 300 m/s) models. Areas well resolved by the checkerboard test strongly match with the DWS contour line. Hence, the tomographic figures are presented with the DWS outline to show the well resolved parts of the map.

RESULTS
Travel-time residual distributions for P-and S-wave velocity models are represented in Figure 3C and Figure 3D, respectively. The mean value of both final distributions (3D model) is closer to zero than the initial distributions (1D model), particularly that of the P-wave. The RMS is reduced by 59% after 12 iterations, from the initial value of 0.175 s to final 0.072 s (yellow stars in Figure 3B). These results demonstrate that the 3D velocity model corresponds to a significant improvement in the earthquake location quality. Figure 4 shows horizontal slices at 2, 3, 4, 5, 6, 7, 8, and 10 km depth of the V P (left panels) and V S (right panels) models. The corresponding maps of the derived V P /V S ratio, ranging from 1.48 to 2.08, whose well resolved area is contoured by the DWS corresponding to the V S model, are shown in Supplementary Figure S4. The sectors resolved by the DWS are contoured by a dark gray line. The V P values span from 3.6 to 6.9 km/s, while V S values are between 2.1 and 3.9 km/s. Each horizontal slice shows dark grey circles that represent seismic events above and below 500 m from the selected depth (labeled in the upper right corner of the maps). Due to the event and station configurations, our results are well resolved mainly in the central part of the study area. Starting from depth of 2 km, the well resolved area becomes broader down to 7 km depth and then it shrinks at larger depths. Thereby, we do not interpret anomalies located at depth greater than 10 km.
Between 2 and 5 km depth, a strongly homogeneous low V P (4.5-4.6 km/s) anomaly appears (see label a in Figure 4). This low velocity anomaly becomes narrower with increasing depth until it disappears at depth of 6 km in the SW sector of the resolved area (15.9°E, 39.85°N). At the same depths, V S maps show a low velocity anomaly (see label b in Figure 4) in the SW corner of the resolved area (15.9°E, 39.87°N). A low V P (5.2-5.4 km/s) and low V S (2.9-3.1 km/s) anomaly appears in the central sector of the map. It deepens from NW to SE between 5 and 8 km depth, following the direction of the western cloud of seismicity (label c in Figure 4). In addition, the eastern cloud of seismicity (see label d in Figure 4) is marked by low V P (5.2-5.4 km/s) and mediumto-low V S (3.0-3.3 km/s).
Between 6 and 8 km depth, the central part and eastern side of the map are characterized by a high V P anomaly (6.7-6.9 km/s, see label e in Figure 4). At 6 km depth, this anomaly intercepts the western seismogenic volume of the sequence, characterizing it almost completely. Average values of V S have been found in the same location (3.3-3.5 km/s). The western side of the western cloud of seismicity (see label f in Figure 4) is marked by a high V P (6.5-6.6 km/s) and the highest values of V S (3.7-3.9 km/s). These values of V S also appear at 8 km depth, but no V P anomalies are evident at the same depth. Figures 5, 6 show the E-W and N-S V P (left panels) and V P /V S (right panels) cross sections, respectively. These cross sections are 50 km wide and include depths between 2 km a.s.l and 13 km b.s.l. They are also projected on the map located on top of Figure 5 with red lines. The range of values of V P /V S spans from 1.45 to 2.1. All the anomalies in terms of the achieved P-and S-wave velocity models are retrieved and shown at depth also in terms of V P /V S .
The low V P anomaly found between 2 and 5 km depth results also in a low V P /V S anomaly (1.55-1.65, see label a in Figure 5 AA', BB'; Figure 6). The shape of this low V P and low V P /V S anomaly in the E-W direction adheres to the slope of the southern edge of the western fault plane that was the most active during the 2010-2014 seismic sequence. From both directions, it is quite clear that as the low V P , low V P /V S body emerges at around 5 km depth, the seismicity of the western cluster ends. This seismicity basically occurred in the deeper layer characterized by high V p (V P 6.0 km/s on average), as shown by the horizontal slices. We also find out that this high V P leads to high V P /V S values (V P /V S > 1.9, see label e in Figure 5 BB', CC' and Figure 6 FF', GG' HH').
On the other hand, the resulting V P /V S for the medium surrounding the eastern cluster of seismicity, responsible for the M L 4.3 earthquake, shows average values (1.65-1.75, see label d in Figure 5, BB'; Figure 6, HH'). Figure 6 EE' and FF' show that the low V P anomaly, deepening from 5 to 8 km depth, is also characterized by low V P /V S (1.65-1.75, see label c). The slope of this low V P anomaly follows part of the western earthquake group until about half of the cloud of hypocenters ( Figure 6, FF').
Finally, the western side of the western cluster of seismicity is marked by high V P /V S (>1.85, see label f in Figure 5 BB' and Figure 6, EE').

DISCUSSION
The tomographic images obtained from our analysis show a strong V P increase at an average depth of 5 km, as marked by a black line in Figures 5, 6 (left panels). We interpret this velocity contrast as the lithological transition from the shallower ductile mesozoic sequence of the Apennine platform (AP) to the deeper limestones of the inner Apulian platform (IAP) (Improta et al., 2000;Improta et al., 2017). This interpretation is confirmed by the structural sections inferred for the Mercure-Pollino region (Ferranti et al., 2017). The shallow AP formations include both the areas labeled a and d (Figures 5, 6). Usual V P values of the Apennine formations found in other similar geological setting of the Southern Apennines are around 5.3-5.5 km/s (Improta et al., 2000;Improta et al., 2017). These values have been found on the Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 735340 Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 735340 8 eastern side of the study area (label d in Figures 4-6) and decrease moving toward west (label a in Figures 4-6), as well as V P /V S . This might suggest a different rheologic behavior moving within the same AP. The substantial absence of seismicity within the volume characterized by anomaly a suggests that this ductile and highly fractured volume acts as a barrier, preventing the seismicity generated by the western volume to propagate upward. Combination of low V P and low-to-medium V P /V S indicates absence of significant fluid saturation within the AP crustal layer. It is noticeable, however, that a further reduction of V P /V S comes out around 0-2 km depth ( Figure 6 FF'), but further analysis is needed to fully explain the nature of this anomaly.
The central sector of the study area shows the shallower AP deepening between 5 and 8 km depth in the NNW-SSE direction within the deeper IAP (label c in Figure 5 CC', Figure 6 FF'). We interpret the peculiar shape of this interface as a portion of the AP-IAP fold and thrust interface. From the geological point of view, this result is not surprising, relying on the fact that Southern Italy is floored by the fold and thrust belt in the western and central part (Ferranti et al., 2017). More specifically, regarding the Pollino area, it has been found that during the late contraction stage (Late Pliocene to Early Pleistocene), the thick Apulia crust was involved in the thrust system and then reversely reactivated as normal faults in the actual extensional regime, responsible for the seismic hazard of the area (Ferranti et al., 2017;Filice and Seeber, 2019). For the first time, a three-dimensional tomographic image shows with impressive detail what until now had been hypothesized by geological sections. This low V P -low V P /V S extension of the AP within the IAP (label c) agrees also with the low attenuation volume observed by Sketsiou et al. (2021) between the western and eastern cluster of seismicity of the 2010-2014 Pollino sequence.
High V P and high V P /V S anomalies characterize the IAP west (label f) and east (label e) of the western seismogenic volume of the 2010-2014 Pollino seismic swarm (Figures 4-6). The western f anomaly is confined between 5 and 8 km depth and seems not to reach the western volume of seismicity. Instead, the eastern e anomaly, mostly confined between 7 and 10 km depth, intersects the southern edge of the western seismogenic volume at 6-7 km depth. This means that fluids might be likely involved in the nucleation and development of the western swarm-like sequence. Indeed, similar values of V P /V S have been associated with high pore fluid pressure and fluid migration along active normal fault systems in the central and Southern Apennines (Chiarabba et al., 2009;Amoroso et al., 2014;Baccheschi et al., 2019;Chiarabba et al., 2020). These values of V P and V P /V S provide a Poisson ratio in the range 0.3-0.4 which are marker of water saturated rocks in the investigated crustal volume (Dvorkin et al., 1999;Dutta, 2002). This resolution boundary, due to the shallow seismicity in the swarm area, does not allow us to understand whether the high V P /V S extends to larger depths. We speculate about a possible deeper extension to the SE of the swarm area, relying on the high attenuation volume on the SE side of the map at 12.5 km depth found by Sketsiou et al. (2021), perfectly matching with the eastern high V P /V S found in this work. In the same paper, another high attenuation body has been found west of the study area that can be related to the f anomaly found in this work. The SE sector of the Pollino area has also been identified by  as high scattering and high absorption area in a broad frequency band. These patterns depicted in their 2D maps have been interpreted as a complex fluid-filled connected fault network, where the largest and deepest structures are located SE of the target area (Cinti et al., 2002;Cheloni et al., 2017). Sketsiou et al. (2020) using a novel technique confirm the results in terms of absorption anomalies location. The two seismogenic volumes involved in the 2010-2014 Pollino sequence belong to two different geological and geophysical settings. The eastern seismogenic volume is entirely included in the low-to-average V P /V S AP, while the western volume is included within the high V P and high V P / V S brittle IAP, along the IAP-AP fold and thrust interface ( Figure 5, BB'; Figure 6, FF', HH'). Different geological settings featured the seismicity of these two volumes. The former sequence developed as a classic mainshock-aftershock sequence, following the M L 4.3 earthquake on May 28, 2012, while the latter developed as combination of swarm-like and mainshock-aftershock sequence (Passarelli et al., 2015). Regarding the western seismogenic volume, Napolitano et al. (2021) applied the master-slave relative location technique (Got et al., 1994) to clusters of earthquakes with local magnitude ranging between 0.6 and 2.7. Their results allowed to reconstruct with high detail the most active fault of the western seismogenic volume. In this work, we computed the absolute location of the cluster reference events using the 3D velocity model obtained from our tomographic inversion and added the hypocenters of the 423 earthquakes analyzed by Napolitano et al. (2021). Figure 7 shows the hypocenters of relative located cluster events (red circles) and the hypocenters of events used in this paper (grey circles). These hypocenters lay within the same seismogenic volume and depict a 10 km × 4 km wide fault system, located between 4 and 8 km depth, orientated NW-SE, and deepening towards SW with a dip of around 40-45°, as found by Napolitano et al. (2021). We argue that the structural barrier represented by the AP did not allow the rupture to further widen, preventing it from generating events larger than M5.0. The likely presence of pressurized fluids stopped by a structural barrier is also suggested by Pastori et al. (2021). FIGURE 7 | Hypocenter locations from three different views (panels A-C) of the western seismogenic volume achieved in this work (grey circles) and cluster events (red circles) by , achieved relocating the reference event using the new 3D velocity model.

CONCLUSION
The detailed knowledge of the crustal structure of an area considered as seismic gap, but affected in the recent past by a seismic sequence, is fundamental to evaluate and mitigate the seismic risk. The widest seismic gap in Italy, the Pollino area, recently experienced a long-lasting seismic sequence, characterized by both mainshock-aftershock and swarm behavior. Despite the fact that many studies have been recently carried out about this area, a detail of the crustal structure involved in the sequence and of the surrounding area has not been provided yet. The use of manually picked Pand S-wave arrival times of local earthquakes allowed to obtain three-dimensional images of the crustal volume involved in the 2010-2014 Pollino seismic sequence and near surrounding crust with unprecedented detail. The images, well resolved within the first 10 km depth, point out the main lithological units and lateral variations of velocity within this seismic gap area. Although the crustal structure of the Pollino area is mainly composed of two large overlapping geological units, namely, the Apennine platform and the underlying inner Apulian platform, our tomographic results highlight the complexity achieved by them in this area. Between 5 and 8 km depth in the central sector of the study area, the Apennine platform deepens in the NNW-SSE direction within the inner Apulian platform, forming what we interpreted as a fold and thrust, a common feature in the Southern Apennines framework. Our images clearly show that the two seismogenic volumes of the last Pollino sequence occurred in two different geological settings. The eastern cluster of seismicity, characterized by an aftershock sequence following the M L 4.3 earthquake, entirely developed within the low V P and low-to-average V P /V S Apennine platform. The larger western cluster occurred at the boundary of the Apennine platform-inner Apulian platform fold and thrust. The southern edge of this cluster is characterized by high V P and high V P /V S , whose values provide a Poisson's ratio indicative of water saturated rocks. In addition, cluster of similar waveforms, analyzed in a previous paper by Napolitano et al. (2021), have been re-located with the new 3D model and added to the hypocenters located in this work. These repeaters occurred along or very near the same western fault plane well depicted in this work: a 10 km × 4 km wide structure, located between 4 and 8 km depth, orientated NW-SE, and deepening towards SW with a dip of around 40-45°. The structural barrier represented by Apennine platform stopped the propagation of the seismicity toward north and upward.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
FN, AG, and ML performed the data processing and the manual picking of the arrival times. FN and OA conceived this work and worked at the tomographic inversion. FN, SG, OA, ML, and PC were involved in model analysis and interpretation of results. FN, OA, SG, and ML wrote the paper. PC provided funding for this work and coordinated the working group. All co-authors were involved in the review of the manuscript.