The 2011–2014 Pollino Seismic Swarm: Complex Fault Systems Imaged by 1D Refined Location and Shear Wave Splitting Analysis at the Apennines–Calabrian Arc Boundary

In the years between 2011 and 2014, at the edge between the Apennines collapsing chain and the subducting Calabrian arc, intense seismic swarms occurred in the Pollino mountain belt. In this key region, <2.5 mm/yr of NE-trending extension is accommodated on an intricate network of normal faults, having almost the same direction as the mountain belt. The long-lasting seismic release consisted of different swarm episodes, where the strongest event coinciding with a M L 5.0 shock occurred in October 2012. This latter comes after a M L four nucleated in May 2012 and followed by aseismic slip episodes. In this study, we present accurate relocations for ∼6,000 earthquakes and shear-wave splitting analysis for ∼22,600 event-station pairs. The seismicity distribution delineates two main clusters around the major shocks: in the north-western area, where the M L 5.0 occurred, the hypocenters are localized in a ball-shaped volume of seismicity without defining any planar distribution, whilst in the eastern area, where the M L 4.3 nucleates, the hypocenters define several faults of a complex system of thrusts and back-thrusts. This different behavior is also imaged by the anisotropic parameters results: a strong variability of fast directions is observed in the western sector, while stable orientations are visible in the eastern cluster. This tectonic system possibly formed as a positive flower structure but as of today, it accommodates stress on normal faults. The deep structure imaged by refined locations is overall consistent with the complex fault system recently mapped at the surface and with patterns of crustal anisotropy depicting fractures alignment at depth. The possible reactivation of inherited structures supports the important role of the Pollino fault as a composite wrench fault system along which, in the lower Pleistocene, the southward retreat of the ionian slab was accommodated; in this contest, the inversion of the faults kinematics indicates a probable southward shift of the slab edge. This interpretation may help to comprehend the physical mechanisms behind the seismic swarms of the region and defining the seismic hazard of the Pollino range: nowadays a region of high seismic hazard although no strong earthquakes are present in the historical record.

In the years between 2011 and 2014, at the edge between the Apennines collapsing chain and the subducting Calabrian arc, intense seismic swarms occurred in the Pollino mountain belt. In this key region, <2.5 mm/yr of NE-trending extension is accommodated on an intricate network of normal faults, having almost the same direction as the mountain belt. The long-lasting seismic release consisted of different swarm episodes, where the strongest event coinciding with a M L 5.0 shock occurred in October 2012. This latter comes after a M L four nucleated in May 2012 and followed by aseismic slip episodes. In this study, we present accurate relocations for ∼6,000 earthquakes and shear-wave splitting analysis for ∼22,600 event-station pairs. The seismicity distribution delineates two main clusters around the major shocks: in the north-western area, where the M L 5.0 occurred, the hypocenters are localized in a ball-shaped volume of seismicity without defining any planar distribution, whilst in the eastern area, where the M L 4.3 nucleates, the hypocenters define several faults of a complex system of thrusts and back-thrusts. This different behavior is also imaged by the anisotropic parameters results: a strong variability of fast directions is observed in the western sector, while stable orientations are visible in the eastern cluster. This tectonic system possibly formed as a positive flower structure but as of today, it accommodates stress on normal faults. The deep structure imaged by refined locations is overall consistent with the complex fault system recently mapped at the surface and with patterns of crustal anisotropy depicting fractures alignment at depth. The possible reactivation of inherited structures supports the important role of the Pollino fault as a composite wrench fault system along which, in the lower Pleistocene, the southward retreat of the ionian slab was accommodated; in this contest, the inversion of the faults kinematics indicates a probable southward shift of the slab edge. This interpretation may help to comprehend the physical mechanisms behind the seismic swarms of the region

INTRODUCTION
In a long period between 2011 and 2014, intense seismic swarms developed at the edge between the Calabrian Arc and the Apennines (Govoni et al., 2013;Passarelli et al., 2015).
This area is a crucial node at the rim of the ionian subduction (Totaro et al., 2014;Chiarabba et al., 2016;Palano et al., 2017). Pollino massif has been retained as the major seismic gap ( Figure 1C) of peninsular Italy . The lack of large events (M W ≥ 6.5) in historical catalogs for a 100 km-long section concurred with this hypothesis in a framework of a unique and continuous seismic belt first introduced by the pioneering intuition of Omori (1909). Therefore, seismic unrest of past years became a great concern for the management of short-term hazard.
According to the gap hypothesis, seismic catalogs did not report large earthquakes in historical time in the Pollino area, while they are present in the north and south confining regions (Tertulliani and Cucci, 2014;Rovida et al., 2016). Few seismic events with a magnitude of likely lower than six are documented, including the M W 5.6 "Mercure" earthquake in 1998 (Brozzetti et al., 2009). Although historical earthquakes are not documented, the seismic hazard map reports a 10% chance of exceeding 0.225 g in the next 50 years (Gruppo di Lavoro 2004). In fact, despite some discordance on fault geometry, paleoseismological studies found evidence for slip occurring in the past 10 kyrs and consistent with M > 6.0 earthquakes (Michetti et al., 1997;Cinti et al., 2002;. These findings suggested the Pollino (P) fault ( Figure 1B) as an active source of major earthquakes and this strong potential allowed for it to be included as a "debated seismogenic source" in the database of Individual Seismogenic Sources (DISS Working Group 2018; DOI:10.6092/INGV.IT-DISS3.2.1). Recent studies (Passarelli et al., 2015;Cheloni et al., 2017), reporting mixed seismic/ aseismic strain release for the 2012 seismicity in the Pollino gap, promote a long recurrence-time for strong magnitude earthquakes.
In this work, we produce accurate locations of seismic events that occurred during the swarm in the years between 2011 and 2014, because of the employ of a quite dense network composed of both temporary and permanent seismic stations. Later on, we discuss the main seismological features of the seismicity, focusing on its evolution in space and time, and provide refined earthquake locations to investigate the presence-and image the geometry-of the fault system in a 1D velocity model. To refine our knowledge of the geometry of the fracture field and the structure of the crust, we also compute crustal anisotropic parameters from the splitting of local S-waves.
Seismic anisotropy is manifested on a 3-component waveform through the shear wave splitting phenomenon.
Seismic shear wave energy, when it enters into an anisotropic medium, is divided into two components having orthogonal polarization directions traveling at different velocities.
Usually, the fast direction (φ) term is used to indicate the polarization direction of the fastest wave and delay time (δt) to identify the lag between the fast and slow components.
The main causes of crustal anisotropy are represented by highly foliated metamorphic rocks, layered bedding in sedimentary formations, or preferentially aligned fracture or microcracks.
Crustal anisotropy, induced by the alignment of vertical fluidfilled cracks, is often used to discriminate the state of stress in the rock volume since in many cases cracks are preferentially aligned by the local active horizontal stress field. Moreover, the δt value is proportional to the thickness of the sampled anisotropic volume and therefore the strength of the anisotropy. The above assumption well describes the theory of the Extensive Dilatancy Anisotropy (EDA) model proposed by Crampin (1993).
The Pollino area shows a wide range of geological units from sedimentary carbonatic to metamorphic rocks such as a complex structural setting. Generally, a totally sedimentary contest, characterized by horizontal layering interfaces, is explained with a vertical transverse isotropic (VTI) symmetry. Instead, the presence of near-vertical fractures aligned in a preferred orientation, as predicted by the EDA model, makes a change this symmetry into horizontal transverse isotropy (HTI) and makes the medium orthorhombic (Tsvankin, 2001;Watson and van Wijk, 2015).
S-wave splitting parameters are tightly related to the geometries and intensities of the fracture field in the area and the presence of fluids during swarm activities in the upper crust, above the earthquake clusters. We investigate the distribution in space of these parameters looking for the possible physical mechanisms behind the seismic swarm and for understanding how the prolonged seismic activity should be interpreted in evaluating the seismic hazard of the region.

GEOPHYSICAL AND GEOLOGIC OUTLINES
The study area is placed in the transition sector between the Apennines and Calabria arc (Figure 1), where one of the last fragments of the former Tethys Ocean is subducted at depth. Across the lineament, referred to as the "Pollino line" (Van Dijk et al., 2000, and references therein) the geological characteristics at the surface change rapidly from the carbonate platform units to the San Donato metamorphic core.
The subduction derives from the ionian oceanic plate sinking below the Calabrian Arc-Southern Tyrrhenian Sea and it represents a portion of the fragmented tectonic edge between Frontiers in Earth Science | www.frontiersin.org April 2021 | Volume 9 | Article 618293 two slowly converging macro-plates: Eurasia and Africa (Lucente et al., 2006, among many others, and references therein). Several seismological studies well document the lithospheric structure of the area (e.g., Piana Agostinetti et al., 2008Agostinetti et al., , 2009Di Stefano et al., 2009;Totaro et al., 2014), and the geometry of the subduction (Chiarabba et al., 2008, and references therein). The Calabrian Arc-Southern Tyrrhenian Basin system is characterized by E-W extensional tectonics, despite the N-S slow convergence between Eurasia and Africa macro-plates. Since the late Miocene, the Calabrian Arc slab was subject to a quick rollback, moving at a rate of 5-6 cm/yr from E to SE, which is greater than the ∼1-2 cm/yr rate of convergence between Europa and Africa (Faccenna et al., 2004). Despite this, in the late Pleistocene, subduction and rollback slowed down and are probably advancing at <1 cm/yr (D'Agostino and Selvaggi 2004). Seismicity and tomographic images of the slab may point out the southward lateral migration of the slab boundary, progressively following the subduction zone lateral reduction caused by the slab detachment process (Orecchio et al., 2015).
The study area structural setting (Figure 1) is indeed a consequence of the different tectonic phases. During the Middle Miocene-Late Pliocene, the Verbicaro and the Pollino Units compressional structures were dislocated by a set of WNW-ESE-oriented regional strike-slip faults characterized by left-lateral kinematics (Ghisetti and Vezzani 1982;Van Dijk et al., 2000). The presence of lateral step in these faults occasionally causes the formation of positive and negative flower structures in the region that, in certain conditions, could be reactivated.
Geodetic measurements show today a continuous extensional area that runs along the Southern Apennines ridge and reaches the Pollino region, which is subject to NE-SW extension (D'Agostino and Selvaggi, 2004); the extension rate seems to reduce from the Southern Apennines to the Calabria-Lucania boundary region (D'Agostino et al., 2013). These observations suggest that the Pollino area is affected by an increase of tectonic strain and deformation that results in a composite system of active normal faults mainly oriented near-parallel to the Apenninic chain (Brozzetti, 2011).
Recently Brozzetti et al. (2017) mapped in detail the active faults at the Calabria-Lucania border ( Figure 1A) finding both NW-SE trending SW dipping and N-S striking E dipping normal fault sets. This complex fault system is located between two major normal tectonic structures, already known in this area as the "Rimendiello-Mormanno" Mercure Basin (RM-MB) fault system to the Northwest and the Castrovillari (CV) fault to the South ( Figure 1B). These are cataloged as individual seismic sources in the Italian database of the Individual Seismogenic Sources (DISS Working Group 2018). and dipping at about 45°. On the whole, focal mechanisms are coherent with a NE-SW extension observed in the area (Totaro et al., 2013).
The seismicity mainly occurred near the Mormanno village (western larger dashed ellipse in Figure 2). On May 28th, 2012, at about 5 km eastward of the first main cluster of seismicity a shallow M L 4.3 event occurred ( Figure 2). The seismic activity remained concentrated in this area (eastern smaller dashed ellipse in Figure 2) until early August. Thereafter, seismicity moved again westward to the Mormanno activated area, where several events with a magnitude larger than 3.0 preceded the strongest earthquake (M L 5.0) that occurred on October 25th, 2012. The seismicity rate persisted high for some months, but magnitudes did not exceed M L 3.7. At the beginning of 2013, the seismic activity started to decrease (inset in Figure 3), remaining relatively low until June 2014, when a magnitude M L four Frontiers in Earth Science | www.frontiersin.org April 2021 | Volume 9 | Article 618293 earthquake occurred in the eastern cluster, giving rise to a sudden increase of the seismic rate in the surrounding area. The seismicity rate kept decreasing during 2015, although remaining above the background seismicity threshold (inset in Figure 3). There is evidence that this seismic sequence behaves as a swarm more than as a mainshock-aftershock seismic sequence, in fact, this preliminary hypothesis is supported by a more accurate and detailed analysis performed by Passarelli et al. (2015); they found that 75% of earthquakes may be related to a transient forcing and only 25% could be interpreted as aftershocks, We show in Figure 3 the temporal evolution of the seismic sequence in term of seismicity rate to investigate different behaviors of the two clusters. The histograms of the daily number of events and the cumulative number of events for the whole area (green line in FIGURE 3 | Temporal evolution of the swarm: black histogram defines the daily number of events, the green line is the cumulative number of events in the whole region, the blue line and the red line are the cumulative number of events for the western and eastern cluster, respectively. Only the blue line trend seems to follow a mainshock-aftershock behavior.
FIGURE 4 | Map of temporary and permanent seismic networks deployed in the area during the seismic activity. In the upper insets an example of installation is reported; in the lower panel is shown the GFZ seismic array configuration. In the map, the two main-shocks are also reported as white stars.
Frontiers in Earth Science | www.frontiersin.org April 2021 | Volume 9 | Article 618293 5 Figure 3) show an increase of the seismicity months before the mainshocks, with a seismic rate that increases and decreases at varying pace.
We also believe, when considered the two main clusters separately, that the events clustered inside the western area have swarm behavior (blue line in Figure 3); whilst the events occurred in the eastern region seem to follow a mainshockaftershock trend (red line in Figure 3).

Temporary Seismic Network
During these four years of intense seismic activity, a group of temporary INGV seismic stations (INGV Seismological Data Center, 2006) and Deutsches Geo Forschungs Zentrum (GFZ) seismic stations and array (Passarelli et al., 2012) were installed, for the purpose of enhancing the earthquakes identification capability of the Italian National Seismic Network there in real-time ( Figure 4, see also Govoni et al., 2013 for details) and to refine the accuracy of the locations using also off-line data. The temporary seismic network was composed of a variable number of seismic stations during the entire seismic sequence, from eight stations in late 2011 to about 20 stations in late 2012, in order to refine the hypocentral locations of small-magnitude earthquakes and ensure a lower magnitude detection threshold of the network (see also De Gori et al., 2014, for the timing of stations deployment). In May 2015 the temporary network deployed in the Pollino region was removed keeping only two stations to maintain a lower magnitude detection threshold in the area with respect to the regional standard.

REFINED EARTHQUAKES LOCATION
The analysis of the seismic sequence greatly benefited from data provided by the temporary seismic stations installed soon after the sequence started (Figure 4), many of which contributed to the real-time monitoring (Govoni et al., 2013). Earthquakes were originally located by analysts on seismic surveillance duty with the HYPOINVERSE-2000 code (Klein, 2002) using a 1D velocity model with two layers at 0-11 km (Vp 5.0 km/s) and 11-38 km (Vp 6.5 km/s) depth over a half-space at depth >38 km with Vp 8.0 km/s. These locations (Figure 2), available in real time on the ISIDe website, were used as the initial dataset.
For this study, we carefully reprocessed the earthquakes reading P-and S-phases arrivals on the three-component digital recordings at local stations from October 2011 to March 2013. As the first step, we located all the earthquakes using the program Hypoellipse (Lahr, 1989) and a starting Vp model taken from Chiarabba et al., 2005, adopted for the whole Italian region.
The average Vp/Vs in the crustal volume where seismicity developed has been estimated through the modified Wadati method (Chatelain, 1978). Since seismicity mainly occurs at depths of 5-10 km, most of the observations consist of upgoing ray-paths, thus the retrieved Vp/Vs (1.81, see Figure 5A) is strongly representative of the focal volume.

Locations With Minimum 1D Velocity Model
In order to improve earthquake locations, we computed an ad hoc 1D velocity model by applying the inversion scheme introduced by Kissling et al. (1994) and implemented in the code VELEST. To obtain a stable 1D reference model, we selected a subset of events that meet the following selection criteria: rms <0.5, azimuthal gap <200°, location errors within 2 km, and at least seven P-and 3 S-wave readings.
We then inverted P-and S-wave arrival times, from a subset of 808 selected earthquakes, to simultaneously compute hypocenter solutions and the minimum 1D velocity model parameterized by horizontal layers with Vp values varying with depth. The damping parameter was chosen by running different inversions as the best trade-off between data variance reduction and model complexity. The optimum model ( Figure 5B) was achieved after 10 iterations, obtaining a final rms of 0.12s with a variance reduction of 48%.
Below a thin layer of sediments cover, with a thickness of 2 km and Vp of 4 km/s, the P-wave velocity ( Figure 5B) displays values characteristic for consolidated sedimentary rocks like Mesozoic evaporites and carbonates, in agreement with the velocity values found by Improta et al. (2017), in an adjacent area.
The retrieved relatively high Vp/Vs values (1.8) is often observed in Mesozoic carbonate rocks of the Apennines and can be related to fracturing and fluid saturation (e.g., Improta et al., 2014).
Finally, we use this minimum 1D model to relocate the whole set of earthquakes through the Hypoellipse code (Lahr, 1989). We obtain a final set of 2,323 well located earthquakes; the performance of the adopted approach, in terms of rms and location errors can be evaluated in Figure 6.
Results of the analysis described above are shown in Figure 7, where map and cross section views of the seismic sequence are provided. The relocated earthquake catalog represents a sensible improvement with respect to that shown in Figure 2 in terms of Wadati (1933) diagram computed on all the earthquakes with up-going ray-paths to the stations. For each event, DTp and DTs are the differences between P-and S-phase arrival times, respectively, at couples of stations. The Vp/Vs value is the slope of the straight red line, while R is the correlation coefficient of the linear regression (B) Starting (dark gray) and final (red) P-wave 1-D velocity model computed using the VELEST code (Kissling et al., 1994).
Frontiers in Earth Science | www.frontiersin.org April 2021 | Volume 9 | Article 618293 geometrical definition of the activated fault system, whose characteristics are here recognizable. The analyzed earthquakes describe two main areas of seismicity: Western and Eastern clusters; in the first one, the M L 5.0 is located instead in the second one the M L 4.3 occurred. The cross-sections well delineate two different behaviors: in the northern cross-section ( Figure 7B) the hypocenters define a ballshaped volume without imaging any planar lineament. The distribution of the M L >3.0 best-recorded events in the area (yellow stars in Figure 7B), also agrees with more diffuse seismicity rather than depicting a planar one.
On the contrary, in the southern cross-section ( Figure 7C), hypocenters define at least three dipping fault planes. Those planes include the one potentially responsible for the M L 5.0 event, which was imaged by the seismicity hypocenters also before the occurrence of such mainshock. This plane has been reconstructed with a similar geometry also by other recent studies (Brozzetti et al., 2017;Napolitano et al., 2021). But with respect to these latter works, thanks to the ∼2,300 events analyzed and relocated, we were also able to depict an antithetic E-dipping fault plane located to the West of the main fault plane, and another NE-dipping, on the East, related to the M L 4.3 event, possibly the source of this second mainshock.

Relative Locations With hypoDD
To further improve locations, we applied the double-difference (DD) relocation algorithm using the hypoDD code (Waldhauser and Ellsworth, 2000;Waldhauser, 2001). The algorithm FIGURE 6 | Histograms summarizing the main outputs of the localization procedure, in terms of frequency of observations for each of the following parameters (from top left to bottom right): number of observations for P phases; number of observations for S phases; residual rms (s); horizontal location errors (km); vertical location errors (km); azimuthal gap between stations (°).
Frontiers in Earth Science | www.frontiersin.org April 2021 | Volume 9 | Article 618293 minimizes the residuals between observed and calculated traveltime differences for pairs of earthquakes at common stations by iteratively adjusting the difference between the hypocenters. This condition holds if the hypocentral distance between earthquakes is small compared to the source receiver separation. We used P and S wave arrival times picked from three component seismograms, the starting locations derived above and the minimum 1D velocity model already described. Doubledifferences are generated to link each event to several neighbors, so that all the events are connected and the adjustment to hypocenter parameters is simultaneously determined. In our case, all the selected events are strongly linked to the others forming one main cluster. We performed 10 iterations optimizing the damping parameter to obtain good conditioning and convergence of the solution. The iterations were grouped into three sets for better control of the inversion procedure and to achieve a consistent convergence. The final adjustments are in the order of meters.
The DD locations show a seismicity distribution similar to that obtained with the minimum 1D model, but emphasizes some aspects of the geometry like the opposite dipping of the two main faults (Figure 8). The northern part of the main cluster still does not present a clear planar geometry.

SHEAR-WAVE SPLITTING
From ∼6,000 events recorded at temporary (T07*) and permanent stations deployed in during the seismic sequence, we analyzed 22,623 event-station pairs looking for crustal anisotropy to evaluate the shear wave splitting parameters: fast direction polarization (φ) and delay time (δt). Seismicity used in the anisotropic analysis is mainly concentrated under Mormanno village between 5 and 10 km depth, with magnitudes ranging from 0.2 to 3.0. To avoid complications on the S-wave arrival, due to the P-wave coda of a strong events that could make difficult the splitting analysis, we decided to discard events with a magnitude greater than 3.0. We used Anisomat+ (Piccinini et al., 2013) a set of MatLab scripts able to retrieve automatically φ and δt from the seismic recording of local earthquakes. We tuned this code by (1) Corner frequencies of cf1 1 Hz and cf2 12 Hz of the four poles 2 pass Butterworth filter; (2) PRE variable set on 0.15 s time before the S-picking in which the analysis window starts; (3) DUR variable, representing the time after the S-picking, is computed using the formula DUR C 1/cf2 where C is a coefficient, in this case set on 2.5 to ensure including at least a complete filtered shear wave cycle; (4) the total duration of the analyzed signal is 0.35 s (PRE + DUR) and encompass the S-wave arrival.
To ensure that the observed shear wave energy pertains only to the horizontal plane, the code performs a number of tests. In detail, the code evaluates i) the S-to-P ratio, the ratio between the amplitude of the horizontal components sum (RMSs) and the amplitude of the vertical component (RMSp) computed on the analysis window as RMSs/RMSp and it is set to be greater than four; ii) the geometrical incidence angle computed as ic sin-1(Vs/Vp) is set to be lower than 45° (Booth and Crampin, 1985), where Vs and Vp are the velocities of S and P waves, respectively.
Following Bowman and Ando (1987) Anisomat + calculate the cross-correlation between the two horizontal components. This allows measuring the similarity of the pulse shape between the S-waves and their delay time (δt). These two time series in fact have similar shapes, mutually orthogonal oscillation directions, and travel with different velocities.
Horizontal components of the ground motion are then rotated (with steps of 1°from 0 to 180) and for each step, the crosscorrelation function is stored. Finally, a two-dimensional matrix is obtained and after a grid search, the code estimates the couple of lag time and rotation angle which maximizes the crosscorrelation coefficient.
These two quantities represent the delay time between the slow and fast shear wave (δt) and the polarization azimuth (φ) of the incoming fast shear wave. Associated errors are defined by a criterion that reflects the shape of the cross-correlation matrix and how rapidly the matrix grows around the observed maximum value. In practice, the goodness of the estimation is related to the difference in delay time and in degrees between the cross-correlation maximum and 95% of the parameter values themselves. The mean errors for fast directions and delay times are 14.7°and 0.008 s, respectively (see Table 1 for mean values at each station). At the end the code returns an ASCII output file composed by 21 fields described in the Supplementary Material. For further details, the reader should refer to Piccinini et al., 2013. Seismicity used in the anisotropic analysis is mainly concentrated under Mormanno village between 5 and 10 km depth, with magnitudes ranging from 0.2 to 3.0. To avoid complications on the S-wave arrival, due to the P-wave coda of a strong event that could make difficult the splitting analysis, we decided to discard events with a magnitude greater than 3.0. After the Anisomat + analysis and its quality controls, in order to discuss a high-quality anisotropic dataset, we admitted only those results with a cross-correlation coefficient CC ≥ 0.75, which allowed us to obtain a total of 3,727 fast single measurements (δt > 0.02 s) and 2,040 null results (δt ≤ 0.02 s). A measure is defined null when the S wave does not show splitting; this happens, in an anisotropic medium, when the initial S-wave polarization is coincident to the slow or the fast axis, and furthermore, nulls do not constrain the delay time.
We present the fast and null directions as interpolated values on a regular grid-box which nodes are 2-km spaced (see Supplementary Figure S1 for details) evaluated with the Tomography Estimation of Shear wave splitting and Spatial Averaging (TESSA) program (Johnson et al., 2011). In each block, the rose diagram of fast values and the associate mean fast direction are displayed, providing i) the standard deviation is <30°, ii) the standard error is <10°, and iii) at least 10 rays pass through the grid-box. Furthermore, TESSA weights the data according to the fast errors estimated by means of Anisomat+ and to the weighting regime selected, in this case, we chose a 1/d regime, where d is the length in km between the grid-box and the station. The insets of Figure 9, the orange and green rose diagrams are frequency plots representing splitting fast and null measurements, respectively, for the whole area. Prevalent NW-SE fast directions are observed: the average fast direction is N124°(see Table 1 for mean values at each station) while the null directions show two main peaks one strikes N127°(considered as the fast direction) and the other roughly orthogonal to it (considered as the slow direction).
In the crustal volume sampled by S-waves, the different mean fast directions observed lead us to suppose the presence of intense and complex anisotropic anomalies. The prevalent φ is both parallel to the major fault structures and consistent with the local stress field being parallel to the NW-SE S Hmax , so the ambiguity on the probable crustal anisotropic source is not solved: structural anisotropy or aligned microfractures opened by the active stress field?
Fast orientations delineate three main areas North-East, West and South-East: 1) NNW-SSE mean direction (stable and coherent for fast and null measurements); 2) WNW-ESE mean direction (variable and not coincident for fast and null measurements) and 3) from NW-SE to E-W mean direction (stable and coherent for fast and null measurements), respectively.
In the north area of the activated fault system, the relocated hypocenters define a diffuse ball-shaped volume of seismicity instead of fault planes. This behavior could have influenced the anisotropic parameters, in fact, the strong variability of fast directions, noticed in the western area, can indicate the absence of dominant structures and the existence of a complex diffusivity system where fluids can flow. Instead, in the eastern part where the hypocenters delineate a main SW-dipping fault plane with an antithetic plane to the west and another dipping to NE on the east side ( Figures 7C, 8C), we observed stable and coherent φ with respect to the major structures, in this case, we can assume that anisotropy is structurally controlled as predicted by Zinke and Zoback (2000) structure-related model. Delay time singular values oscillate between 0.024 up to 0.248 s ( Figure 10A) with a total average value of 0.06 s (see Table 1 for mean values at each station). These values agree with those found in other regions of Italy along the Apennine Piccinini et al., 2006;Pastori et al., 2009;2019), or in the crust of other regions in the world, as on the Karadere-Düzce fault (Peng and Ben-Zion, 2004). We will interpret the results in terms of delay time normalized for the hypocentral distance (δtn) because we suppose that the anisotropy measured at the station is just the sum of the anisotropic structures crossed by the seismic ray (e.g., Crampin, 1991;Zhang et al., 2007).
Normalized delay time values are represented in Figure 10B as the geographical distribution of interpolated single measurements on a 2-km spaced grid, similar to that used for the fast measurements averaging. The higher culmination of δtn (from 0.007 to 0.01 km/s) is visible in the Northeast sector that represents the hangingwall of the fault that generated the M 5.0 event. This value allows us to confirm that the principal source of crustal anisotropy, in this area, is represented by the structural control, as already seen by the φ pattern.
The lower values of δtn (from 0.003 to ∼0.005 s/km) are visible in the external west and southeast areas. These values allow us to suppose the presence of a complex diffusivity fluid system where fluids can migrate along a relatively high-permeability hydraulic path but cannot be gathered. The existence of this network of fractures and/or porous fluid-filled is also invoked, i. e, by Della Vedova et al. (2001) to interpret the unusual low geothermal gradient of the region, by Barberi et al. (2004) and Totaro et al. (2014) to explain an high Vp/Vs value observed in the shallowest crust and by Napolitano et al. (2020) to justify the scatter and absorption anomalies.

DISCUSSION
The Pollino gap area is located at the edge of the southern Apennines where a broad strike-slip shear zone separates two lithospheric blocks (the Apennines and the Calabrian arc) with different evolution in the Plio-Pleistocene (Spina et al., 2009;Orecchio et al., 2015;Palano et al., 2017). This lithospheric discontinuity favored the retreat of the ionian slab during the lower Pleistocene and consists of a set of N120-trending faults, which comprises the Pollino fault and the Amendolara fault in the ionian off-shore (Ferranti et al., 2014).
In Late Pliocene-Quaternary the ancient thrust structures are crosscutting by extensional faults that forming intermountain basins along the Apennines axis (Ferranti et al., 2017), as well as FIGURE 9 | Map of the spatial average for 3,727 fast (orange bars) and 2040 null (green bars) measurements. Three sectors are recognized: West) even if a WNW-ESE mean direction is observed, fast and null orientations are variables and not coincident; North-East) a NNW-SSE mean direction is visible, fast and null orientations are stables and coherent; South-East) a trend from NW-SE to E-W mean direction is delineated from both fast and nulls measurements. The upper right insets show the total rose diagrams for both fast and null measurements.
Frontiers in Earth Science | www.frontiersin.org April 2021 | Volume 9 | Article 618293 the Mercure and the Crati basins respectively to the North and South of the Pollino fault. Common tectonic models are based on the continuity of deformation between the southern Apennines and Calabrian arc (i.e., Valensise and Pantosti, 2001); lateral heterogeneities across the Pollino mountain range instead suggest that the strikeslip shear zone permits a kinematic decoupling between the extending Apennines and that retreating Calabrian block (Chiarabba et al., 2016), as also signaled by the imprint of toroidal flow at the ionian slab edge on GPS data .
The Quaternary normal active fault systems segmented by the tectonic structures inherited from the past and the complex interaction between the activated fault planes evidenced by the refined 1D and DD relocations can be considered as the cause of the great variability defined by the anisotropic pattern. , red dashed lines with arrows represent the regional scale strike-slip episodes related to geodynamic processes, yellow stars identify the strongest events, and the main fault systems are also reported.
Frontiers in Earth Science | www.frontiersin.org April 2021 | Volume 9 | Article 618293 Recent studies tried to model the geometry of the activated fault system in this area. Napolitano et al. (2021) used a smaller portion of events (422 in 27 clusters) with respect to our dataset (∼2,300 events). Their locations were concentrated mainly on the western cluster and depicted a structure, having a similar geometry to the one that we consider to be probably responsible for the M L 5.0. On the contrary, they did not identify the antithetic fault, that we assume likely causing the M L 4.3, because their data did not sample the eastern cluster. Brozzetti et al. (2017) used both seismological and geological data to infer the geometry of activated faults. Their solutions agree with the geometry of the M L 5.0 reconstructed by our data, on the other hand, to define the fault causative of the M L 4.3, their locations showed an almost vertical alignment of earthquakes and considering the geological surface expressions they preferred a fault dip to the NE.
In our study, we compare refined 1D and DD locations and anisotropic results to better constrain the geometries of the activated fault system. The S-waves fast polarizations show a general rotation from roughly WNW-ESE (sub-parallel to CaF) to NNW-SSE (sub-parallel to the ROCS active extensional systems) to WNW-ESE (sub-parallel to the Pollino line) (orange bars in Figure 11A). Such orientations are also observed in null results, except for the Western area (green bars in Figure 11A); as previously explained we observe a null direction if the initial S-wave polarization is coincident to the source structures of the anisotropy. After these considerations, we can affirm that for the western sector, fast direction polarizations are probably more related to aligned microfractures opened by the active stress field, while nulls reflect major structure orientations confirming a structural control of anisotropy.
All these orientations well-delineate the regional scale structural complexities, such as the strike-slip episodes related to geodynamic processes and those at the small-scale represented by the local stress variations.
In general, the interaction between the irregular geometries of the fault systems could favor multisegment ruptures and/or the reactivation of inherited structures with the inversion of the kinematics due to the present stress regime. These segments take part in the generation of a complex diffusivity system where fluids can be channelized and carried toward heavily fractured rock volume where in some cases they can be gathered.
Overall, we found in the western sector, the volume where most of the seismicity occurs, a higher value of δtn. This crustal volume represents, under a structural point of view, the hangingwall of the main fault where the M 5 is supposed to be enucleated. Similarly, the 1D and DD accurate relocations depict a SW-dipping main fault plane but also an antithetic plane to the west. This area can be imaged as an ancient positive flower structure reactivated in the current extensional regime, as also considered by Totaro et al. (2015) to explained their focal mechanism solutions and high-precision locations using the DD method in the Mercure Basin.
The main faults, that border the flower structure, could play the role of a preferential pathway in which fluids are channelized and accumulated on the top of the heavily fractured rock volume and laterally confined by these structures (Figures 11B,C).
Instead, in the southeast sector, the 1D and DD relocation possibly represents a NE-dipping plane where the fewest seismicity occurs. In this area, the anisotropic pattern of fast directions recognizes a dominant WNW-WSW trend and lower value of δtn. These findings suggest that fluids migrate preferentially along hydraulic pathways having a high percentage of permeability, such as the hypothesized secondary fault, up to the free surface not allowing a constant pervading the crustal volume by fluids.

Seismogenesis
In the hangingwall of the activated fault systems (western sector), our anisotropic evidence and the excess rate of seismicity during the sequence are indicative of countless fluid-saturated cracks. Napolitano et al. (2020) through a study based on seismic scatter and absorption confirmed the presence of a highly fluid-filled micro-fractured seismogenic volume around the main activated structures during the swarm. Passarelli et al. (2015) established that the temporal and spatial behavior of seismicity in the western sector follows prevalently a swarm-like consistent with a transient, external forcing. Even Cheloni et al. (2017) evidenced a transient slip, around the mainshock, approximately 4 months before and gone on about one year.
Transient slip associated with a high percentage of fluid-filled fractures could have reduced the effective normal stress facilitating fault slip, partially inhibiting the development of large rupture (Cheloni et al., 2017;Yoshida and Hasegawa, 2018), increasing the rate of seismicity but limiting the maximum magnitude of events.
Furthermore, also lithology might affect the permanence of fluids, in fact, the carbonate rocks forming the hangingwall have an important function in both preserving and hosting the microfracture system responsible for the observed anisotropic pattern. This permanent presence of fluid could reduce the long-term fault strength, promoting aseismic fault slip and causing important weakening along the existing fault planes. In this contest, looking at the structures obtained by the 1D and DD relocations we can suggest the reactivation of an ancient positive flower structure in the actual stress regime.
We, therefore, explain the space-time-magnitude pattern of the seismic sequence, and the excess rate of seismicity, as governed by the preservation of the filled-fluid micro-fracture network in the seismogenic volume that followed transient slip.

Constraints to Hazard
Persistent swarms in the Pollino seismic gap area created repeated concerns in public and seismic risk stakeholders, for uncertainties on possible forecasting scenarios. Most of the concerns derived from the perception that the area was seismically silent for a long time and possibly prone to a large event. The analysis of the 2011-2014 swarm yields therefore precious information that might help to constrain hazard models.
The seismic gap idea, first expressed by Omori (1908), was motivated by the lack of large earthquakes since at least 10 Kyr . In the continuous extensional fault system model (Valensise and Pantosti, 2001), a 60 km long gap encloses the Pollino mountain range and partially the upper Crati valley. The faults studied by paleo-seismological trenches in the southern border of the Pollino area evidenced large slip earthquakes older than 1 Kyr (Galli et al., 2008;Cinti et al., 2002;Michetti et al., 1997), reinforcing the idea that large earthquakes are expectable. The absence of large earthquakes in the historical catalog has been associated with a combination of low strain rate and mixed seismic/aseismic strain release that could increase the recurrence time of large magnitude events (Cheloni et al., 2017).
Under these circumstances, the persistence of a pervasive fluid system network in the pop-up structure accounted for increasing the stability of the faults, partially inhibiting the generation of large earthquakes. In any case, the reactivation of pre-existent thrusts and back-thrusts of the pop-up would promote the segmentation of the seismic gap in fragments with a maximum length of 10-12 km.
Although shear wave splitting analysis does not measure directly the activity of a fault, our results might help us to delineate the seismotectonic settings of the Pollino region.
In the western sector, where the highly fractured rock and the excess rate of seismicity are observed, we can suppose a reduction of shear-strength along CRFS and ROCS faults ( Figures 11B,C), probably representing the boundaries of the flower structure. This reduction in turn can generate earthquakes also for a low level of stress accumulation.
In the eastern sector, instead, the lower values of δtn and fewer seismicity rates are concentrated; this area coincides with the upward continuation of the surface expression of CAS (Cinti et al., 2015). The absence of a strong pervasive fluid system network that reduces the effective normal stress facilitating fault slip, around the intersection between CAS fault might suggest that this eastern segment could be more prone to an unstable deformation and a large rupture.

Differences in Crustal Response to Slab Tearing Between Apennines and Calabria
Comparing Figures 7B,C, 8B,C strikingly illustrates the different behavior in the release of seismic energy between the Northern half and the Southern half of the study area. While a more standard seismicity clustered along defined planes is present in the Southern half ( Figures 7C, 8C), diffuse seismicity is found in the North (Figures 7B, 8B). Those two phenomena have been observed worldwide, but rarely in such close proximity. We speculate that such rapid spatial change in crustal behavior is caused by a different response of the Southern Apennines and the Calabrian block to the tearing of the ionian slab. Slab tearing is likely to generate diffuse deformation in the overriding plate (i.e., Maesano et al., 2020;Chiarabba and Palano, 2017;Gutscher et al., 2016). The Apennines slab has been torn in several places (Rosenbaum et al., 2008) and the tearing process has shown to affect the entire crustal volume . Such "crustal response" is likely to interest a wide area rather than a focalized rock volume on top of the tear fault (Rosenbaum et al., 2008), depending on the local rheological properties of the crust (Boncio et al., 2007;Polonia et al., 2011). Here the tearing of the ionian slab, revealed at depth by passive seismic (Rosenbaum et al., 2008;Orecchio et al., 2015), probably interested in different ways the juxtaposed two crustal domains, focusing a higher level of deformation on the Southern Apennines side, and, thus, weakening the entire crustal volume. This interpretation is supported by the distribution of earthquakes in a diffuse area (NW of Mormanno) indicating the inability of cumulating stress on relevant tectonic structures (i.e., km-scale planar faults) and corroborated by the misoriented fast polarization directions with respect to the major structures, given by the complex fluid-filled fracture network. Opposite observations in the Southern portion (planar faults, where main M L 5.0 and M L 4.3 events occurred) where fast directions reflect the strike of the major structures testifying to the different behavior of the crustal response to slab tearing, on the two sides of the boundary between Apennines and Calabria.

CONCLUSION
The seismic swarms occurred within the gap area, along the Pollino shear zone. Relocated earthquakes reveal that seismicity developed on two main NNW-trending and opposite dipping clusters with primary events being the M L 4.3 and the M L 5.0 on 28th May and 25th October 2012, respectively.
The geometry drawn by the refined 1D and DD locations delineates the activated fault system having the following characteristics: a) it is composed at least of three fault segments with a length <10 km; b) it is located in the shallower 10 km of the crust; c) it is oriented roughly NW-SE. The dip of the main fault segment, where most of the seismicity occurs such as the M L 5.0 earthquake, is toward SW while, the two other minor segment faults dip toward NE. Although with a lower resolution also Totaro et al. (2014) pointed out a similar geometry of the main fault segment.
According to the data reported in DISS database, we could place the activated fault system between the two major seismogenic sources present in the area: the Rimendiello-Mormanno Mercure Basin fault system (RM-MB) and the Castrovillari Fault (CV). The geometry of the RM-MB, constrained by geological data, is about comparable with the trend of the recent seismicity, but its dip is to the NE. Instead, CV source has a similar orientation (N20W) and a western dip, as reported by the analyzed seismicity, but its northern border match to the southernmost epicentres. Finally, our results in terms of fault orientations delineated at depth by 1D and DD refined locations could be connected to the mapped fault systems reported by Brozzetti et al. (2017).
The reactivation of an ancient positive flower structure in the actual stress regime is supported by the fast direction polarization rotation from NNW-SSE in the middle of the extensional basin toward WNW-ESE in the area where the Pollino fault is supposed to be. As imaged by location earthquakes this structure is composed of a set of thrust and back-thrust strikes mainly NW-SE.
The interaction at a depth between these faults could generate a pathway for fluids and in certain conditions, i.e. in the presence of a structural barrier, reach the overpressure status. As observed from anisotropic results, the presence and mobility of channelized fluids in the crust promote multisegment ruptures and/or the reactivation of inherited structures only if favorable oriented with the present-day stress regime.
Back thrust reactivation contributes to the segmentation of the 60 km long seismic gap in fragments with a maximum length of 10-12 km. In this vision, the subdivision could partially restrict the generation of large earthquakes in this area.

DATA AVAILABILITY STATEMENT
The data of used stations are from the Italian National Seismic Network (doi.org/10.13127/SD/X0FXnH7QfY). The waveform data of earthquakes used are available at EIDA database (European Integrated Data Archive http://eida.rm.ingv. it/).

AUTHOR CONTRIBUTIONS
MP wrote the manuscript, analyzed the data and evaluated the anisotropic parameters; LM coordinated the working group; PDG, AM and DL made picking analysis and 1D locations; AG installed the temporary stations in the field, made data available and helped in the discussions; FL, NA, DP and CC made picking analysis and discussed the results; MM coordinated the fieldwork and made picking analysis; PDL made picking analysis and evaluated anisotropic parameters; RDG, MA and AN made picking analysis; LP discussed the results.

FUNDING
This study has benefited from funding provided by the Italian Presidenza del Consiglio dei Ministri -Dipartimento della Protezione Civile (DPC). This paper does not necessarily represent DPC official opinion and policies.