**Contents**


Reprinted from: *Remote Sens.* **2022**, *14*, 2217, doi:10.3390/rs14092217 ................. **189**


## **About the Editors**

## **Chisheng Wang**

Chisheng Wang is currently an associate professor with the School of Architecture and Urban Planning, Shenzhen University, Shenzhen, Guangdong, China. He received a B.S. degree from Beijing Normal University, Beijing, China, in 2007; an M.S. degree from the Institute of Applied Remote Sensing, Chinese Academy of Science, Beijing, in 2010; and a Ph.D. degree from the Department of Land Surveying and Geo-Informatics, Hong Kong Polytechnic University, Kowloon, Hong Kong. His research interests include disaster and infrastructure monitoring, InSAR, point cloud processing, and photogrammetry.

## **Daqing Ge**

Daqing Ge is currently a Senior Engineering (Professorship) with the China Aero Geophysical Survey and Remote Sensing Center for Land and Resources, Beijing, China. He was born in Shaanxi, China. He received a B.S. degree in engineering of surveying and mapping from the China University of Mining and Technology (CUMT), Xuzhou, China, in 2002; an M.S. degree in photogrammetry and remote sensing from the China University of Mining and Technology (Beijing) (CUMTB), Beijing, China, in 2005; and a Ph.D. degree from the China University of Geosciences (Beijing) (CUGB), Beijing, China, in 2013. He has served as the Principal Investigator of several Chinese national projects. His research interests mainly include radar interferometry technology, geohazards monitoring applications, and the development of SAR satellite application systems.

## **Guohong Zhang**

Guohong Zhang is a research professor in the Institute of Geology, China Earthquake Administration. He received his PhD degree in Geophysics from Institute of Geology, China Earthquake Administration. His current research interests include the history of earthquake ruptures using joint inversion of InSAR and seismic data, co- and post-seismic crustal deformation using InSAR technology and its modeling, the active deformation of the continents, and the rheology of the lithosphere. He has more than 80 journal and conference publications.

## **Wu Zhu**

Wu Zhu is currently a professor with School of Geology Engineering and Geomatics, Chang'an University, Xi'an, China. He received an M.Sc. degree in geodesy and geodynamics from Chang'an University, Xi'an, China, in 2010 and a Ph.D. degree from the Department of Land Surveying and Geo-Informatics (LSGI), Hong Kong Polytechnic University, Kowloon, Hong Kong, in 2015. His research interests include the estimation and correction of ionospheric effects on interferometric synthetic aperture radar (InSAR) and analyses of InSAR-measured surface deformation associated with geophysical hazards.

## **Siting Xiong**

Siting Xiong is currently an associate researcher with the Guangdong Laboratory of Artificial Intelligence and Digital Economy (SZ), Shenzhen, China. She received a B.Sc. degree from the School of Geodesy and Geomatics, Wuhan University, Wuhan, China, in 2011; an M.Sc. degree from the Department of Earth and Space Sciences, Peking University, Beijing, China, in 2014; and a Ph.D. degree from the Department of Space and Climate Physics, University College London, London, U.K., in 2019. She was a post-doctoral researcher with the College of Civil and Transportation Engineering, Shenzhen University, Shenzhen, China, from 2019 to 2022. Her research interests include image processing, radar subsurface mapping, interferometric SAR, and planetary radar data processing. She has more than 30 journal and conference publications.

*Article*

#### **Rupture Models of the 2016 Central Italy Earthquake Sequence from Joint Inversion of Strong-Motion and InSAR Datasets: Implications for Fault Behavior**

**Chuanhua Zhu 1,2,3, Chisheng Wang 1,2,\*, Xinjian Shan 3, Guohong Zhang 3, Qingquan Li 1,2, Jiasong Zhu 1,4, Bochen Zhang 1,4 and Peng Liu 5**


**Abstract:** We derived the joint slip models of the three major events in the 2016 Central Italy earthquake sequence by inverting strong-motion and InSAR datasets. *b*-values and the historic earthquake scarp offset were also investigated after processing the earthquake catalog and near-field digital elevation model data. The three major events gradually released seismic moments of 1.6 × 10<sup>18</sup> Nm (Mw 6.1), 1.5 × 10<sup>18</sup> Nm (Mw 6.1), and 1.1 × 10<sup>19</sup> Nm (Mw 6.7), respectively. All the ruptures exhibit both updip and along-strike directivity, but differ in the along-strike propagation direction. The high *b*-value found beneath three mainshock hypocenters suggests possible fluid intrusions, explaining the cascading earthquake behavior. The cumulative surface scarp from past earthquakes shows rupturing features that are consistent with the 2016 earthquake sequence, suggesting a characteristic fault behavior. Under the assumption of the Gutenberg–Richter law, the slip budget closure test gives a maximum magnitude of Mw 6.7 and implies the seismic hazard from the largest event has been released in this sequence.

**Keywords:** source inversion; rupture model; strong motion; InSAR; seismic hazard

## **1. Introduction**

5

Earthquakes have sustainedly caused a large number of casualties and damages worldwide. Seismological research relies heavily on the inversion for the earthquake source, including the source location, fault geometry, slip distribution, rupture directions, etc. Using that, we can explore fault behaviors and seismic hazards [1–4]. As such, inversion for earthquake sources has long been one of the most popular topics after a destructive earthquake. Benefitting from the algorithm progressive and data coverage, joint inversion of seismic and geodetic data have been successfully applied for detailed source rupture processes of earthquakes worldwide, e.g., the 1999 Mw 7.1 Duzce, the 2011 Mw 9.0 Tohoku, and the 2012 Mw 7.6 Nicoya, Costa Rica, earthquake, etc. Additionally, it has been proved to be able to provide a stabler and higher resolution source model than a single-data inversion [5–8].

On 24 August 2016 (UTC 01:36, local time 03:36), a destructive earthquake (Mw 6.2) occurred in central Italy (the Amatrice earthquake). The USGS reported that the earthquake originated at a depth of 4.4 km, with a normal faulting mechanism. The epicenter was

**Citation:** Zhu, C.; Wang, C.; Shan, X.; Zhang, G.; Li, Q.; Zhu, J.; Zhang, B.; Liu, P. Rupture Models of the 2016 Central Italy Earthquake Sequence from Joint Inversion of Strong-Motion and InSAR Datasets: Implications for Fault Behavior. *Remote Sens.* **2022**, *14*, 1819. https:// doi.org/10.3390/rs14081819

Academic Editor: Salvatore Stramondo

Received: 20 March 2022 Accepted: 8 April 2022 Published: 10 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

1

located at 42.70◦N, 13.23◦E, between the towns of Norcia and Amatrice. According to the official figures of the Protezione Civile, this event caused the death of 297 people, with 234 of the casualties occurring in the town of Amatrice. Two months later, two major earthquakes were triggered northwest of the Amatrice earthquake on 26 October (19:18 UTC) and 30 October (19:18 UTC), with individual moment magnitudes of 6.1 (Visso earthquake) and 6.6 (Norcia earthquake). There were no reports of serious injuries in the October earthquakes.

The ground motion associated with the Amatrice earthquake was recorded by many geodetic and seismic instruments, including space-borne synthetic aperture radar (SAR) sensors and ground seismometers. These datasets can be used to constrain a fault slip model, which can help us to understand the earthquake mechanics and seismic hazard. To date, source models have been inferred using one or several datasets [9–16]. For example, Tinti et al. (2016) presented the first kinematic model of this event by inverting the waveforms from 26 3-component strong-motion accelerometers [14]. Lavecchia et al. (2016) derived a static slip model constrained by differential interferometric SAR (DInSAR) measurements from several SAR satellites [12]. However, the features of these models differ from each other because different datasets were used in these studies. As InSAR datasets and strong-motion datasets have complementary strengths in earthquake source inversion [17], it is possible to invert both datasets simultaneously for a more comprehensive fault slip model.

The future hazard of a fault after an earthquake is of grea<sup>t</sup> importance for human life and property safety. The 2016 central Italy earthquake sequence is located in a segmented fault system mixed with modern and ancient structures. Conventional seismic hazard models assume earthquake ruptures are controlled by fault segmentation, where the rupture is believed to be unlikely propagate from one segmen<sup>t</sup> to another. In the 2016 sequence, the heterogeneous crust has broken the lateral continuity of the seismogenic fault, prohibiting it to grow to a large devastating earthquake. However, increasing observations, e.g., the 2010 EI Mayor-Cucapah earthquake [18] and the 2016 Kaikoura earthquake [19], show that multi fault segments are possible to simultaneously rupture in a single earthquake. In addition, the segmented faults may become more continuous and mature after many earthquakes repeatedly ruptured, and consequently, inclined to host a larger earthquake. Therefore, the future seismic hazards in this region need to be re-examined based on more observations. Besides the comprehensive rupture model, we obtained as mentioned above, the complete Italian earthquake catalogue and available near-field high-resolution topography enable us to retrieve more fault information including *b*-value, historical events trace, and seismic activity. A detailed interpretation of the causative fault by combining the information can therefore benefit the earthquake hazard evaluation.

In this paper, we first present the rupture processes of the three major earthquakes in the central Italy earthquake sequence by inverting joint datasets, including both InSAR and strong-motion data. We then explore the *b*-values and historic earthquake scarp offsets through processing the earthquake catalog and near-field digital elevation model (DEM) data. Based on these results, we interpret the fault behavior of the earthquake sequence and discuss the future seismic hazard in this area.

## **2. Tectonics**

The 2016 central Italy earthquake sequence took place in the northern and central Apennines of Italy, which were formed due to the collision between the Adria microplate and the Eurasian plate. The collision induced compression in the front arc and extension in the back arc. The back arc of the Apennines chain is thus characterized by a set of NW–SE-striking normal fault systems with a current overall east–west extension rate of 1.5–3 mm/yr [20,21]. The earthquake sequence occurred across several fault systems (Figure 1), including the Mt. Vettore–Mt. Bove Fault system (VBFS) and the Mt. Gorzano Fault system (GFS).

**Figure 1.** Overview of seismicity, faults and SAR data coverage for the 2016 central Italy earthquake sequence. (**a**) Seismotectonic setting. The three major earthquakes in the sequence are denoted by the red stars and beach ball symbols, and the two recent large historical events are colored in black. The dashed black rectangles represent the surface projection of the fault planes adopted in this study. The bold black lines are the two seismogenic normal fault systems, namely, the Mt. Vettore–Mt. Bove Fault (VBFS) system and the Mt. Gorzano Fault (GFS) system. The pink line shows the simplified trace of the preexisting compressional front named the Olevano–Antrodoco–Sibillini (OAS) thrust. Aftershocks are marked by the small dots with different colors. The blue, purple, and yellow points represent the events taking place between 24 August 2016–26 October 2016, 26 October 2016–30 October 2016, and 30 October 2016–8 October 2017, respectively. (**b**) Map view of the SAR data coverage. The black rectangles represent the coverage of the Sentinel-1A and ALOS 2 SAR data for generating the interferograms. The red rectangle denotes the area shown in (**a**). (**c**) Cross section through A-A' shown in Figure 1a. The dashed black line represents the fault plane of the 30 October 2016 Mw 6.7 Norcia earthquake.

The VBFS is a SW-dipping fault system, composed of different splays and segments. The fault system scarps are exposed along the SW foothills of Mt. Vettore, Mt. Porche, and Mt. Bove, with a length of ~27 km. The GFS is a SW-dipping fault, with a ~26 km long fault scarp developing along the foothills of Mt. Gorzano. These two faults are segmented by existing tectonic structures inherited from pre-Quaternary compressional tectonics [22]. A ~3 km thick layer in which small events and some large extensional aftershocks occur is found below the seismogenic fault, limiting the seismicity to the first 8 km of the crust [10].

A series of moderate earthquakes have repeatedly struck this area in the last 400 years. The largest one occurred near Norcia town in 1703, with a magnitude of 6.8. The closest large historical event was the 1639 Mw 6.2 earthquake that took place near Amatrice town. According to instrumental records, two Mw > 6 earthquakes have also struck the area in modern times. One was the 1997 Mw 6.0 Colfiorito earthquake that occurred to the northwest, and the other was the 2009 Mw 6.3 L'Aquila earthquake to the south. The 2016

earthquake sequence occurred in a seismic gap which is located between the areas hit by the 1997 Colfiorito earthquake and the 2009 L'Aquila earthquake.

## **3. Rupture Model**

*3.1. Data*

We used four SAR image pairs to measure the ground displacement relative to the 24 August M w 6.1 Amatrice earthquake, the 26 October M w 5.9 Visso earthquake, and the 30 October M w 6.5 Norcia earthquake (Figure A1 and Table A1). To better serve the rupture model inversion, each selected SAR pair only covered one event mentioned above. Among the image pairs, three were from the Sentinel-1 satellite and the other was from the Advanced Land Observing Satellite (ALOS) satellite. Many institutions and researchers have created interferograms for the Amatrice and Visso events using Sentinel-1 images, and their results are quite similar. In this study, we directly used the InSAR interferograms from the European Space Agency's InSARap program for the analysis (provided free online at http://insarap.org/, last accessed 20 March 2022). For the 24 August event, the deformation pattern is characterized by two NNW–SSE striking main distinctive lobes with different shapes, but having almost the same maximum negative line-of-sight (LOS) deformation of around 20 cm. Regarding the 26 October event, the interferogram reveals an ear-shaped deformation pattern, striking NNW–SSE, similar to the 24 August event.

The 30 October M w 6.5 Norcia earthquake generated much larger ground displacement than the previous two events, making the InSAR measurements more challenging. Therefore, we adopted the L-band ALOS 2 SAR image pair, which has better resistance to phase incoherence, to retrieve the coseismic deformation. A two-pass technique was used to process the data with Gamma software. The TanDEM-X DEM with a 12-m resolution was used to remove the topographic components in the interferogram. A baseline refinement step was carried out to remove the ionospheric disturbance. Finally, we obtained the displacement map after unwrapping the interferogram by the use of the minimum cost flow (MCF) algorithm. The InSAR interferogram reveals a similar ear-shaped deformation pattern to the 26 October event, with around ~90 cm maximum negative LOS displacement. To reduce the quantity of InSAR data, we first applied uniform downsampling to reduce the displacement map to a size of ~500 × 500, and we then used the equation-based quadtree downsampling algorithm to select ~1000 LOS measurements from each interferogram [23].

In addition to the InSAR datasets, we also collected three-component strong-motion datasets of the three earthquakes for joint inversion [24]. These data were obtained from the Italian National Accelerometric Network operated by the Istituto Nazionale di Geofisica e Vulcanologia (INGV) and the Rete Accelerometrica Nazionale (RAN). We selected strongmotion stations with good azimuthal coverage and epicentral distances of less than 50 km. Under these criteria, we separately adopted strong-motion data from 30, 24, and 19 stations for the 24 August, 26 October, and 30 October earthquakes (Figure A2). To better represent the longer period features of strong motion, the velocity waveform data were used for the inversion. We first removed the mean offset and instrument response from the original accelerogram, then filtered the local site effects, and finally, integrated the accelerogram in time. The frequency of each time-series velocity waveform was resampled to 2 Hz to reduce the computational burden during the joint inversion. More details on the inversion strategy and setting are provided in the Supporting Information Texts S1 and S2 [9,25–31].

#### *3.2. Slip Models*

We tested three different inversion scenarios: InSAR data inverted alone, strongmotion data inverted alone, and a joint inversion with both geodetic and seismic data (Figures A2 and A3). In the three events, the InSAR-only model differs from the strongmotion data-only model in slip distribution. However, the joint inverted slip distributions seem to make a compromise, absorbing the characteristics from both models. Overall, the joint slip model is closer to the InSAR-only results, which is consistent with the fact that the near-field InSAR data have a better ability to constrain the fault slip pattern. The InSAR

and strong-motion prediction from the joint model fits quite well with the observations (Figures A2 and A7, Figures A8 and A9).

The joint model of the 24 August earthquake shows two separate major slip concentrations with a maximum slip of 0.76 m and 0.72 m, respectively, locating at depths between 5 km and 3 km (Figures 2 and A4). The characteristics of the rupture model are in general accordance with previous models [13,14], exhibiting a normal faulting mechanism, bilateral rupture directivity, and a relatively fast rupture velocity. A more detailed comparison with previous models is provided in supporting information Text S3 [9–14,16,32]. Assuming a shear modulus of 30 GPa, the overall seismic moment of the two fault segments is 1.6 × 10<sup>18</sup> Nm, equivalent to a moment magnitude of M w 6.1. The slip pattern is mostly constrained by the near-field InSAR data, but the relative far-field strong motion still fits quite well. During the first 6 s, the majority of the seismic moment was released. The moment rate rapidly increased in the initial stage and reached 3.9 × 10 17 Nm/s at 2.8 s, and then decreased rapidly with time. The rupture process took place in a relative manner, propagating gradually from the epicenter to distant patches, and no delayed slip patches were observed.

**Figure 2.** The inverted joint slip models for the 2016 central Italy earthquake sequence. (**a**) Distributions of slip at depth for the Visso (left) and Amatrice (right) earthquakes. The yellow stars denote the start point for the ruptures. (**b**) Distributions of slip at depth for the Norcia earthquake. (**<sup>c</sup>**–**<sup>e</sup>**) The moment rate functions of the Visso, Amatrice, and Norcia earthquakes, respectively.

For the 26 October earthquake, our joint model suggests a single elongated normal faulting slip concentration in the northern segmen<sup>t</sup> of the VBFS (Figures 2 and A5). The rupture started from the southeast part of the fault plane and propagated mostly unilaterally toward the northwest. The moment rate rapidly increased at the initial stage, reaching 2 × 10 17 Nm/s at 2.1 s, but the relatively high moment rate (>1 × 10<sup>17</sup> Nm/s) lasted for about 4 s, and then decreased rapidly. This event released a seismic moment of 1.5 × 10<sup>18</sup> Nm, corresponding to a Mw 6.1 event.

The 30 October earthquake was the largest event in the earthquake sequence. The inverted rupture model suggests a normal faulting mechanism, consistent with the moment tensor solution, as well as the long-term behavior of the VBFS system. Two slip patches were found, the larger one located in the north with a maximum slip of 2.8 m at a depth between 4 km and 6 km, and the other one peaking at 2.3 m at a depth between 5 km and 7 km (Figures 2 and A6). It can be noted that the latter slip patch is almost below the northwest slip concentration of the 26 August event. The slip history shows that the largest moment rate reached 1.8 × 10<sup>18</sup> Nm/s at 4.2 s, and almost all of the coseismic moment was released in the first 7 s. The overall seismic moment was 1.1×10 19 Nm, equivalent to a Mw 6.7 event.

#### **4. Mapping of the** *b***-Values and Scarp Offsets**

## *4.1. b-Values*

The Gutenberg–Richter (G–R) law is the commonly used statistical model when describing the size distribution of earthquakes. In G–R law, the number (N) of earthquakes having a magnitude ≥M follows a logarithmic relationship with magnitude M, expressed as log10 N = a − bM, where a and b are constants. Parameter b describes the occurrence ratio of small to large earthquakes. The variance of the *b*-value is thought to be related to local conditions, e.g., stress applied to the material, the strength heterogeneity of the material, the crack density, and the thermal gradient [33]. Among these potential factors, applied stress is usually cited. It is believed to have a negative linear relationship with the *b*-value, which has been proved by a number of laboratory experiments and earthquake observations. In this section, we investigate the *b*-value variation of the 2016 central Italy earthquake sequence area, attempting to reveal the potential stress heterogeneity of the seismogenic fault. ZMAP software [34] was used to calculate the *b*-values. An event catalog covering the 12 recent years of 3 April 2005–8 October 2017 with a depth below 30 km was downloaded from the INGV (http://cnt.rm.ingv.it/en, last accessed 20 March 2022) and adopted as the data source (Figure 3a).

The results give a cut-off magnitude (MC) for this region of 1.3. The average *b*-value for this sequence is estimated to be equal to 1.03 ± 0.02, with a 90% goodness-of-fit level. The average *b*-value is slightly larger than the global mean value of 1.0, which is consistent with previous findings of normal fault-related earthquakes having relatively high *b*-values compared with strike-slip and reverse-slip faulting mechanisms [35]. We further mapped the spatial distribution of the *b*-value in a 0.05◦ × 0.05◦ grid (Figure 3a). Each grid selected 300 neighboring events and required at least 50 events above the local value of MC. The estimated *b*-values vary from 0.55 to 2.05. The three events are found to be located in the low *b*-value region, suggesting high-stress regimes. High *b*-values are also observed to the southwest of the 24 August earthquake hypocenter. This event happened in the connecting area between the Gorzano fault and Vettore fault, where the Olevano–Antrodoco–Sibillini (OAS) thrusting structure is thought to intersect. The complex structure may generate a highly fractured rock mass, which is thought to correspond to large *b*-values.

**Figure 3.** *b*-value and cumulative vertical displacements for the 2016 central Italy earthquake sequence. (**a**) Seismic activity, *b*-value distribution map, and cross section using the INGV earthquake catalog covering the period from 3 April 2005 to 8 October 2017. The pink dashed rectangles indicate the location of the *b*-value swath profile. The red, blue, and black dashed lines indicate the 60% maximum slip area of the Amatrice, the Visso, and the Norcia earthquakes. (**b**) Cumulative vertical displacements (red arrows) along the VBFS/GFS fault systems, and an example of scarp offset measurement. Red arrows represent vertical offsets of ruptures, and they are orthogonal to fault strikes.

We also mapped the depth variation of the *b*-values across the rupture zone (Figure 3a), which was performed in a 1 km × 1 km grid on the swath profile, with 300 neighboring events and a minimum of 50 events above MC. The cross section clearly reveals a low *b*-value layer above 10 km, suggesting that high differential stress exists in the first 10 km of the upper crust. The lower part has low stress possibly due to the presence of fluids in the rock matrix, and the energy is released by a series of small events, as observed. However, these features are not held to the southeastern end of the GFS, where the number of fault branches increases followed by rotating counterclockwise from NNW–SSE to NWW–SEE trending [36], indicating that the magnitude and direction of tectonic stress changed considerably.

#### *4.2. Scarp Offsets*

Scarps are universal along the GFS and VBFS systems. Measuring the scarp offset through topography analysis is a common approach in paleoseismology for slip rate estimation. The diversity of the offset values along the fault trace can, therefore, represent the rupture history on the fault. We adopted the 12 m resolution TanDEM-X DEM obtained before the earthquake sequence to extract the scarp offsets along the seismogenic fault trace. Guided by the active fault map reported in Falcucci et al. (2016) [37], we identified the fault trace through visual analysis of high-resolution Google imagery. We then exacted the elevation profile from the TanDEM-X DEM across the fault trace. Two parallel lines orthogonal to the fault strike were generated by least-squares fitting on each side of the fault. As the height profile across the fault is complex, we designed an interactive approach where the user can manually select a profile section with a relatively constant slope for linear fitting. The distance between the lines is considered as the vertical offset of the scarp. As the long-term interseismic slip may not generate a near-fault scarp, we assume that the measured topographic offsets were contributed by large historic earthquakes.

We located 28 sites in the Google imagery and TanDEM-X DEM where obvious scarp features can be observed (Figure 3b and Table A2). The scarp features distribute along the surface trace of the Norcia and Visso earthquakes and to the south of the Amatrice earthquake, while no obvious scarp is found on the south fault segmen<sup>t</sup> of the Amatrice earthquake. This agrees with the results reported in Falcucci et al. (2016) [37], i.e., there is no evidence at the surface of late Quaternary fault activity in this area. The scarp features are consistent with the mechanisms of the 2016 central Italy earthquake sequence. Vertical offsets in the height profile are obvious at these sites, while no horizontal displacement is visible in the high-resolution satellite imagery. The average of the measured offsets is about 1.4 m, with the largest value exceeding 3.8 m, which is approximately 1.5 times larger than the maximum slip found for the 30 October Norcia earthquake.

## **5. Discussion**

#### *5.1. Fault Behavior*

Historical earthquakes in the central Apennines of Italy show an obvious space-time clustering behavior, where one main shock triggers a series of subsequent events in a relatively short period. Two recent earthquakes nearby this sequence share the same cascading character. To the north, the 1997 M w 6.0 Colfiorito earthquake triggered six M > 5 events in 20 days. To the south, the 2009 M w 6.3 L'Aquila earthquake started a strong sequence of aftershocks. The 2016 sequence is a typical cascading event in the central Apennines, which offers us another opportunity to investigate the fault behavior in this region.

One notable feature in our joint slip models is the complex multi-fault segments and heterogeneous coseismic slip distribution. Several slip concentrations with different maximum slips are located in different fault patches. Such heterogeneity is also evidenced by spatial *b*-value mapping, where the estimated *b*-value varies largely in different locations (Figure 3a). The slip distributions of three events connect quite well, with few overlaps or gaps, implying an almost complete release of the accumulated interseismic energy. An

exception is the conjunction location between the two fault segments of the 26 August event. Such discontinuity is thought to be related to the inherited compressional thrust fault. The interaction between the inherited thrust fault and the active normal fault controls the seismicity very obviously, by which the aftershocks are divided into two clusters, separately distributing on each side of the inherited thrust fault. Fault segmentation is a common features found in earthquakes worldwide, where the fault is composed of discrete segments divided by geometrical discontinuities [38]. It normally acts as a structural control on earthquake magnitude and rupture progress. The central Apennines is dominated by a mix of modern and ancient structures, resulting in heterogeneous crust and segmented faults. Comparing some mature faults which have up to 10 million years of history (e.g., most faults in the San Andreas fault system), the modern extending faults in the central Apennines starting from about a half million years ago are still very young. In an immature fault system, it is mechanically difficult for rupture propagate from one segmen<sup>t</sup> to another. The moderate magnitude and slip heterogeneity as observed in our joint slip model is thus a result of the segmentation and immaturity of the seismogenic fault system. Research on the 2016 Kaikoura earthquake [19] warns that multi-fault segments are likely to rupture simultaneously to generate a large earthquake, but its seismogenic faults within the southern island of New Zealand are relatively mature. The immature fragmented fault systems in the central Apennines may still favor moderate magnitude earthquakes at the current stage.

Normally, it is believed that the heterogeneous fault slip distribution is mainly controlled by the fault strength variation [39]. As reported by the joint slip model, several asperities were found in this fault system, with the middle one corresponding to a much larger fault strength. From the observation of the *b*-value mapping, we can also find that low *b*-values are associated with these strong patches. The *b*-value represents the relative occurrence of large and small events, where a low value indicates a larger proportion of large earthquakes [35]. The *b*-value is normally believed to be related to the stress condition, as large earthquakes are thought to be caused by highly loaded stress. Beneath the three mainshock hypocenters, we can observe a high *b*-value region associated with low strain energy (Figure 3a). It is probably due to the presence of fluids in the rock matrix, which is also proposed to exist in the nearby 2009 L'Aquila earthquake [40]. The hypothesis of deep fluid intrusion may also give a good explanation for the space-time clustering behavior of earthquakes in the central Apennines [41]. Pore fluid pressure diffusion after an earthquake can dramatically reduce the shear strength of adjacent fault segments, and consequently induce a cascade of multiple events on them even they are not fully loaded.

Another notable feature of this earthquake sequence is the rupture directivity of the slip models. The retrieved rupture progress shows that the nucleation of three major events all started at the bottom of the fault and then propagated to the upper patches. The along-strike unilateral propagation is also very obvious. The Amatrice and Visso events propagated mostly toward the NNW, and the Norcia event was toward the SSE. The preference for unilateral propagation is observed in many earthquakes, and a potential explanation is the fault segmentation [42], as the earthquake prefers to propagate unilaterally along strike until reach discontinuities. The up-dip rupture direction might result from the material property contrast along with the depth. As denoted above, high *b*-value is observed at the deeper depth from the *b*-value cross section, suggesting the existence of fractured rock mass saturated with fluids over there. The rupture may thus prefer to initiate at the more compliant down-dip part of the rupture zone.

Laboratory experiments have shown that non-uniform normal stress characterized by a high amplitude single-stress asperity favors the occurrence of strong characteristic microquakes, which share similar locations, magnitudes, and return periods [43]. Both rupture model and *b*-value mapping sugges<sup>t</sup> heterogeneous stress in terms of several asperities in this fault system, and the middle one associated with the M w 6.7 Norcia earthquake has the largest amplitude. We can expect that the M w 6.7 Norcia earthquake was a characteristic earthquake in this fault system, and that it may rupture again within a

certain period. An effective way to test the long-term characteristic behavior of an active fault is through a comparison between historical fault scarps and recent earthquakes [44]. We can refer to such independent observations to verify this characteristic earthquake assumption. In this case, we obtained the cumulative scarp features of historic earthquakes through measuring the near-fault topographic offset. The scarp characteristics agree quite well with the joint slip model in this fault system, where the vertical displacement is obvious in the seismogenic fault of the Norcia event and is several times larger than the coseismic slip. This suggests that long-term constant strong patches locate in this area, and this fault may have displayed characteristic slip behavior.

#### *5.2. Slip Budget Closure Test*

If we assume a characteristic fault behavior in this fault system for this earthquake sequence, an important issue related to the seismic hazard assessment is whether the Norcia earthquake is the largest earthquake in this fault and whether there is any unruptured fault segmen<sup>t</sup> that has the potential to generate a larger earthquake. For a characteristic earthquake, it is possible to derive information about the maximum-magnitude earthquake and its return period from the earthquake catalog by assuming that the frequency magnitude follows G–R law. The relationship between the maximum magnitude ( M*w*) and its corresponding frequency (N*max*) can be expressed as [45]:

$$\log \mathbf{N}\_{\text{max}} = -\frac{3}{2} \mathbf{M}\_{\text{w}} - 9 + \log \mathbf{M}\_0 + \log \left[ a \left( 1 - \frac{2b}{3} \right) \right] \tag{1}$$

where, *α* is the fraction of transient slip that is seismic, *b* is the *b*-value in G–R law, and M0 is the seismic moment deficit.

The slip budget closure test can therefore be run with the maximum earthquake frequency equation and G–R law. This test has been found to be quite successful in several faulting systems, such as the longitudinal valley fault in Taiwan and the Sumatra Megathrust fault [45]. It is thus interesting to examine the fault slip budget in this complex normal faulting system under frequency–magnitude law.

We tested the seismicity in the area between the 1997 Colfiorito and 2009 L'Aquila earthquakes. We assumed a 75 km length (L) and 16 km width (W) fault plane and half-area (S*L* = 0.5×L × W) is locked according to the observed distribution of the earthquakes. The long-term slip rate (V) was set to be 2 mm/yr, referring to previous studies [20,21]. Using a shear modulus (G) of 30 GPa, we obtained a rough seismic moment deficit of 7.06 × 10 16 Nm/yr, according to the simplified equation M0 = *G*S*LW*.

We combined three catalogs for the frequency–magnitude linear fitting. The INGV catalog includes a wide range of earthquakes with different magnitudes, but the span period is only from 2005 to 2022. The USGS catalog has a longer cover period from 1950 to 2022, but is incomplete for small events. The DBMI catalog contains historical events between 1005 to 2014, but also lacks small and recent earthquakes (http://emidius.mi.ingv.it/DBMI11/, last accessed 20 March 2022). As earthquakes in the central Italy often occurred in terms of cascading earthquake storm, the earthquake frequency in a short-period earthquake catalog is significantly biased. In order to obtain a more comprehensive catalog, we formed a combined DBMI-USGS catalog, where earthquakes before 2014 are from DBMI and events between 2014 and 2017 are from the USGS. The three catalogs give a similar *b*value, while the *a*-values differ a lot. This is because the seismic rates are overestimated in the previous two catalogs. Therefore, we prefer to use the combined DBMI-USGS catalog. If we assume that the seismic moments are released partly seismically with *α* = 0.8, estimated according to the area proportion of asperities on fault planes, the return period line given by Equation (1) will intersect with the frequency–magnitude line at point (M w = 6.7, Lg N = –3.02), suggesting a maximum magnitude of 6.7 and a predicted return period of 1064 years (Figures 4 and A10).

**Figure 4.** Seismic slip budget closure test on the 2016 central Italy earthquake sequence. The red line represents the G–R law fitted by the DBMI-USGS catalog. The blue line shows the return period of the maximum magnitude event given in Equation (1). The yellow star depicts the intersection of the blue line and red line, giving a maximum magnitude of 6.7 and a return period of 1064 years.

We can also directly estimate the earthquake recurrence interval T by T = S/V, as suggested by Shen et al. (2009) [46], where S is the mean coseismic slip on the fault segmen<sup>t</sup> and V is the secular slip rate. Our joint slip model suggests 1.8 m average slip on the slip concentration for the Norcia event (~0.16-m and ~0.24-m average slips for the Amatrice and Visso events). Together with the predetermined ~2 mm/yr slip rate [20,21], we can obtain a recurrence interval of around 900 years for the Norcia fault segmen<sup>t</sup> (~80- and ~120-year recurrence intervals for the Amatrice and Visso fault segments). This is in general agreemen<sup>t</sup> with the slip budget closure test, implying a low seismic hazard for this area in the future. However, we should also note that both results are rough estimations because they rely on many assumptions, such as the characteristic earthquake, the frequency–magnitude law, the secular fault slip rate, and the catalog completeness. Meanwhile, although the return period for a maximum-magnitude earthquake is long (~1000 yr), such normal faulting areas usually have relatively high *b*-values, which means that the return period of an earthquake drops sharply with the decrease of the magnitude, according to the magnitude–frequency law. It should also be noted that the risk of smaller earthquakes (Mw 5–6) still exists, because they can rupture again with a much shorter return period.

## **6. Conclusions**

This paper has provided a complete rupture history of the three main shocks in the 2016 central Italy earthquake sequence. The 24 August Amatrice earthquake occurred on two fault segments divided by an inherited compressional thrust fault, with a total seismic moment of 1.6 × 10<sup>18</sup> Nm (Mw 6.1). The 26 October and 30 October events both ruptured in the VBFS system, releasing seismic moments of 1.5 × 10<sup>18</sup> Nm (Mw 6.1) and 1.1 × 10<sup>19</sup> Nm (Mw 6.7), respectively. The *b*-value mapping reveals a complex and non-uniform stress condition in this area. The complex segmented faults and heterogeneous coseismic slip distribution sugges<sup>t</sup> an immature fault behavior over there. Possible fluid intrusion implied by high *b*-value at deep depth may be the cause of the rupture of multi-fault segments in a short period. The cumulative surface scarp displacement from past earthquakes shows faulting that is consistent with the 2016 earthquake sequence, suggesting a characteristic fault behavior. Under the slip budget closure test for characteristic earthquakes, we obtained a maximum magnitude of 6.7 in this area. It is likely that the maximum earthquake has already been triggered in this sequence, and the seismic hazard from large earthquakes (M 6+) in this region will be at a low level for a long period (~1000 years).

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/rs14081819/s1, Text S1: inversion strategy; Text S2: inversion setting; Text S3: slip model comparison.

**Author Contributions:** Conceptualization, C.W.; methodology, C.Z., C.W. and G.Z.; software, C.Z. and C.W.; validation, X.S., Q.L. and J.Z.; formal analysis, B.Z. and P.L.; writing—original draft preparation, C.Z. and C.W.; writing—review and editing, C.Z., C.W., B.Z. and P.L.; supervision, Q.L. and J.Z.; funding acquisition, C.Z., C.W. and Q.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Key Research and Development Program of China (Grant No. 2019YFC1509203), the National Natural Science Foundation of China (Grant No. 41974006), the Shenzhen Science and Technology Program (Grant No. RCYX20210706092140076 and 20200812164904001), the state Key Laboratory of Earthquake Dynamics (Grant No. LED2021B05), and the Lhasa National Geophysical Observation and Research Station (Grant No. NORSLS20-08).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The S1A InSAR interferograms were obtained from the European Space Agency's InSARap program (http://insarap.org/, last accessed 20 March 2022). ALOS-2 data are copyright of JAXA (2016), and they have been collected through the RA6 project (ID: 3078). The strong-motion data were obtained from the INGV and the RAN (available from the Engineering Strong-Motion database website: http://esm.mi.ingv.it/, last accessed 20 March 2022).

**Acknowledgments:** The authors would like to thank the European Space Agency, JAXA (2016), the INGV and the RAN for providing SAR images and strong-motion data. The Multi-Data Source Modeling and Inversion Toolkit was used for the strong-motion data processing and mapping (https://github.com/dmelgarm/MudPy, last accessed 20 March 2022). Some of the figures were generated by Generic Mapping Tools (GMT 4.5 from http://gmt.soest.hawaii.edu/, last accessed 20 March 2022). We are grateful to the three anonymous reviewers for constructive comments, that helped us to improve this manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

**Appendix**

 **A**

**Figure A1.** SAR Interferograms, predictions, and residuals based on the joint slip models derived in this study for the Central Italy earthquake sequence. The S1A interferograms are wrapped to 2.8 cm, and the ALOS interferogram is wrapped to 11.45 cm. The red beach ball in each interferogram represents the location and focal mechanism of the corresponding event.

**Figure A2.** Distribution of the strong-motion stations used to retrieve the slip models of the Amatrice (**a**), Visso (**b**), and Norcia (**c**) earthquakes.

**Figure A3.** Comparison of slip models using different data constraints. (**a**) The Amatrice earthquake. (**b**) The Visso earthquake. (**c**) The Norcia earthquake.

**Figure A4.** Rupture snapshots in a time interval of 1 s for the Amatrice earthquake.

**Figure A5.** Rupture snapshots in a time interval of 1 s for the Visso earthquake.

**Figure A6.** Rupture snapshots in a time interval of 1 s for the Norcia earthquake.

**Figure A7.** Comparison of strong-motion velocity records (black) and the synthetic seismograms (red) predicted by the joint slip model of the Amatrice earthquake.

**Figure A8.** Comparison of strong-motion velocity records (black) and the synthetic seismograms (red) predicted by the joint slip model of the Visso earthquake.

**Figure A9.** Comparison of strong-motion velocity records (black) and the synthetic seismograms (red) predicted by the joint slip model of the Norcia earthquake.

**Figure A10.** Distribution of seismic activity along the VBFS/GFS systems from the INGV, USGS and DBMI-USGS catalogs. The dashed black rectangles represent the surface projection of the fault planes adopted in this study. The bold red lines are the two seismogenic normal fault systems, namely, the Mt. Vet-tore–Mt. Bove Fault (VBFS) system and the Mt. Gorzano Fault (GFS) system.

**Table A1.** SAR scenes used for generating the coseismic interferograms.



**Table A2.** Scarp offset measurements through topography analysis.
