A High‐Resolution Seismic Catalog for the Southern Apennines (Italy) Built Through Template‐Matching

The incompleteness of earthquake catalogs is a well‐known issue caused by our technical limitation in detecting the small‐to very small‐magnitude seismicity falling near or below the background seismic noise. The detection of small‐magnitude events is fundamental for improving our knowledge of geometry and kinematics of seismogenic sources and the spatio‐temporal characteristics of seismicity, thus leading to better models for seismic hazard. Template‐matching (TM) is a powerful technique that, based on similarity measure (cross‐correlation) of seismic waveforms, allows to detect hidden earthquakes that are similar to known events (called templates). The high computational effort often limits such technique to small areas and for short time frames (less than 1 year). In this work, we present the first application of template‐matching at regional scale for the Italian Peninsula, focusing on the Southern Apennines. We use about 3,600 high‐quality events as templates, scanning 6‐year long continuous recordings (2009–2014), at more than 180 stations of the INGV network. About 20,000 new events are found, showing a comparable quality to the template catalog in terms of hypocentral solution, reaching a decrease of the magnitude of completeness of about one unit. To highlight the improved quality of the TM catalog, we report two main examples regarding the Sannio‐Matese area, where TM allowed us to unravel relevant details on the spatio‐temporal distribution of the local seismicity.


Introduction
Despite the theoretical and instrumental advances of the last decades in the field of earthquake seismology, the initiation process of earthquakes remains poorly understood.In general, small-magnitude events may precede a mainshock, then followed by a series of aftershocks.In other cases, the seismicity shows a swarm-like behavior, with the earthquake rate gradually building up and then decreasing toward the long-term background activity.Depending on the case, the earthquake nucleation and triggering are explained with a variety of processes (or their combination) such as static stress change due to fault slip (Stein, 1999) and dynamic stress change (Felzer & Brodsky, 2006).Fluids can also play a role, causing an increase in hydrostatic pressure (Wei et al., 2015) or aseismic slip due to fluid migration (De Barros et al., 2020).However, our comprehension of the earthquake mechanism is still inadequate and mainly hampered by the technical limitations in the detection of small-to very small-magnitude seismicity.This constitutes the vast majority of the whole seismic activity, as predicted by the Gutenberg-Richter law for the frequency-magnitude distribution (Gutenberg & Richter, 1944).However, such small-magnitude events fall below the detectability threshold of conventional short-term/long-term average algorithms (Allen, 1982) used for real-time earthquake monitoring, and therefore all catalogs at our disposal are inherently incomplete.An effective method that is employed for the detection of hidden seismicity and the improvement of catalog completeness is template-matching (TM).Based on cross-correlation measures, the algorithm uses the recorded waveform of known earthquakes (templates) to identify highly similar, previously undetected events near or below the usual detectability limit (Frank et al., 2014;Gibbons & Ringdal, 2006;Shelly et al., 2007).Template-matching has been widely used across different scales and for a variety of applications, ranging from the characterization of fluid-induced seismicity (Shelly et al., 2013(Shelly et al., , 2016)), detection of repeating earthquakes (Chamberlain, Boese, & Townend, 2017;Chamberlain et al., 2014;Kato et al., 2012;Shelly et al., 2007), and for catalog enhancement from local (Cabrera et al., 2022;Essing & Poli, 2022;Vuan et al., 2018) to regional scale (Ross et al., 2019).It allows us to sharpen our image of the three-dimensional structure of seismogenic faults at depth (Ross et al., 2017;Shearer, 2002), contributing to better seismo-tectonic models and improved estimation of the overall seismic hazard in a certain area.However, large-scale applications of template-matching are still hindered by the high computational burden, mainly due to the massive I/O operation on the continuous recording.Attempts have been made to speed up the computational operation, including the optimization of the cross-correlation algorithms (Beaucé et al., 2018) and/or the parallelization on CPU and GPU infrastructures (Mu et al., 2017;Ross et al., 2019).Moreover, with the recent advent of artificial intelligence in several scientific fields, machine-learning algorithms have been developed for phase detection and association (Mousavi et al., 2020;Shi et al., 2022;W. Zhu & Beroza, 2019), and successfully employed for creating highresolution catalogs at sub-regional scale (Liu et al., 2020;Tan et al., 2021;L. Zhu et al., 2019).This work aims to produce a new, enriched earthquake catalog at regional scale, covering the entire Southern Apennines chain (Italy), an area that has been struck by several large earthquakes (M > 6-7) in the past 300 years.About 4,000 high-quality events between 2009 and 2014 were used as templates for scanning the continuous recordings of the Istituto Nazionale di Geofisica e Vulcanologia (INGV) seismic network in the same period.Computations are parallelized on a CPU cluster leading to the detection of about 20,000 new high-quality events, more than a 5fold increase with respect to the initial template catalog.In the following sections, we discuss the characteristics of the starting catalog, the employed template-matching methods, and its application.Then we discuss the characteristics of the retrieved seismicity in terms of the quality of the hypocentral solutions and its spatio-temporal distribution.Finally, we provide two main examples of application, the first regarding the 2013-2014 Matese sequence, and the second focusing on the seismicity distribution in the Sannio area.

Building the Template Catalog
A high-quality template catalog is a necessary prerequisite for the successful employment of the templatematching technique.In fact, the cross-correlation function highly fluctuates with time lag, therefore only precise P and S picks of the template's arrivals maximize the chances of detecting similar and coherent signals that can be indicative of a new, previously undetected, seismic event.In this work, we decide to focus our effort on the 2009-2014 time period.Starting in 2009 the INGV seismic network underwent a substantial improvement within the CESIS-INGV project, resulting in a substantial improvement in terms of quality (i.e., switching to 3-component instruments) and quantity of its sensors.As an example, between 2009 and 2014 the total number of permanent seismic stations increased by ∼50% (from about 120 to more than 180).The year 2014 has been chosen as the upper limit for our investigation as it was the maximum extent of the CLASS catalog (Latorre et al., 2023) when this work was initiated (see next paragraph for more details).A total of 183 (mostly) permanent stations (Figure 1) are deployed in the area of interest, which includes a portion of the Central Apennines, the Irpinia region, the Calabrian Arc and the NE sector of the Sicily Island.The spatial arrangement of seismic stations is rather inhomogeneous, with a lower station density (inter-station distance of ∼20-30 km) in areas where the seismic hazard is known to be moderate or low (e.g., eastern side of the Apennines and the Apulia region).
The Italian Seismic Bulletin (Bollettino Sismico Italiano, BSI) reports about 22,000 events in the study area between 2009 and the end of 2014.Due to the uneven station coverage, varying azimuthal gap, and uncertainty assigned to the P and S arrivals, only a portion of such events are worth considering as candidate Geochemistry, Geophysics, Geosystems templates, maximizing the chances to find new detections and maintaining a manageable computational expenditure.For the selection of the candidate templates, we employ the newly released CLASS catalog (Latorre et al., 2023) that consists of all events cataloged in the BSI across the Italian Peninsula from 1981 to 2014, re-located with a probabilistic algorithm and a 3D velocity model (Di Stefano & Ciaccio, 2014).First, we selected all events at a maximum depth of 60 km, having a horizontal and vertical error on the hypocenter location less than 2 and 4 km, respectively, with an RMS < 0.8 s and an azimuthal gap <270°.
The selection leads to a total of ∼9,000 events, with a magnitude of completeness (M C ) of 1.9.We construct an ad-hoc 1D velocity model (see Section 2.3) to be used for determining the hypocenter locations of the known seismicity (and later also for the detections) retaining only the events honoring the criteria (see Table 1) for a high-quality earthquake location.This selection led to a total of 3,658 events (M C = 1.9) that are finally employed as templates.We anticipate that almost half of the template events refer to the swarm sequence of the Mt.Pollino (Calabria region), started in the second half of 2012 and peaked with an M W 5.5 mainshock (Totaro et al., 2015).The remaining events concentrate around some well-known seismogenic structures in the Southern Apennines (DISS, 2021).

Running the Template-Matching Algorithm
For this study we employed EQcorrscan (Chamberlain, Hopp, et al., 2017), a Python-based, open-source code designed for the detection of repeating seismicity through template-matching.The algorithm uses normalized crosscorrelation that, in time-domain, at the time sample y for a template t with n samples, is: with t being the mean value of the template signal and d being the mean of nsample long chunk of continuous data.Employing multi-channel data, all the cross-correlation functions of each channel are stacked, taking into account the relative time delay between channels as in the templates.The case in which the stacked cross-correlation function exceeds a certain threshold is indicative of a possible new detection, therefore the P (only on the vertical component) and S picks (on the horizontal ones) are stored.
We choose a threshold of the stacked cross-correlation function that is 7.5 times the Median Absolute Deviation (MAD).Such value resulted from a fine-tuning test between a value of 7 (many events declared per minute, suggesting the detection of many false positives) and 8 (almost no event declared even for a long period).Regarding the cross-correlation value, we imposed a lower limit of 0.65 at each channel.We assign an error to each phase-pick of a detection event based on the value of the original pick uncertainty of its template, linearly weighted on the basis of the cross-correlation value (Valoroso et al., 2023).Being xcorr the cross-correlation, e T the original error on the template pick (commonly 0.1, 0.3, 0.6, or 1 s), we calculate the error on the detection pick e D as follow: Before comparison, both templates and continuous waveforms were subjected to identical processing steps, consisting in mean and linear trend removal.Given our intent to detect hidden low-magnitude seismicity commonly characterized by higher corner frequencies, waveform data were band-pass filtered between 4 and 15 Hz.By appropriately adjusting the minimum allowable cross-correlation and the length of the template, one could aim to find an almost exact replication of a known earthquake.In this regard, we opted for a less conservative approach that allowed a wider range of dissimilarity between a new detection and its template in terms of size, focal mechanism, and hypocentral location.For this reason, we also opted for a short template length, amounting to 0.5 s for the P and S arrivals, plus a pre-pick of 0.1 s (Valoroso et al., 2023).
To optimize the amount of continuous data to scan and decrease the computational burden, we initially discarded all stations at more than 40 km from the template event.This decision was based on the assumption that small-magnitude events occurring at larger epicentral distances could be hardly detected.However, we finally increased such distance to 80 km, since a test with a subset of 1,000 templates showed a 5-fold increase in the total number of new detections with respect to the limit of 40 km that was initially chosen.In fact, by allowing a larger epicentral distance we also increase the chance to detect new events with an adequate amount of arrivals, even in areas with poor station coverage and high magnitude of completeness, such as those in the external part of the Apennines (where seismic hazard is lower).Given the large amount of data to be scanned, we employed the code on a cluster infrastructure using 224 cores.The total computational expenditure amounted to about 6 × 10 5 CPU-hours.

Data Selection, Earthquake Location and Magnitude Estimation
The use of template-matching over the 6 years-long study period led to ∼10 6 triggered detections.We selected only those with (see Table 1) (a) at least four P and four S picks and (b) five or more recording stations, to discard all the events with poor station coverage and an insufficient number of S-wave arrivals for an adequate constraint on the hypocentral depth (Husen & Hardebeck, 2011).The earthquake locations have been determined through the probabilistic, non-linear code NonLinLoc (Lomax et al., 2000) using an ad-hoc, average 1D velocity model.Figure 2 shows the velocity profiles that have been considered.Four 1D models (Matrullo et al., 2013;Pastori et al., 2021;Trionfera et al., 2020;Valoroso et al., 2009) were determined by applying robust inversion schemes (Kim et al., 2006;Kissling et al., 1994) on high-quality travel-time catalogs obtained by dense local networks.Others are based on the results of local earthquake tomography and seismic refraction surveys carried out in the Sannio-Irpinia and Basilicata regions (Improta & Corciulo, 2006;Improta et al., 2014Improta et al., , 2017)).All models show a rather high degree of similarity, and for most models, the V P / V S ratio is constant with depth, within 1.8 and 1.9 across all models.However, to evaluate how the dispersion of the resulting 1D model may affect the uncertainty of hypocentral parameters, we randomly selected and located 1,000 events using 20 randomly sampled models from the mean velocity profile and ±1σ (see Figure 2).The majority of events show deviations for the origin time (<0.1 s) and hypocenter location that are within the uncertainties that typically arise from the errors in the readings of the P and S arrivals.This demonstrates the robustness of our approach when employing a unique, mean and smooth velocity model for the entire study area.
For the estimation of earthquake size, we employ a local magnitude (M L ) scale based on the maximum, half peakto-peak amplitude (A in mm) of the horizontal components.The maximum amplitude is extracted after linear and mean trend removal of the seismic waveforms, deconvolution for the instrument response, simulation of the WA seismometer (sensitivity: 2800), and high-pass filtering (f > 1 Hz).Being R the hypocentral distance in km, M L is calculated following Hutton and Boore (1987): Additional information is in Supporting Information S1 (Text S1).

Results and Discussion
Since the aim of the study is to provide a large, high-resolution, comprehensive catalog, we combine the newly made detections, the templates, and the BSI catalog.For the sake of homogeneity, the hypocenter locations are determined with the same method and 1D velocity model described in the previous section.Moreover, local magnitudes are calculated for detections and re-calculated for templates and BSI using the same method.shows the average 1D model for P, S wave velocities, V P /V S ratio and their dispersion around the mean value (±σ).

The High-Resolution Catalog
The template-matching catalog consists of ∼20,000 new events, a 5-fold increase with respect to the starting 3,658 templates.To evaluate the quality of detections, we compare the quality parameters of the hypocenter solution with those of the templates and the remaining events composing the BSI catalog (Figure 3).Regarding the error on the horizontal component (panel a), the template-matching catalog shows a slightly skewed distribution with a median value of 0.98 km, as opposed to the well-peaked distribution for the templates with 0.59 km as the median value.The two catalogs differ the most when the vertical error (panel b) is considered.The two medians show a difference of ∼0.76 km (1.19 vs. 1.95 km, for the templates and detections, respectively) with a distribution for the detections being again rather skewed.Such discrepancy can be attributed to the substantially different amount of phases (and thus stations) used for obtaining the hypocenter solution for the two different data sets (see panel d).In fact, templates' hypocenters are constrained on average with 60% more arrivals (10 vs. 16) with respect to the detections.As shown in panel (d), the template-matching detections also suffer, consequently, from a slightly higher azimuthal gap (121°vs.96°).Overall, considering the significant amount of picks available for the template events and their better azimuthal coverage, it is remarkable that the events obtained by templatematching show only a slightly higher uncertainty of the hypocenter solution.Such homogeneity of the quality of the two data sets is a solid basis for their combined use in any further analysis for which well-constrained hypocenters are a prerequisite.Considering all the evaluated parameters, the BSI catalog shows an intermediate quality between the templates and the detection.This is easily expected since the templates are, in fact, a highquality subset of all the available events in the study area.Finally, panel (f) in Figure 3 shows the value of the mean cross-correlation of the phase picks for each detection.The majority of events are comprised in the 0.7-0.8range, while ∼1% of the data shows a high degree of similarity with their templates (cross-correlation >0.9).
In Figure S2 of the Supporting Information S1, we plot the mean crosscorrelation between each detection and its template against the two other main predictors for waveform similarity: the difference in magnitude and the inter-event horizontal distance.Panel (b) clearly shows that, within a 5 km distance, both highly similar (cross-correlation >0.9) as well as dissimilar (cross-correlation <0.7) events are found.However, regardless of the distance, almost the totality of high-correlation events are found within 10 7 s (∼4 months), which appears to be a rather sharp threshold beyond which waveform coherence starts to heavily degrade.In Figure S2 of the Supporting Information S1 (panel a), we show the same plot as before except for another predictor of waveform similarity, being the magnitude difference between each detection and its template.As expected, events with the highest correlation show a similar magnitude to their template, with waveform sharply degrading when the magnitude difference increases.The plot shows two wellseparated clusters: (a) the first refers to detections usually within two units of magnitude from their templates, occurring within minutes or a few days from their template (e.g., 10 6 s ≈ 11 days) showing a medium to high value of cross-correlation; (b) the latter (in the order of 10 3 events, be aware of the distortion effect of the logarithmic scale) consists of events significantly smaller (ΔM L < 3), highly dissimilar to their template, and up to 3 years apart from these.Regarding this particular cluster, it is worth noting that the low cross-correlation and large difference in magnitude of its events might raise the concern that such events might be false positives, a possibility that is problematic to discern for such low-magnitude events close to or below the noise level.However, the observation of a rather sharp boundary at an inter-event distance of ∼10 7 s (beyond which the detection of highly similar events becomes unlikely) can be mainly explained as the result of temporal clustering of seismicity.In addition, we speculate that also timevarying characteristics of the propagating medium can play a role in the loss of waveform coherence.For example, significant variations in such properties can occur on a small temporal scale (few months) due to cyclical/seasonal processes (e.g., fluid circulation induced by precipitation) and also on a larger scale due to tectonic or long-term fluid diffusion processes, or a combination of both.For example, time-lapse tomography spanning several years in the Irpinia region (De Landro et al., 2022) has highlighted substantial variation in the V P , V S and V P /V S velocity structure, in the upper as well as in the middle crust.Multi-annual variations of the groundwater level in the shallow karst aquifers have been invoked to explain such changes in the velocity structure, as well as in the crustal deformation and seismicity rate (D'Agostino et al., 2018).
We can quantify the impact of the detection of small-magnitude seismicity in terms of completeness magnitude (Mc), as shown in Figure 4.For this purpose, we employ three different methods (Max Curvature, BestCombo, and MC90) as in Wiemer (2001).In coherence with the small number of selected events, the template catalog has Mc = 1.6-1.9,which is slightly higher than that referred to the entire BSI catalog (Mc = 1.3-1.6).By adding the detection made through template-matching, Mc drops to 0.6-1.5, testifying the improvement in terms of resolution of the resulting earthquake catalog.Such value for Mc can be regarded as rather satisfactory given the characteristics of this study.Only a higher number of templates, a longer time span for the continuous recording and less stringent criteria for the selection of the earthquake location would have assured lower values in magnitude of completeness.

Spatio-Temporal Analysis
The map in Figure 5 (left) shows the earthquake locations determined using the same 1D velocity model (see Section 2.3) for the templates (red), detections (black), and the BSI catalog (blue).Despite the large overlapping of events, it remains evident that template-matching allowed the recovery of a considerable amount of previously missing seismicity.As expected from a method based on measures of waveforms' similarity, such new seismicity falls in the close vicinity of the templates (less than 5 km for 97% of the detections).This is consistent with previous works showing the decay of waveform coherence for increasing distances between templates and detections.As an example, a synthetic sensitivity study in Schaff (2010) shows that detections with an x-corr >0.7 (83% of detections in our case) can be found up to 10-13 km from their template.In a real case scenario, such distance is reduced due to subsurface heterogeneities, differences in event magnitude, fault rupture, focal mechanism, and duration.A map of the density of detections per km 2 is given in Figure 5 (panel b), highlighting the areas with the most prominent recovery of new seismicity, ranging from a few hundred (Matese sequence) to more than 1,000 events per km 2 in the area of the Mt.Pollino sequence.
Templates and detections show the same overall distribution at depth (see panel b in Figure 6) and in both cases, most of the seismicity is confined within 10 km depth.Panel (a) in Figure 6 shows the mean depth of seismicity of the detections in the entire study region, in a 2 × 2 km bin.The map highlights the shallow to very-shallow seismicity (depth ∼5 km) of the large and long-lasting Mt.Pollino sequence, that contributes to the pronounced peak of the depth distribution shown in panel (b) at around 5 km depth.Overall, the depth of the seismicity does not seem to follow a particular spatial pattern, except for the northern part of the study area.Here, in fact, abundant seismicity is found along the axis of the Southern Apennines chain, as well as in its external part toward the foreland domain.Interestingly, earthquake locations are found at increasing depths moving eastwards, thus suggesting a deeper brittle-to-ductile transition and potentially thicker seismogenic layer.Moreover, such variation in hypocentral depth has a magnitude substantially larger than topographic undulation (10-15 km vs. 1-1.5 km) that can be thus considered negligible.Chiarabba and De Gori (2016) suggested that the variable distribution of hypocentral depths along the Apennines can be related to areal changes in heat flow, an explanation that would likely hold for this area as well.The recovery of deep and previously missed seismicity in the external domain of the Southern Apennines can have a direct implication in the evaluation of the seismic hazard of the area, due to the relevant potential for a large magnitude usually associated with the large seismogenic thickness (see Section 3.3).
In Figure 7 (top) we show the cumulative number of events versus time, for detections, templates, and the BSI catalog.The three curves depart from each other as early as in the beginning of 2009, with a sharp increase in the seismicity rate during 2011 in the template-matching catalog, apparently unrelated to any major seismic event.
From the second half of 2012, we observe the highest increase in the number of events due to the highly productive, swarm-like Mt.Pollino sequence that peaked with two M5 events (Cheloni et al., 2017;Passarelli et al., 2015;Totaro et al., 2013Totaro et al., , 2015)).Such increase is visible in the template and BSI catalog, but even more prominent in the detection catalog.From the beginning 2013, while the seismicity rate appears to drop toward the background value, template-matching recovers a large amount of previously missing seismicity at an unprecedented rate with respect to the preceding years.In contrast, such an increase in the recorded seismicity is only mild in the template and BSI catalog, despite the upgrade of the INGV seismic network (∼120 vs. ∼180 active stations from 2009 to 2015).The reason for this likely lies in operational criteria that are currently used for the construction of the official BSI catalog (only events with M > 1.5 are routinely revised by operators), hampering the location of the small magnitude seismicity (Di Maro et al., 2022).This work demonstrates that an a-posteriori analysis of data through template-matching if accompanied by the classical STA/LTA methods commonly used for real-time detection, can lead to a dramatic improvement in the recovery of seismicity.The bottom panel in Figure 7 shows the temporal evolution of magnitudes for templates (red) and detections (black) and the BSI catalog (blue).Most of the detections are within the 0-1 magnitude range.Particularly for the post-2012 period (after the Mt.Pollino sequence) where the seismicity rate is rather uniform, the detections appear to form a background, low-magnitude seismicity unrelated to the M > 3 events occurring in the same period.
In the following sections, we focus on the Sannio-Matese area and the hidden seismicity recovered through template-matching.We show the contribution of our enhanced catalog in the characterization of the local seismicity, in light of the known seismotectonic context.For the sake of completeness, in both cases, we integrated the templates' and detections' catalogs with all the events from the BSI, whose hypocentral location has been determined using the method and the velocity model described in Section 2.3.Then, all events have been relocated with double-difference technique (Waldhauser & Ellsworth, 2000) to improve the imaging of possible fault structures through the well-constrained hypocenter locations.Additional information about the doubledifference location is described in Texts S2 and S3 and related figures of the Supporting Information S1.

The Sannio-Matese Area
The Sannio-Matese area is located at the transition between the central and southern part of the Apennines chain, the east verging fault-and-thrust belt accreted since Miocene due to the subduction of the Adriatic plate (Faccenna et al., 2001(Faccenna et al., , 2014;;Gattacceca & Speranza, 2002).Since the Quaternary, the Apennines belt has experienced a NE-SW directed extension, mainly accommodated by NW-SE striking, newly formed normal faults, or by reactivation of inverted thrust faults inherited from the past collisional phase.As shown in Figure 8 (panel a), this extensional belt is characterized by several M > 6 historical earthquakes, thus posing a substantial seismic hazard (Rovida et al., 2020).According to geological and paleoseismic surveys, the main seismogenic sources are NWtrending normal fault systems bounding the southwestern and northeastern sides of the Matese massif (Boncio et al., 2022): the Aquae Iuliae fault to the west, causative source of the 1349 M6.7 earthquake (P. A. C. Galli & Naso, 2009) and the Bojano fault system to the east, source of the 1805 M 6.6 earthquake (P.Galli et al., 2017;Ferrarini et al., 2017).Moving to SE, an M7 earthquake struck Benevento in 1688, the main municipality in the Sannio region (Figure 8, panel a).In the database of Italian seismogenic sources, the 1688 destructive earthquake is associated with a NW-SE striking and NE dipping extensional structure, but large uncertainties remain regarding the geometry of the causative fault (DISS, 2021), as well as of the sources of the M6+ earthquakes that hit the southern Sannio in 1702, 1732, and 1962 (Figure 8, panel a).

The 2013-2014 Matese Seismic Sequence
On 29 December 2013 a M W 5 event occurred beneath the Matese massif at 17 km depth (see panel a in Figure 8).It was followed by an aftershock sequence that rapidly migrated toward SE (see panel c of Figure S4 in Supporting Information S1).About 22 days later a M W 4.2 earthquake struck in the vicinity of the mainshock (De Gori et al., 2014).In both cases, the fault plane solution suggested a rupture along a NW-striking, steep normal fault that likely dips SW-wards (Ferranti et al., 2015).Afterward, the overall seismicity decayed and completely ceased by February 2014.While the past background seismicity and weak sequences (M < 3) recorded in the Sannio-Matese are confined in the upper crust (e.g., Bisio et al., 2004, among others) the 2013-2014 sequence was unusually deep (15-20 km depth).Given the peculiar finger-like distribution of the seismicity and the type and pattern of amplitude attenuation, this sequence has been related to fluid over-pressure and migration of deep melts rather than tectonic stress release (Di Luccio et al., 2018).The upward migrating aftershocks formed two separate  Westaway (1987).As black dots, the double-difference relocated events detected in this study through template-matching, in addition to their template and other events in the Bollettino Sismico Italiano (BSI).(b) Cumulative number of events in the area of the Matese sequence.The events recovered through template-matching (black) reconstruct the background seismicity of the area, otherwise seismically silent according to both the template (red) and BSI catalogs (blue).clusters, laterally separated by a 1.5 km-wide, dike-like aseismic volume that was interpreted as indicative of a stagnating intrusive body (Di Luccio et al., 2018).In the area of the 2013-2014 sequence (green box in Figure 8, panel a), a total of 312 events are cataloged in the BSI (108 of which have been used as templates in this study).Interestingly, the area appeared seismically silent in the years preceding the 2013-2014 sequence.However, by applying the template-matching technique we obtained a substantial enrichment of the seismic catalog during the 2013-2014 sequence, as well as in previous years.As shown in the cumulative plot of events in Figure 8 (panel b), from 2009 until the M5 mainshock of December 2013, TM allowed us to detect +100 events, as opposed to only four, good quality events, cataloged by INGV in the same time period.Interestingly, less than half of the detections belong to a small seismic burst in the first half of 2009 (see Figure 8, panel b), for which no record exists in the official catalog.The templates that led to such new detections belong to the 2013-2014 sequence, implying an inter-event temporal distance of the template-detection pairs of about 4.5 years.Interestingly, despite the large temporal distance, the detections show a cross-correlation around 0.8, implying a rather high degree of resemblance.This represents a remarkable observation, in light of the generalized dissimilarity between template and detections that is shown to grow in a few months in the rest of the catalog (see Section 3.2 and Figure S2 in Supporting Information S1).
After the minor seismic burst in 2009, template-matching has revealed a rather constant-rate seismicity (Figure 8, panel b) which represents, a previously unnoticed, small-magnitude background activity.The AB profile (Figure 9, panel a) shows the depth distribution of the 2013-2014 Matese sequence (in black) as opposed to the pre-sequence events obtained through TM, occurring as early as 2009.Overall, these events spread across the entire zone where the Matese sequence later occurred, with a noticeable clustering of earlier events (i.e., 2009) in the SE sector and just below the hypocenter of the 2013 mainshock.The pre-sequence events occurring in 2013 are spread in the entire area as well, with some clustering that concentrates in the vicinity of the mainshock.However, the analysis of the space-time distribution of this pre-sequence seismicity does not show any peculiar migration pattern.Strictly considering the sequence period, template-matching recovered several aftershocks that are double in comparison to the official catalog.While this shows a sharp decline in the seismicity rate, implying a ceasing seismic activity in the second half of 2014, template-matching has recovered more than 100 events in the same period.Only at the end of 2014, the number of detections appear to reach a plateau that approaches the level of pre-sequence seismic activity.
Owing to the absence of recorded seismicity before the Matese sequence (see the red and blue cumulative curves in Figure 8, panel b), in their work Di Luccio et al. (2018) excluded the hypothesis of a steady accumulation of magma, upholding the scenario of sudden and episodic dike intrusion and subsequent seismic burst.In light of our results that have unraveled a constant, small-magnitude, background seismicity in the mainshock zone preceding the Matese sequence, we suggest a scenario that is typical of volcanic environment (McNutt & Roman, 2015;Traversa & Grasso, 2010).It generally consists of alternating distribution of volcano-tectonic earthquakes between episodic bursts, with a constant-rate seismicity during the inter-eruptive phases.In this regard, it is instructive to compare the black curve in panel (b) in Figure 8 with the long recording of cumulative volcanotectonic events of the Mt.Etna in Traversa and Grasso (2010), showing a clear pattern of steady background activity and sudden surges in seismicity, corresponding to dike emplacement and subsequent eruptions.An alternative explanation for the episodic seismicity pattern observed in the Matese area relies on the presence of deep-seated fluids.In several areas of the Apennines, relevant emissions and/or circulation of deep-seated CO 2rich fluids have been related to the generation and modulation of low-magnitude seismicity and aseismic deformation (e.g., Valoroso et al., 2017).Other works report on how alternating phases of constant-rate seismicity and mainshock-aftershock sequences can be caused by transient variations of pore fluid pressure at depth (Lucente et al., 2010;Marzorati et al., 2014;Miller et al., 2004).The seismicity pattern during the Matese sequence (Figure S4 in Supporting Information S1) shows a clear along-strike (NW toward SE) migration of the seismicity, suggesting the involvement of fluids.Also, the clear alternation of fast seismic bursts with more quiescent phases (i.e., vertical stripes in Figures S4c and S4d of the Supporting Information S1) might be evidence of aseismic slip deformation triggering part of the seismicity (De Barros et al., 2020).In conclusion, our new data can stimulate a better understanding of such peculiar seismic sequences.As an example, we suggest that the cluster of pre-sequence micro-seismicity in the area of the mainshock could provide evidence of a preparatory phase.Moreover, the retrieval of the pre-sequence, background seismicity might clarify the possible role of fluid circulation and/or other hydrological contributors in modulating the magnitude and spatio-temporal characteristics of seismicity in the Matese area, similar to what is observed for other seismogenic zones in the Apennines (D'Agostino et al., 2018).

Seismicity Distribution in the Sannio Area
The southern Apennines features two different seismotectonic domains: the narrow extensional belt along the range axis and deep shear zones under the external Apennines toward the foreland in the east (Ciaccio et al., 2021;Di Luccio, Fukuyama, & Pino, 2005;Di Luccio, Piscini, et al., 2005;Fracassi & Valensise, 2007;Pino et al., 2008).In the former domain, moderate to large earthquakes show extensional mechanisms striking NW-SE (1980M6.9 Irpinia earthquake, 1998M5.6 Castelluccio, 2013-2014 M5.0 Matese sequence, 2011-2014 Pollino sequence) and the seismogenic layer is about 10-15 km thick (Chiarabba et al., 2005;Latorre et al., 2023).With the only exception of the Matese seismic sequence described in the previous section, this observation is clearly evidenced in Figure 6 which displays a clear concentration of events at a depth around 5 km, followed by a sharp decline below 10 km.In contrast, the external Apennines domain has been characterized by moderate strike-slip earthquakes (1990M5.7 Potenza, 2002M5.9 Molise, Pondrelli et al., 2006) along E-W structures at mid-to lowercrustal depths within the Apulian plate.However, such a regional seismotectonic zonation seems inappropriate for the region encompassing the southern Sannio and northern Irpinia, which was struck by M6+ earthquakes in 1702, 1732, 1962 (see Figure 8, panel a).For the 1962 sequence, published focal solutions for the two M ∼ 6 shocks are debated and vary from NW-striking transtensional mechanism Westaway (1987) to E-W striking pure strike-slip mechanism Vannoli et al. (2016), although their hypocenters have been located at upper crustal depths (<10 km).In September 2012, a M W 4.2 earthquake occurred in this region at a depth of ∼10 km (http://iside.rm.ingv.it/event/1335371) and published focal solutions varies again from a NW-trending transtensional mechanism (Scognamiglio et al., 2006) to a pure strike-slip mechanism trending E-W (Adinolfi et al., 2015).Extensional to strike-slip rupture mechanisms also characterize the local background seismicity recorded in 2014-2016(De Matteo et al., 2018).Based on these observations, the Sannio-Irpinia boundary has been regarded as a complex transitional zone where the seismotectonic setting changes both in the SW-NE direction as well as in-depth (Adinolfi et al., 2015;De Matteo et al., 2018), with an extensional regime confined in the first 15 km, superimposed to a transcurrent one at mid-lower crustal depth that characterizes the external part of the belt.Our new template-matching catalog provides valuable information on the complex seismotectonic setting of this region.The CD profile in Figure 9 (panel b) that runs SSW-NNE from the axial to the external portion of the Apennines reveals a clear NE-wards deepening of seismicity.In the southernmost part of the profile (about 10 km NE of the city of Benevento), seismicity is confined in the upper crust (<15 km depth) and appears to concentrate in two distinct depth levels, 5-9 km and 11-15 km.In the outer belt, two remarkable vertical clusters are found at substantially greater depth (20-25 km), consistent with the overall deepening of seismicity shown in Figure 6, when moving toward the external Apennines.The template-matching catalog displays, for the first time, that the eastward deepening of the seismicity is extremely rapid: the cut-off of the background seismicity increases from 15 to 25 km depth in about 10-12 km distance (Figure 9, panel b), considering that the profile D is not exactly perpendicular to the belt axis.If we take into account our re-location of the 2012 M W 4.2 earthquake, this is coherent with the western seismicity confined in the upper crust.For both the mid-crustal clusters, seismicity persisted during the whole 2009-2014 interval, illuminating two vertical 2018)).In conclusion, previous works and our new results, taken collectively, indicate that an abrupt deepening and change in the tectonic style, from normal to strike-slip, occurs in the NE direction at the border between the southern Sannio and the northern Irpinia region.Coherently, the comparison of our template-matching catalog with the stress field inversion of De Matteo et al. ( 2018) favors an extensional rupture mechanism for the debated 1962 M6 and 2012 M w 4.2 earthquakes.The emergence of two vertical, km-scale structures with strike-slip dextral kinematics, at mid-crustal depth in the Apulian plate has a significant implication for seismic hazard.Such structures may represent silent seismogenic faults that are analogs to the dextral E-W shear zones that ruptured SE-wards and NE-wards of the Sannio-Irpinia region, such as the 1990 M5.7 Potenza and the 2002 M5.9 Molise deep sequences.

Conclusion
We have produced the first high-resolution seismic catalog at regional scale for the Southern Apennines (Italy), employing a template-matching (TM) technique, detecting about 20,000 new seismic events.Located with the same smooth, average 1D velocity model, the templates and the new detections show comparable quality of the hypocenter solution.The events retrieved by template matching have magnitudes confined in the 0-1 range, reaching values as low as 0.5, thus leading to an improvement of the magnitude of completeness that, depending on the method employed, shows a decrease of about one unit with respect to the template catalog.
To show the potential of the template-matching catalog for the understanding of the spatio-temporal seismicity distribution and its possible impact on the assessment of the seismic hazard, we focus on the Sannio-Matese area.
(i) Regarding the 2013-2014 Matese sequence, TM successfully recovered a substantial amount of previously undetected seismic activity, revealing a sudden surge of seismicity more than 3 years before the Matese sequence.Additionally, we retrieved the small-magnitude, constant-rate, background seismicity occurring near the M5 mainshock and subsequent 2013-2014 sequence.Such background activity was not reported in the starting catalogs.In fact, this area was considered seismically silent for the time preceding the 2013-2014 sequence.(ii) Located at the border of the Southern Apennines belt axis toward the outer belt, the Sannio area shows a rapid deepening of the seismic activity unraveled by template matching.Moving from SSW to NNE (from the internal to the external portion of the Apennines belt), hypocenters concentrate between 20 and 25 km, rather than 5-10 km.Such observation, combined with evidence on the focal mechanism of past events, suggests that the Sannio area can be regarded as a complex area characterized by a sharp thickening of the seismogenic layer and a change of the deformation style, from pure extensional to transtensive/strike-slip.

Figure 1 .
Figure 1.Map of the study area (Southern Apennines, Italy), with the major seismogenic faults (associated with events with M > 5.5) as reported in the DISS database (https://doi.org/10.13127/diss3.3.0).Red triangles represent the +180 permanent stations from the INGV seismic network employed for this study.

Figure 2 .
Figure 2. 1D velocity models in terms of V P (a) and V P /V S ratio (b) that have been used in this study to build an average, 1D model (black solid line).Panel (c and d) shows the average 1D model for P, S wave velocities, V P /V S ratio and their dispersion around the mean value (±σ).

Figure 3 .
Figure 3. Quality parameters of the hypocenter solutions for the template-matching catalog (black), the templates catalog (red), and the Bollettino Seismico Italiano (BSI, blue): (a) horizontal error, (b) vertical error, (c) root-mean-square error (RMS) on arrivals (d) number of seismic phases, (e) azimuthal gap.(f) Mean crosscorrelation of picks for all detections.

Figure 4 .
Figure 4. Cumulative frequency-magnitude distribution for the templates catalog (red), BSI (blue).The distribution of the full catalog (black) is obtained by adding the detection to the BSI catalog.

Figure 5 .
Figure 5. (a) Map of the events of the Bollettino Sismico Italiano (BSI) in blue, templates (red), and template-matching detections (black).The orange stars represent the M L > 3.5 events that occurred between 2009 and 2014.The gray boxes highlight the location of the two main seismic sequences in the same time period.(b) Density plot of the detection 1 × 1 km bin.Black lines indicates major seismogenic faults (M > 5.5) as in the DISS database (https://doi.org/10.13127/diss3.3.0).

Figure 7 .
Figure 7. (a) Cumulative curves of templates (red), events cataloged in the BSI (blue), and detections (black).Orange stars indicate the events with M > 3.5 occurring in the study area between 2009 and 2014.(b) Distribution of magnitude versus time of events for templates, BSI, and detections.As in the panel above, events with M > 3.5 are superimposed.

Figure 8 .
Figure 8.(a) Map of the Sannio-Matese area.Beachballs indicate the focal mechanism of the three main events that occurred in the area between 2009 and 2014.The 2013-2014 Matese sequence is highlighted with a green box while in yellow the Sannio area is indicated.Superficial traces of capable faults are in red (ITHACA Working Group, 2019).Blue squares indicate the proposed epicenters for large historical earthquakes (e.g., 1805 and 1688).The focal mechanism of the 1962 M6.1 earthquake is fromWestaway (1987).As black dots, the double-difference relocated events detected in this study through template-matching, in addition to their template and other events in the Bollettino Sismico Italiano (BSI).(b) Cumulative number of events in the area of the Matese sequence.The events recovered through template-matching (black) reconstruct the background seismicity of the area, otherwise seismically silent according to both the template (red) and BSI catalogs (blue).

Figure 9 .
Figure 9. (a) Depth profile along the AB trace (see Figure 8, panel a) cutting the 2013-2014 Matese sequence.The pre-sequence events are color-coded according to the year of occurrence.Aftershocks are indicated in black.The red stars represent the mainshocks of the sequence.(b) Depth profile along the CD trace in the Sannio area, note the substantial deepening of the overall seismicity.
structures covering more than 5 km in depth.A comparison with the 2014-2016 catalog of M < 3 events in De Matteo et al. (2018) unravels that these two vertical clusters were also active in the following years.Moreover, the focal mechanisms of the best-located events of De Matteo et al. (2018) in the two deep vertical clusters show a clear E-W trending, strike-slip dextral mechanisms, whereas the seismicity located NE of Benevento and confined in the upper crust have NW-SE trending extensional kinematics (see Figure 4 in De Matteo et al. (

Table 1
Thresholds Used for the Selection of Both Templates and New TM-Detections DIAFERIA ET AL.