Machine‐Learning‐Based High‐Resolution Earthquake Catalog Reveals How Complex Fault Structures Were Activated during the 2016–2017 Central Italy Sequence

On 24 August 2016, a moment magnitude (⁠$Mw$⁠) 6.0 earthquake struck near the town of Amatrice (Fig. 1), resulting in the death of 299 people (Italian Civil Protection). It was followed two months later by an $Mw$ 5.9 earthquake near the town of Visso on 26 October, before culminating in an $Mw$ 6.5 earthquake near the town of Norcia four days later (Fig. 1). The Norcia earthquake is the largest earthquake in Italy since the 1980 $Mw$ 6.9 Irpinia earthquake (Boschi et al., 2000). Subsequently, four $Mw$ 5.0–5.5 events occurred in the Campotosto area on 18 January 2017. This sequence of moderate‐to‐large earthquakes along an 80 km long normal‐fault system devastated the surrounding towns and villages and left 20,000 people homeless (Galli et al., 2016). Seismic sequences characterized by the occurrence of multiple $Mw>5$ earthquakes seem to be a characteristic of the central Apennines. For instance, the 1703 sequence around the Norcia region involved multiple $M>6$ earthquakes within 20 days (Boschi et al., 2000). More recently, the 1997 Colfiorito sequence involved six $Mw>5$ earthquakes occurring over two months (Chiaraluce et al., 2003), whereas the 2009 L’Aquila sequence involved five $Mw>5$ earthquakes occurring over one week (Valoroso et al., 2013). The Colfiorito and L’Aquila sequences occurred immediately to the northwest and southeast of the 2016–2017 Amatrice–Visso–Norcia (AVN) seismic sequence (Fig. 1).

Understanding the factors controlling the evolution of earthquake sequences is critical for hazard assessment and operational earthquake forecasting (Mancini et al., 2019). Because the 2016–2017 AVN sequence is well recorded by a dense network of seismometers (Fig. 1), it has been the focus of numerous studies. The routine catalog produced by the Italian National Institute of Geophysics and Volcanology (INGV; Italian Seismological Instrumental and Parametric Database [ISIDe] Working Group, 2007) contains $∼82,000$ events related to this sequence over a 1 yr period (Fig. 2). Chiaraluce et al. (2017), analyzing a subset of $∼26,000$ events, inferred that the largest earthquakes nucleated at the base of a southwest‐dipping normal‐fault system segmented by crosscutting compressional structures and bounded at $∼8km$ depth by a shallowly east‐dipping shear zone (SZ). Vuan et al. (2017) analyzed the spatiotemporal evolution of $∼3000$ earthquakes in the eight months before the Amatrice earthquake and found that most events were located within the SZ, with the seismic activity in the SZ being asynchronous compared with that in the shallower upper crust. They inferred this to reflect slip on the SZ loading the overlying normal‐fault system, which results in subsequent unlocking. Walters et al. (2018) interpreted the seismicity distribution for the entire AVN sequence, together with geodetically estimated slip distributions, to suggest that the rupture extent and termination of the three largest events were controlled by intersecting subsidiary faults that acted as structural barriers. However, detailed geometry of the various synthetic, antithetic, and cross faults, and their role during the largest events’ coseismic ruptures are still under debate (Cheloni et al., 2017, 2019; Chiaraluce et al., 2017; Scognamiglio et al., 2018; Walters et al., 2018; Improta et al., 2019; Michele et al., 2020).

High‐resolution catalogs that are complete to small magnitudes are crucial for mapping 3D fault structures and probing earthquake nucleation and triggering processes. In this article, we present an enhanced earthquake catalog for the AVN sequence. We constructed the catalog using a deep‐neural‐network‐based phase picker (Zhu and Beroza, 2019) combined with high‐precision double‐difference relative relocation (Waldhauser and Ellsworth, 2000). The resulting 1 yr catalog contains 900,058 events, more than a 10‐fold increase compared to the INGV routine catalog and is complete for earthquakes of $ML>0.3$ (Fig. 2). Our catalog reveals rich details regarding the complex fault structures that were activated during the AVN sequence.

We use 1 yr of continuous data from available stations within 80 km of the epicentral area, starting from 15 August 2016. These include permanent stations that are part of the Italian National Seismic Network and stations temporarily deployed shortly after the Amatrice event by the INGV and the British Geological Survey (Moretti et al., 2016; see supplemental material available to this article for details). From 10 September 2016, we have a stable network of 139 stations until the end of our catalog period, with average station spacing in the epicentral region of $∼8km$ (Fig. 1).

We use PhaseNet (Zhu and Beroza, 2019), a deep‐neural‐network‐based picker, for detecting earthquakes, and picking their P‐ and S‐arrival times (Figs. S1 and S2). We then associate the phase picks with individual events using the Rapid Earthquake Association and Location package (Zhang et al., 2019). Using a subset of 5000 events with at least 120 associated picks, we utilize Velest (Kissling et al., 1994) to estimate optimal 1D P‐ and S‐velocity models with station corrections (Fig. S3), starting with the 1D velocity models in Chiaraluce et al. (2017). We then locate all events with the HypoInverse software (Klein, 2002) using the optimal 1D velocity models before relocating events with at least four P picks and seven total picks using the HypoDD double‐difference method (Waldhauser and Ellsworth, 2000). Our final catalog contains 900,058 events (Fig. 2), and local magnitude (⁠$ML$⁠) was estimated for these events. We further evaluate the PhaseNet pick uncertainties by comparing pick‐derived and cross‐correlation‐derived relative arrival times for a subset of events with similar waveforms (Fig. S4). We find that the estimated pick uncertainties are comparable with previous studies (e.g., Waldhauser et al., 2020). Finally, we estimate relative location errors for each event using bootstrapping (Waldhauser and Ellsworth, 2000). The horizontal minor and major axes, and the vertical projections of the 95% confidence ellipsoids (Fig. S5) have median values of 36, 56, and 87 m, respectively. Details about the catalog production and uncertainty estimation procedures can be found in the supplemental material.

INGV’s routine catalog has a magnitude of completeness (⁠$Mc$⁠) of 2.3 (Mancini et al., 2019), when estimated using the goodness‐of‐fit method and requiring that 95% of the data can be modeled by a power‐law fit (Wiemer and Wyss, 2000). Using the same method, we estimate $Mc=0.3$ for our catalog; however, the magnitude–time plot clearly shows time‐varying $Mc$ (Fig. 2a). The detection of smaller magnitude events deteriorates immediately after the largest earthquakes when the seismic data are saturated by the mainshock and aftershock codas before recovering to the background detection level with time (Fig. 2a, inset), similar to what was observed for the INGV routine catalog (Mancini et al., 2019). In addition, the detection rate of the smallest events shows daily fluctuation (Fig. 2a, inset). Looking at events that occurred at least six months after the Norcia earthquake, we find that there are about 50% more events around midnight compared with around noon (Fig. S6a). The minimum $ML$ of events that occurred around midnight are also, on average, about 0.3 smaller compared with events that occurred around noon (Fig. S6b). These observations, which are not apparent in INGV’s routine catalog (Fig. S6c,d), suggest that the smallest events detected by PhaseNet are probably close to the background noise level, and that the detection threshold fluctuates with diurnal changes in cultural noise level. The time‐varying $Mc$ will have to be carefully accounted for during subsequent analysis of temporal variations in $b$‐value.

Using the maximum‐likelihood method that accounts for binned magnitudes (Utsu, 1966), we obtain a $b$‐value of 0.86, with the $b$‐value increasing when using larger $Mc$ values. This is because the frequency–magnitude curve shows a deviation from linearity at large magnitudes (Fig. 2b), which likely reflects a change in $ML–Mw$ scaling (Malagnini and Munafò, 2018). We apply the scaling relation derived by Grünthal et al. (2009) for Europe to convert $ML$ to $Mw$ for our catalog (see the supplemental material). This derived $Mw$ is more suitable for subsequent analysis of $b$‐value variations, because the frequency–magnitude curve obeys linearity for the entire magnitude range above $Mc$ (Fig. 2b).

Our enhanced catalog provides an opportunity to further illuminate the geometry of complex fault structures that were activated during the AVN sequence, as well as their role in the development of the sequence (Cheloni et al., 2017, 2019; Chiaraluce et al., 2017; Scognamiglio et al., 2018; Walters et al., 2018; Improta et al., 2019; Michele et al., 2020). The large number of events means that simply displaying all earthquakes as point clouds in spatial plots, as is often done for smaller catalogs, is no longer the best way to display the seismicity for discerning complex fault structures. Therefore, we first apply the hierarchical density‐based spatial clustering of applications with noise (HDBSCAN) algorithm (McInnes et al., 2017) to filter out diffuse seismicity and thus accentuate structures delineated by high‐density clusters of seismicity (Fig. S7). The algorithm is an extension of the density‐based spatial clustering of applications with noise algorithm (Ester et al., 1996), which has previously been utilized to identify seismicity clusters to delineate faults (e.g., Schoenball and Ellsworth, 2017) but with the added advantage of not having to arbitrarily set a single search radius for the clustering process (see the supplemental material for details). We then make map‐view and cross‐section (including events within 1 km of the profiles) plots of the high‐density seismicity clusters color‐coded by time to illuminate the complex fault structures activated during the AVN sequence (Figs. 3 and 4).

The cross section near the $Mw$ 6.0 Amatrice earthquake hypocenter shows that its aftershocks delineate a southwest‐dipping plane that aligns with mapped surface rupture (Fig. 3b). However, the aftershocks concentrate down‐dip of the mainshock hypocenter, and the map‐view and along‐strike cross‐section plots show a $∼15km$ long section with lower aftershock density (Fig. 3a,f). This is consistent with finite‐fault inversion results in which the main slip patches are up‐dip of the mainshock hypocenter and extend $∼15km$ along strike (e.g., Tinti et al., 2016; Chiaraluce et al., 2017; Cirella et al., 2018; Walters et al., 2018), with previous studies having observed that the main slip region is relatively devoid of aftershocks (e.g., Improta et al., 2019; Michele et al., 2020). In comparison, a large number of aftershocks delineate a subhorizontal feature at $∼8–10km$ depth (Fig. 3b,f), which has previously been reported (e.g. Chiaraluce et al., 2017; Improta et al., 2019; Michele et al., 2020) and inferred to represent a buried thrust related to the rise of the Apennines (Carannante et al., 2013), an east‐dipping low‐angle normal fault (Lavecchia et al., 2017) or the brittle portion overprinting a zone characterized by strain localization and eventually ductile deformation (Chiaraluce et al., 2017). We will subsequently refer to this feature as a SZ.

To the north, the aftershocks delineate both a southwest‐dipping plane and an additional northeast‐dipping plane that align with mapped surface fault traces (Fig. 3c). The northeast‐dipping plane has previously been inferred to be an antithetic fault that hosted a large aftershock (e.g., Chiaraluce et al., 2017; Improta et al., 2019; Michele et al., 2020). In this region, the SZ appears thinner and dipping eastward gently (Fig. 3c). In map view, it can be seen that later aftershocks extended northward of the mapped surface normal‐fault traces and delineate a series of lineation that are westward of but parallel to the Olevano–Antrodoco–Sibilini (OAS) thrust surface trace (Fig. 3a). Taking a cross section perpendicular to the ramp portion of the OAS, the aftershocks delineate a feature that is dipping northwest, with dip angle decreasing with increasing depth (Fig. 3e). Improta et al. (2019) previously observed a similar aftershock alignment and inferred the structure to be an inherited thrust that was reactivated as a normal fault, but our enhanced high‐precision catalog provides a better constrain on the geometry of this structure. We find that the shallowest portion at $∼2km$ depth is more steeply dipping than the presumed low‐angle dip of inherited compressional structures, and its surface projection is $∼4km$ northwest of the OAS surface trace beneath the Mt. Vettore normal‐fault surface trace (Fig. 3e). Nevertheless, it is possible that this structure represents the OAS thrust ramp at depth, with the shallowest portion not illuminated by aftershocks having a shallower dip angle that connects to the surface trace. This geometry would be consistent with the Pizzi et al. (2017) geological interpretation. Alternatively, this structure could be another northwest‐dipping fault (Walters et al., 2018), such as a splay of the OAS thrust ramp. We note that the seismicity delineating this feature forms a relatively narrow band (Fig. 3e) in comparison with seismicity delineating the main normal‐fault planes (Fig. 3b–d) and, thus, is unlikely to represent the projection of the apparent dip of nearby southwest‐dipping normal faults. The western extension of this structure is also apparent in the along‐strike cross‐section plot (Fig. 3f, outlined by dashed ellipse), bounding the northernmost extent of the Amatrice earthquake’s aftershock zone. This structure appears to merge into and segment the SZ (Fig. 3e,f), with the SZ north of this structure located at a shallower depth than in the south (Fig. 3f).

The cross section near the $Mw$ 5.9 Visso earthquake hypocenter shows that its aftershocks delineate a southwest‐dipping plane that aligns with mapped surface rupture (Fig. 3d). The aftershocks concentrate up‐dip of the mainshock hypocenter and in the SZ, with the main slip region inferred from geodetic inversion (Walters et al., 2018) also relatively devoid of aftershocks (Fig. 3a,d). East and westward of the main fault plane, there are steeply dipping structures that are delineated by the Amatrice earthquake’s later aftershocks (Fig. 3d), with the structure to the east possibly the northern extension of the potentially reactivated thrust ramp (Fig. 3a).

The cross section near the $Mw$ 6.5 Norcia earthquake hypocenter shows that its aftershocks delineate a southwest‐dipping plane that aligns with mapped surface rupture and terminates into the SZ (Fig. 4b). The aftershocks also delineate a northeast‐dipping structure in the hanging wall at $<3km$ depth, which has previously been inferred to be an antithetic fault (Chiaraluce et al., 2017; Improta et al., 2019) that potentially ruptured coseismically (Walters et al., 2018). There is also a hint of another antithetic fault in the footwall. However, both antithetic structures were already seismically active before the Norcia earthquake (Fig. 4b).

In January 2017, a sequence of $Mw>5$ earthquakes activated seismicity in the Campotosto region. The Campotosto fault, which was activated during the 2009 L’Aquila sequence (Valoroso et al., 2013), is reactivated (Fig. 4e). However, the seismicity generally delineates faults that appear complementary to faults activated during the L’Aquila sequence (Fig. 4a,d,f). This is also the case when comparing seismicity from the AVN sequence with those from the 1997 Colfiorito sequence (Fig. 1). This suggests that earthquake triggering primarily occurs on critically stressed faults, and there is still a significant likelihood that future earthquake sequences can occur on unmapped blind faults that lack recent seismicity even in regions where there has been recent seismic activity.

The role of inherited thrusts in controlling the evolution of the largest earthquakes’ coseismic ruptures, in reason of their oblique intersection with younger normal faults and distribution of the more competent lithologies, is still under debate (Cheloni et al., 2017, 2019; Scognamiglio et al., 2018; Walters et al., 2018; Michele et al., 2020; Barchi et al., 2021). For the Amatrice earthquake, its two main slip patches (Tinti et al., 2016; Cheloni et al., 2017; Chiaraluce et al., 2017; Cirella et al., 2018; Walters et al., 2018) are located on either side of the OAS surface trace in map view. Its immediate aftershock zone and inferred rupture plane, illuminated by a zone that is relatively devoid of aftershocks, also crosses the OAS surface trace in map view (Fig. 3a). We also observe a dearth of seismicity in the top 6 km of the crust between the OAS and the Gran Sasso (GS) thrusts (Fig. 5), consistent with Chiaraluce et al. (2017), in spite of the $∼30$ times more events included. These observations appear to suggest that the OAS and GS thrust ramps dip at high angles and segment the normal faults into two distinct structures, but appear inconsistent with the hypothesis that the OAS thrust ramp acted as a barrier to the northward rupture propagation of the Amatrice earthquake (e.g., Pizzi et al., 2017; Improta et al., 2019). However, cross‐section plots along strike (Fig. 3f) and perpendicular to the OAS thrust ramp (Fig. 3e) show a northwest‐dipping structure that marks the northernmost extent of the Amatrice earthquake’s aftershock zone. If this structure is part of the OAS thrust ramp, it would imply that the OAS thrust ramp has a shallow dip near the surface but steepens at depth, hence at depth it marks the northernmost extent of the Amatrice aftershock zone, but in map view the Amatrice aftershock and rupture zones extend across its surface trace. Alternatively, this structure could be a different northwest‐dipping fault (Walters et al., 2018), such as a splay of the OAS thrust ramp.

Whether the coseismic rupture of the Norcia earthquake involved an ancillary northwest‐dipping fault plane is still under debate (e.g., Cheloni et al., 2017, 2019; Scognamiglio et al., 2018). Interestingly, the northwest‐dipping structure we observed to have been activated during the Amatrice earthquake aftershock sequence (Fig. 3e,f) was not as strongly activated during the Norcia earthquake’s immediate aftershock sequence (Fig. 4c,f). This would be consistent with the structure accommodating coseismic slip of the Norcia earthquake, if the inverse relation between slip and aftershock density observed for the Amatrice and Visso earthquakes (Fig. 3a,b,d,f) also applies here. Nevertheless, resolving whether inherited thrust structures accommodated coseismic slip and acted as barrier during the AVN sequence, and if so which structures and how many were involved, will require further detailed studies that precisely coregister the relative locations of seismicity, coseismic slip, and geological structures in 3D.

We also observe that the extent of the SZ activated after the AVN earthquakes (Figs. 3–5) generally coincides with the mainshocks’ rupture extent. This suggests that seismic activity on the SZ responded passively to triggering by the shallow normal‐faulting earthquakes and did not play an active role in controlling the progressive activation of different shallow normal‐fault segments during the AVN sequence. This differs from what was proposed by Vuan et al. (2017) for seismic activity before the start of the sequence, suggesting that coseismic deformation during the AVN sequence might have changed the coupling relationship between the SZ and the shallow normal faults.

Overall, we find that applying the HDBSCAN algorithm to first filter out diffuse seismicity better accentuates fault structures compared with simply plotting all events as point clouds (Figs. S8 and S9). Our enhanced catalog also more clearly delineates complex fault structures, especially the geometry at depth of a northwest‐dipping fault that intersects the normal faults obliquely, when compared with the most complete, double‐difference relocated version of the INGV catalog published to date (Figs. S10 and S11), which consists of $∼34,000$$ML>1.5$ events (Michele et al., 2020). Although the main fault planes delineated by the seismicity in our catalog appear principally planar, dipping southwest at $∼45°$⁠, and not listric (Figs. 3b–d and 4b), we find that the fault planes delineated by seismicity are relatively thick compared with the earthquakes’ estimated relative location uncertainties, with some showing structures suggesting varying fault orientations within the deformation zone (e.g., Fig. 4d). There also appears to be a large number of diffuse seismicity, for instance, $∼30%$ of seismicity was unclustered by the HDBSCAN algorithm (Fig. S7d). Although the proportion is dependent on the chosen parameter (see the supplemental material), the parameters were manually tuned to pass the eye test in terms of isolating high‐density clusters (Fig. S7). In addition, this region is populated by numerous kilometer‐long fault segments (e.g., north of 43° N), many of which are subvertical (e.g., Figs. 3d and 4d). Ongoing work incorporating waveform cross‐correlation aims to investigate how the complexities of the deformation zones relate to properties such as fault age and roughness, as well as the inherited structures.

The investment in supplementing the permanent network with additional seismic stations shortly after the $Mw$ 6.0 Amatrice earthquake has provided a dramatically improved view of the evolution of a complex and protracted earthquake sequence. Our enhanced catalog revealed additional insights on how complex fault structures were activated during the AVN sequence, including their detailed geometry and distribution, their role during the rupture of the largest events, and how they relate to recent nearby earthquake sequences. Future studies further leveraging this enhanced catalog will hopefully lead to better understanding of the controlling factors behind these fault‐zone complexities and how they impact earthquake triggering and fault activation.

The earthquake catalog presented is available at the Zenodo dataset repository (doi: 10.5281/zenodo.4662870). Data from the permanent network (Network code: IV) are available at the European Integrated Data Archive (EIDA; https://www.orfeus-eu.org/data/eida/, last accessed March 2021). Data from temporary stations deployed by British Geological Survey (Network Code: YR; doi: 10.7914/SN/YR_2016) are available on the IRIS Data Management Center.

Seismic waveform data used in this study were downloaded from the Incorporated Research Institutions for Seismology (IRIS) and the National Institute of Geophysics and Volcanology (INGV) Data Centers. The Global Centroid Moment Tensor Project database was searched using www.globalcmt.org/CMTsearch.html (last accessed December 2020). The INGV seismic bulletin is available at http://cnt.rm.ingv.it/ (last accessed December 2020). Details about the catalog production, uncertainty estimation, and seismicity clustering procedures are provided in the supplemental material.

The authors declare no competing interests.

The authors thank Ian Main and Danielle Spallarossa for fruitful discussions. This work is supported by the United Kingdom's National Environment Research Council (NERC)‐United States National Science Foundation, Directorate for Geosciences (NSFGEO) funded project The Central Apennines Earthquake Sequence under a New Microscrope (NE/R0000794/1). G. Beroza was supported by Department of Energy (Basic Energy Sciences; Award DE‐SC0020445). Y. J. Tan was also supported by the Direct Grant for Research (Grant 4053483) from the Chinese University of Hong Kong. The deployment of the temporary U.K. British Geological Survey (BGS) stations was enabled by NERC direct funds under an instrument loan number 1067 and 1077 provided by the Geophysical Equipment Facility in collaboration with SEIS‐UK, led by M. Segou.