Centroid Moment Tensor Catalog With 3D Lithospheric Wave Speed Model: The 2016–2017 Central Apennines Sequence

Moment tensor inversions of broadband velocity data are usually managed by adopting Green's functions for 1D layered seismic wave speed models. This assumption can impact on source parameter estimates in regions with complex 3D heterogeneous structures and discontinuities in rock properties. In this work, we present a new centroid moment tensor (CMT) catalog for the Amatrice‐Visso‐Norcia (AVN) seismic sequence based on a recently generated 3D wave speed model for the Italian lithosphere. Forward synthetic seismograms and Fréchet derivatives for CMT‐3D inversions of 159 earthquakes with Mw ≥ 3.0 are simulated using a spectral‐element method (SEM) code. By comparing the retrieved solutions with those from time domain moment tensor (TDMT) catalog, obtained with a 1D wave speed model calibrated for Central Apennines (Italy), we observe a remarkable degree of consistency in terms of source geometry, kinematics, and magnitude. Significant differences are found in centroid depths, which are more accurately estimated using the 3D model. Finally, we present a newly designed parameter, τ, to better quantify and compare a‐posteriori the reliability of the obtained MT solutions. τ measures the goodness of fit between observed and synthetic seismograms accounting for differences in amplitude, arrival time, percentage of fitted seconds, and the usual L2‐norm estimate. The CMT‐3D solutions represent the first Italian CMT catalog based on a full‐waveform 3D wave speed model. They provide reliable source parameters with potential implications for the structures activated during the sequence. The developed approach can be readily applied to more complex Italian regions where 1D models are underperforming and not representative of the area.

. Map view of the study area. Dots represent earthquakes from Michele et al. (2020) while beachballs are focal mechanism solutions for the 136 events (among 159) with at least a "fair" quality (see Section 3) obtained in this study; both are color-coded by depth. Enlarged beachballs with red stars are events with M w ≥ 6.0; beachballs with black stars are events with 5 ≤ M w < 6.0. Black lines represent cross-section profiles. The black dashed box refers to the events described in Section 3.2. Red line represents the Olevano Antodoco Sibillini (OAS) thrust front trace (modified after Centamore and Rossi (2009)) while blue lines represent mapped normal faults (Pucci et al., 2017, and references therein). It is well-known that a high-quality focal mechanism catalog is of crucial importance to obtain good constraints on regional stress field, to assess earthquake hazards, and to understand tectonic processes. The usual procedure to determine MTs for small to moderate earthquakes considers simple local or regional 1D seismic wave speed models (Herrmann et al., 2011;Scognamiglio et al., 2009;Yang et al., 2012). However, some regions are characterized by strong 3D heterogeneities, which can range from local to regional scales (Takemura et al., 2021;Wang & Zhan, 2019), and the adoption of 1D models may cause errors and unstable MT solutions (Hingee et al., 2011;Scognamiglio, Magnoni, et al., 2016;Wang & Zhan, 2020). In some cases, 3D Earth structures have to be adopted to take into account the nonuniform distribution of the stations used to perform the inversion. The effects of the lateral Earth variations may also be of crucial importance when considering small-to-moderate-sized earthquakes with high signal-to-noise ratio (SNR) waveforms only at short periods (Wang & Zhan, 2019;Zhu & Zhou, 2016).
Seismic tomography and full-waveform inversion have made incredible progresses. Thanks to recent improvements in computational resources as well as forward solvers (e.g., Peter et al., 2011), the generation of very accurate 3D models has become feasible. The inclusion of a good 3D wave speed model, when available at the appropriate resolution scale, could be now considered as a reasonable and standard procedure to estimate seismic source solutions.
One of the most important results when using 3D wave speed models for computing MT solutions is to obtain better constraints of the kinematic parameters (Hingee et al., 2011), earthquake size, and depths through a better waveform fit (Covellone & Savage, 2012;Nayak & Dreger, 2018). Takemura et al. (2021) show that a 1D wave speed model systematically overestimates the M w when compared with their local 3D model, while Hjörleifsdòttir and Ekström (2010) and Wang and Zhan (2019) point out that depths are better constrained when considering a 3D instead of 1D wave speed model. Wang and Zhan (2020) conclude that MT solutions from 3D wave speed models help in identifying the geometry and kinematic of the subsurface activated structures. Well-constrained off-shore solutions are also of great importance. Takemura et al. (2018Takemura et al. ( , 2020, computing centroid moment tensor (CMT) inversions of earthquakes along the Nankai Trough (Japan), demonstrate that using a 3D wave speed model improves the source parameter estimate for off-shore earthquakes, especially in terms of dip angles and centroid depths. Finally, an increase of double-couple component is also documented in regions with complex heterogeneous structures when 3D wave speed models are adopted (Covellone & Savage, 2012;Hejrani et al., 2017;Jechumtàlovà & Bulant, 2013;Wang & Zhan, 2019).
In this study, we present reviewed source geometries of the recent AVN sequence as retrieved by MT analyses performed for small-sized to moderate-sized earthquakes and based on a new 3D wave speed model of the Italian lithosphere (Magnoni et al., 2022). Using this sequence in Central Italy, we can take advantage of the large number of M > 3.0 events recorded by numerous well-equipped seismic stations, and of the TDMT catalog needed as input for our calculation. The availability of the new 3D model represents a great step forward to get more reliable MT solutions in terms of estimated source parameters, fault geometry, and faulting styles. Moreover, obtaining 3D source solutions that compares favorably with the TDMT catalog based on a 1D wave speed model calibrated for this area, can support the robustness of our new 3D catalog. In the future, this will also allow us to derive MT solutions for Italian regions where a 1D model is underperforming.
In our work, we introduce a new quality parameter, in order to quantify the goodness of the obtained solutions. This parameter evaluates the quality of the data fit and is calculated by considering the differences between recorded and synthetic data in cross-correlation, amplitude, timing, and length of fitted seconds (Section 2.3). The performance of the solution is no longer defined only by a Variance Reduction (VR, as usually done when assessing the quality of fit) but takes into account other characteristics of the signal.

3D Centroid Moment Tensor Inversion
In this section, we present the data set and procedures used to compute the 3D Centroid Moment Tensor (CMT-3D) solutions, and how we evaluate the solution quality.

Data Set
We examine 159 earthquakes with M w ≥ 3.0 belonging to the studied area (

Inversion
We follow the point-source MT inversion procedure proposed by Liu et al. (2004; https://github.com/UTComp-Seismo/GRD_CMT3D) by using the python code "pycmt3d" (https://doi.org/10.5281/zenodo.56124). Starting from an existing, reliable focal mechanism solution, this CMT inversion technique minimizes a waveform misfit function between data and synthetics by constructing the misfit function variation from numerically calculated Fréchet derivatives with respect to the considered source parameters (Liu et al., 2004).
In order to simulate full synthetic waveforms, we use the spectral-element method (SEM) code SPECFEM3D_ Cartesian (Peter et al., 2011), which allows for accurate simulations in complex heterogeneous media. Considering the TDMT solutions as starting solutions, the synthetic seismograms for the considered events are computed together with synthetic seismograms for perturbed source parameters, which are needed to construct the Fréchet derivatives for the six MT components and the three source location parameters.
Since our goal is to produce an MT catalog based on a reliable Earth structure, the wave speed model used in the simulation code is the 3D model Im25 (Magnoni et al., 2022), which has been recently obtained for the Italian lithosphere by an adjoint 3D full-waveform travel-time tomography. Im25 is a regional-scale model for V p and V s with an unprecedented spatial resolution which, considering the retrieved wave speed values, corresponds to a minimum period of ∼10 s (Figures S1 and S2 in Supporting Information S1). For this model, the quality factor Q is obtained as a linear function of V s , and the values of density ρ are calculated as a quadratic function of V p based on an empirical relationship (Magnoni et al., 2022).
Before pycmt3d, we use the FLEXWIN windowing code (Maggi et al., 2009) in order to select the time windows, on recorded and synthetic seismograms, suitable for the MT inversion. Only earthquakes with at least 10 time windows are inverted (Magnoni et al., 2022). The time window selection is done via user-tunable parameters by imposing threshold values for cross-correlation, amplitude ratio, and time-shift (see definitions in Maggi et al. (2009)) between synthetics and observables. We choose the following requirements of goodness of fit to be satisfied by these quantities: cross-correlation ≥0.8, |amplitude ratio| ≤ 0.8, |time-shift| ≤ 5 s or |time-shift| ≤ 7 s for M L < 3.95 or M L ≥ 3.95, respectively. In addition to the above parameters, we modified the original FLEX-WIN code by adding a new condition: vr ≥ 0, where vr is the window-VR defined as Here, d win and x win are the data and synthetic time series within the window.
In FLEXWIN, we also check the SNR within each time window, by imposing a minimum threshold of 4.0. With this strong constraint, 30 earthquakes among 159 are excluded from the analysis. For these events, we lower the SNR threshold for the time windows to 3.0 and carefully verify the goodness of the selected windows.
The code pycmt3d performs the MT inversion considering only the selected windows for each pair of real and synthetic seismograms. A zero-trace constraint is always imposed, thus implying the isotropic MT component to be zero.
During the MT inversion, the code perturbs the initial solution to explore the space parameters and uses the synthetics for the perturbed source parameters, simulated with SPECFEM3D_Cartesian, to construct the corresponding Fréchet derivatives. The value of the perturbation for each MT component is chosen as the maximum order of magnitude of the MT components for the given earthquake. For most of the events, latitude and longitude perturbation is 0.18°, and depth perturbation is 8 km. For three shallow events, we reduce the depth perturbation to 3 km to avoid nonphysical solutions, which locate earthquakes above the local topography. For the mainshocks' hypocenters, which are already widely studied, we reduce the latitude, longitude, and depth perturbation to 0.045° and 5 km, respectively. The synthetics for the new CMT solutions are constructed by pycmt3d as a combination of the Fréchet derivatives for the parameter perturbations finally resulting from the inversion.
Using 150 CPU cores on the INGV MERCALLI cluster (Intel Xeon Gold 5218 processors), the total average computation time to produce a CMT-3D solution is ∼600 CPU-hrs, the time for the FLEXWIN window selection is up to 2 CPU-hrs, increasing as a function of the station number, and a negligible time is required for the inversion part as well as the processing.

Evaluation of the Solutions
In order to provide an a-posteriori quantitative estimate of the MT solution quality, we define a parameter τ inspired by the metric proposed by Covellone and Savage (2012). τ quantifies the capability of the source solution of modeling the real data and therefore of giving a reasonable estimate of the seismic source parameters: good solutions correspond to small τ values. For each analyzed event, we compute τ as = ∑∑ .
(2) C τ is a sigmoid-like function such that its value is ∼1 for a number of windows greater or equal to 10, while the value increases for a smaller number of windows here, x is the number of windows, a is the center of the sigmoid function, while b is the width. With this choice, C τ = 2 if the number of windows is 8. We verified that solutions with <10 windows are often unreliable and set the above parameters accordingly. Exceptions are represented by very good values of the other parameters composing τ (see below) that can balance the small number of windows (e.g., 8 or 9) yielding anyway a value of τ < 1. The parameter in Equation 2 is defined as where the bar symbol stands for the average of the considered parameters. Here, 1 − vr is based on the definition in Equation 1, Tshift and CC are the FLEXWIN time-shift and cross-correlation parameters, respectively. To define the parameter dA, we applied some algebra starting from the definition of the amplitude ratio dlnA in FLEXWIN (Maggi et al., 2009) From Equation 6, we note that 1 − dA is directly proportional to the difference of the seismic wave energy. Finally, the parameter fits = 1 − N fitted /N tot , where N fitted is the length of the inverted window and N tot is the total analyzed seconds for each event.
The weights w m of each parameter are calculated as the standard deviation except for fits to which we assign w fits = 0.1.
The averages of the parameters in Equation 4 are computed considering their values in each window of the given event, and are weighted by the duration of the windows. In the cases of dA and Tshift, we also normalize the averages by their maximum possible values. For Tshift, we choose T max = 7, while for dA, we have that the condition |dlnA| ≤ 0.8 (Section 2.2) implies −3.95 < < 0.8, thus we choose dA max = 4. The final form of τ (expressed as a weighted average) as well as the selected thresholds result from the FLEXWIN parameters chosen for the window selection, and from the fact that we consider all the parameters equally important. An example of fit, including the window parameters used to evaluate τ, can be found in Figure S3 in Supporting Information S1.
Based on the well-defined physical meaning of each parameter in τ (Equation 2), we consider an "excellent" fit for values of τ from 0 to 0.25, a "good" fit from 0.25 to 0.50, a "fair" fit from 0.50 to 0.75, and for values above 0.75 an "unresolved" fit.

Results
We calculated a total of 159 CMT-3D solutions for small to moderate earthquakes occurred during the AVN seismic sequence. Estimated moment magnitudes M w range from 3.01 to 6.43 and depths from −0.6 to 15.1 km (see Table S1 and the ".csv" Table in Supporting Information S1 for source parameter values and associated uncertainties). Most of the retrieved solutions show NW-SE normal fault mechanism confirming the wellknown northeast-trending tectonic extension of this portion of the Apennines (Boncio et al., 2004;Carminati & Doglioni, 2012;Chiarabba et al., 2005;D'Agostino et al., 2009, for a review). Nevertheless, strike-slip, transtensional and low-angle normal fault mechanisms exist, confirming the heterogeneous pattern of seismicity ( Figure 1).
Following the newly developed metric τ for MT solution quality estimation, we found that 83.1% of the retrieved solutions show a good or excellent fit between data and synthetics (τ ≤ 0.50, green beachballs in Figure 2), while 14.5% have τ ≥ 1.00 and a poor fit caused by the small number of the selected time windows (<8, displayed as unresolved solutions, red beachballs). The remaining 2.4% solutions show a satisfactory fit (yellow beachballs). The MT unresolved solutions belong to the smallest magnitude earthquakes of the catalog or to events occurred immediately after the Amatrice and Norcia mainshocks. Indeed, the time proximity of a mainshock causes the overlap and interference of phases from the two events.

Comparison With the TDMT Catalog
To further evaluate the reliability of our CMT-3D solutions, we compare source mechanisms and moment magnitudes to those from the TDMT catalog (http://terremoti.ingv.it/tdmt), one of the reference catalogs for the Italian seismicity.
Mechanisms are compared through the Kagan rotation angle, which is defined as the minimum rotation about any axis that is needed to transfer from one focal mechanism to the other (Kagan, 1991). This angle ranges between 0° and 120° and allows us to quantitatively describe the dissimilarity between two focal mechanism solutions. We found an excellent consistency between the two catalogs ( Figure 3a): <7% of the focal mechanisms have a rotation angle larger than 30°, which is generally considered as the threshold for a good agreement (Bernardi et al., 2004). Similarly to the large τ values, most of the larger Kagan angle values are mainly related to solutions with a small number of inverted windows and/or small magnitude events. Kinematically diverging solutions come out with a high τ value, confirming the suitability of the quality parameter introduced for MT quality estimation.
As an additional test, we computed the Omega angle distance (Steinberg et al., 2021) that takes into account the 6 MT components (see Figure S4 in Supporting Information S1). We observe that almost all the events have an Omega angle distance <0.1° confirming the kinematic consistency between TDMT and CMT-3D solutions.  more among different catalogs is often documented and is mainly attributed to the different velocity structures, procedures and frequency range. This result shows that, at least in the same frequency range, the two models and the two procedures match very well, supporting M w estimates obtained independently from each other. This is a noteworthy outcome, since M w is the reference estimate of earthquake size, both for scientific and nonscientific community. The 3D Centroid Moment Tensor (CMT-3D) M w values are also plotted together with their uncertainties (see Text S1 in Supporting Information S1 and the ".csv" Table in Supporting Information S1), an important feature, often disregarded, which contributes to determine the reliability of the estimates. As better shown in Figure 3b inset, only 7 events among the plotted 136 have uncertainties larger than 0.2 (green dots).
For scientists, a realistic and well-constrained M w value is needed for computing ground-shaking scenarios or when dealing with ground motion assessment, because it directly affects the amplitude of the simulated ground motion. Instead, common people perceive the magnitude as an absolute and perpetual value so discrepancies on this parameter could lead to misunderstanding and debates (La Longa et al., 2014;Scognamiglio, Tinti, & Quintiliani, 2016).
We also compared our centroid depths to those derived from TDMT, whose solutions are inferred exploring a range of depths discretized every 1 km around the initial depth ( Figure 3c). We observe a large discrepancy between the two catalogs, probably due to a limited depth resolution of the TDMT procedure, together with the use of a simplified 1D model, although well calibrated for the area. Our results suggest that incorporating the 3D wave speed model in source inversions can improve the focal depth assessment. For sake of readability, we avoided plotting uncertainties on CMT-3D depth values, but they are all reported in the ".csv" Table in Supporting Information.
With an additional analysis, we studied if the retrieved M w values satisfy the scaling relationship between M L and M w recently proposed for Central and Northern Apennines by Malagnini and Munafò (2018). For the 49 coinciding events, by comparing our M w estimates with their M L , computed by using a regionally calibrated attenuation relation (Munafò et al., 2016), we fit within one standard deviation the bilinear regression with the crossover at M L = 4.3 proposed by these authors (Figure 3d, and Text S1 in Supporting Information S1). This result further corroborates the reliability of the obtained M w estimates which consequently behave as the M w values of Malagnini and Munafò (2018).
As a further comparison, in Figure 4, we present the fit of both 3D and 1D models to the real data of the event no. 2 (24 August 2016, 02.33.29, M w 5.3). The figure shows that the CMT-3D synthetics' fit to the real data is often very good and that the 3D wave speed model allows us to add stations excluded by the TDMT procedure. This result stresses the need of including the 3D wave speed model into the MT estimation to get more constrained solutions.

Earthquake Hypocentral Location: Constraining the Tectonic Structures
In order to further constrain the subsurface fault geometries and faulting styles of the tectonic structures activated during the AVN sequence, we compare our new CMT-3D catalog with the double difference relative locations of Michele et al. (2020), the best available at the time of writing. Locations and focal mechanisms from MT solutions are overlapped in map-view, color-coded by depth (Figure 1), but also shown in a longitudinal cross-section, oriented N151°, and seven vertical cross-sections, oriented orthogonally to the longitudinal one (Figures 5-7). Vertical cross-sections include events within 2 km (i.e., ±2 km) from the profiles.
For common events, we find a general good agreement between the epicentral locations. Histograms of longitude and latitude distances are shown in Figure S5 in Supporting Information S1. We observe a slightly displaced location toward northeast for CMT-3D solutions (top panels), but we get a median of 0.0073° for the latitudes and 0.0094° for the longitudes from the distribution of distances in absolute value (bottom panels), indicating an evident good agreement. We highlight that the comparison of the location has been done only for the events with M w ≤ 5.0, in order to make reasonable the comparison between our CMT-3D centroid and hypoDD locations. The map in Figure 1 shows that most focal mechanisms match in color with the background seismicity and coherently characterize the extensional kinematic of the main fault plains and their antithetic structures (Cheloni et al., 2019;Scognamiglio, Tinti, & Quintiliani, 2016Tinti et al., 2016;Walters et al., 2018).  Table S1 in Supporting Information S1. If red waveforms are not present, this means that the TDMT are not capable to fit the real data of the corresponding stations.
Moving from NW to SE along the AVN sequence strike direction, we discover the relationship between aftershocks and focal mechanisms. Aftershocks in cross-section A-Aʹ delineate a synthetic-antithetic fault system (Figure 5a). Focal mechanism locations and the dip angles agree with the event distribution and highlight that this system hosted some of the larger magnitude earthquakes occurred in this area (up to M w = 4.1). The shallower northwestern events nucleated on the antithetic fault, while the deeper ones, between 3 and 6 km, are located close to the synthetic fault. The deepest normal mechanism is the M w 4.09 event, the first M w ≥ 4.0 event occurred in this area. It occurred few hours after the Visso mainshock and is located close to the low-angle fault shear zone (SZ) which has been previously identified by Chiaraluce et al. (2017), Michele et al. (2020), and Lavecchia et al. (2017). This seems to support the hypothesis that "seismic activity on the SZ responded passively to triggering by the shallow normal-faulting earthquakes" (Tan et al., 2021). The lowest dip angle of the considered mechanism is 33° and is related to the SW-dipping plane, that is in the opposite direction with respect to the dip delineated for the SZ by the studies on aftershocks Improta et al., 2019;Michele et al., 2020). This feature has been already observed by Michele et al. (2020) and agrees with the idea that the events nucleating close to the SZ may dislocate higher angle planes as a brittle response to the subhorizontal deformation process. The A-Aʹ section also contains three strike-slip focal mechanisms. The role of the strike-slip kinematic cannot be easily inferred from the aftershocks' distribution, since no vertical event alignment is found even if we change the cross-sections' azimuth according to their strike values. These earthquakes are probably associated to minor structures that accommodate the overall tectonic deformation.
Moving toward southeast, we explore the correlation between MT solutions and the faults where Visso and Norcia earthquakes nucleated. Seismicity in the section crossing the Visso mainshock (Figure 5b) slightly describes a southwest-dipping plane. The portion of the fault around the Visso hypocenter, located at 4.6-km depth, appears almost as a shadow area for the aftershock distribution. This feature could be explained with the commonly observed anticorrelation between aftershocks and main slip release fault area. Indeed, the finite-fault model proposed by Chiaraluce et al. (2017) shows that the slip is entirely distributed on an elongated patch of 8 × 4 km at depths between 3.5 and 6.5 km that well agrees with the lack of aftershocks. This portion of the fault system hosts two CMT-3D solutions, the mainshock, whose centroid location and focal mechanism are perfectly consistent with the fault plane, and a shallower mechanism which could suggest a dip-slip southwest-dipping kinematic for the eastern events aligned in parallel to the mainshock's fault.
In section C-Cʹ, which hosts the M w 6.5 Norcia earthquake, the aftershock distribution images a quite complex geometry of activated faults whose kinematic is explained by the retrieved CMT-3D solutions (Figure 5c). These solutions clearly describe the mainshock's fault plane dipping 47°SW. Two other CMT-3D solutions (immediately on the right of the mainshock) lie along the shallow part of this highlighted plane and confirm the geometry and the kinematic of this structure. It is interesting to note that the two solutions belonging to the main fault plane are related to the seismicity active before the Norcia earthquake. Both in the hangingwall and in the footwall of the mainshock's fault, a clustered seismicity suggests the existence of two NE-dipping antithetic faults and the CMT-3D solutions located close to these faults are consistent with this hypothesis. Once again, the eastern strike-slip retrieved solutions cannot be  Table S1 in Supporting Information S1). associated with vertical aftershock alignments. These events belong to the Amatrice earthquake's aftershocks and demonstrate the temporal and geometrical complexity of the fault system activated during the AVN sequence. Finally, the deeper focal mechanism solution, occurred 2 days before the Norcia mainshock, could be related to the activity of the nearby basal SZ. The role of this low-angle east-dipping SZ is still a matter of debate Improta et al., 2019;Lavecchia et al., 2017;Michele et al., 2020;Tan et al., 2021;Vuan et al., 2017). However, the CMT-3D solution is not geometrically coherent with the SZ, its east-dipping fault plane dips 42°, too much for a low-angle structure.
In order to characterize the kinematic of the faults activated between the Norcia and Amatrice domain, we crosscut the seismicity through the largest Amatrice aftershock (Figure 5d), an area marked by a high-density of CMT-3D solutions. In this section, the aftershocks delineate two opposite dipping structures that overlay the subhorizontal seismicity which is slightly dipping toward the NE, and very well illuminated here. Unfortunately, the MT solutions located on or near this tectonic structures are not consistent to each other and do not help discriminating the kinematic behavior of the structures. On the contrary, the NE-dipping alignment of kinematically coherent MT solutions, outlined by the dashed blue line in Figure 5d, perfectly fits with the antithetic fault suggested by Scognamiglio, Tinti, & Quintiliani (2016) for the Amatrice aftershock and partially depicted by the aftershocks' location. The preferred focal plane dips ∼50°. All these earthquakes are aftershocks of the Amatrice mainshock.
Moving down toward southeast through the Laga Mts. fault system, section E-Eʹ hosts the M w 5.96 Amatrice earthquake. Very few aftershocks image the mainshock's fault; however, the retrieved mainshock's MT solution perfectly fits the hinted 50°SW-dipping plane, as well as the other two displayed CMT-3D solutions (Figure 6a). The lack of aftershocks in the first 7-8-km depth is consistent with the location of the main slip patches of the Amatrice finite-fault source model .
South of the Amatrice ruptured area, two cross-sections (F-Fʹ and G-Gʹ) describe the geometry and the kinematic of the faults involved in the southern termination of the AVN sequence. Section F-Fʹ (Figure 6b) is characterized by a main SW-dipping fault compatible with the dip-slip kinematic shown by the nearby CMT-3D solutions. The shallowest focal mechanisms have a preferred nodal plain dipping between 50° and 60°, while the deepest mechanisms show dip values ranging between 30° and 40°. The observed lowering dip angle seems to confirm the suggestion by Michele et al. (2020) that, in this area, shallow seismicity flattens with depth. The CMT-3D solutions, located in the small cluster of earthquakes on the fault hangingwall, helps to interpret this cluster as a minor antithetic fault. In Figure 6c Finally, section H-Hʹ gives an overall view of the along-strike geometry and kinematic of the tectonic structures involved in the AVN sequence, and of the extension of the three mainshocks' rupture areas (Figure 7). Along this section, the retrieved CMT-3D solutions help to describe the rupture characteristics of the numerous activated

Discussion and Conclusions
We generated a new CMT catalog for the 2016-2017 AVN sequence by inverting three-component waveforms within the recent 3D lithospheric wave speed model for Italy Im25. The CMT-3D catalog contains 159 solutions with moment magnitudes ranging from 3.01 to 6.43 and depths from −0.6 to 15.1 km (see Table S1 in Supporting Information S1 and the extended Table in Supporting Information).
The quality of the solutions is provided by the newly developed parameter τ, which quantifies the capability of the CMT-3D solutions of modeling the real data, therefore giving a reliable estimate of the seismic source parameters. We found that ∼83% of the retrieved solutions shows a good or excellent fit between data and synthetics, while ∼14% shows higher value of τ mainly caused by the small number of selected time windows.
The presented CMT-3D catalog contains the uncertainties on the components of the MT, which translate into uncertainties on the derived parameters allowing us to perform some statistical analysis on the inverted source parameters including the seismic moment, and the time and spatial location parameters.
We compared the M w determined in our study to that derived from the TDMT catalog, which is based on an ad-hoc 1D wave speed model for Central Italy. Differences between the M w estimates are subtle, with almost all the values contained in a very small interval (±0.1).
The M w values have been compared also with Malagnini and Munafò (2018) through the M w − M L regression that confirms an excellent agreement among the M w estimates.
A significant enhancement of the CMT-3D solutions is the inferred events' depth. Discrepancies with respect to the earthquake focal depths obtained by TDMT result from a raw exploration in depth (with a 1-km increment) of that procedure.
As a result, inversions for the location parameters, differently from the TDMT technique, give us the possibility to investigate the features and the role of the structures activated during the sequence, and to compare our location parameters with the most accurate available location catalogs.
Different MT solutions for the same event are often attributed to different velocity structures, procedures and the adopted frequency range. While the latter represents the scale length at which we are illuminating the earthquake process, the first two factors can determine/include epistemic and aleatoric uncertainties. The kinematic consistency among the 1D and 3D catalogs, documented through the small values of Kagan angle, Omega distance, and magnitude discrepancy, allows us to say that for this region both the adopted models and procedures are capable of providing realistic and well-constrained solutions. The CIA model is built ad-hoc for the Central Italy region, and is the best model we currently have for this region. The excellent agreement with CMT-3D solutions found in this study suggests that the resolution of the 3D model will allow us to retrieve reliable solutions also in areas where the CIA model is less performing and not representative (i.e., Sicily, Adriatic region).
The proposed CMT-3D catalog can be used in full-waveform inversion to improve the resolution of 3D Earth models, in studying the temporal and spatial variations of stress conditions, and the depth distribution of seismicity. It can also contribute to explain the complexity of the seismogenic processes active in the Central Apennines and help comprehend the main features of the seismic sequences. The improvement of our understanding on the activated fault systems could have potential implications to mitigate seismic hazard and risk in the area.

Data Availability Statement
Data from regional seismometers are available via FDSN services from EIDA. The INGV Italian Seismic Network earthquake catalog, along with metadata and other ancillary data, such as MTs and focal mechanisms has been used and is available at http://terremoti.ingv.it. The complete solutions of CMT-3Ds are included in the Supporting Information.