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

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 × 1018 Nm (Mw 6.1), 1.5 × 1018 Nm (Mw 6.1), and 1.1 × 1019 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.


Introduction
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][2][3][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][6][7][8].
On 24 August 2016 (UTC 01:36, local time 03:36), a destructive earthquake (M w 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 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 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]. Ã 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 M w 6.2 earthquake that took place near Amatrice town. According to instrumental records, two M w > 6 earthquakes have also struck the area in modern times. One was the 1997 M w 6.0 Colfiorito earthquake that occurred to the northwest, and the other was the 2009 M w 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.

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][26][27][28][29][30][31].

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][10][11][12][13][14]16,32]. Assuming a shear modulus of 30 GPa, the overall seismic moment of the two fault segments is 1.6 × 10 18 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. and strong-motion prediction from the joint model fits quite well with the observ ( Figures A2 and A7-A9).
The joint model of the 24 August earthquake shows two separate major slip c trations with a maximum slip of 0.76 m and 0.72 m, respectively, locating at dep tween 5 km and 3 km (Figures 2 and A4). The characteristics of the rupture mode general accordance with previous models [13,14], exhibiting a normal faulting m nism, bilateral rupture directivity, and a relatively fast rupture velocity. A more d comparison with previous models is provided in supporting information Text S3. A ing a shear modulus of 30 GPa, the overall seismic moment of the two fault segm 1.6 × 10 18 Nm, equivalent to a moment magnitude of Mw 6.1. The slip pattern is constrained by the near-field InSAR data, but the relative far-field strong motion s quite well. During the first 6 s, the majority of the seismic moment was released. T ment rate rapidly increased in the initial stage and reached 3.9 × 10 17 Nm/s at 2.8 then decreased rapidly with time. The rupture process took place in a relative m propagating gradually from the epicenter to distant patches, and no delayed slip p were observed.  For the 26 October earthquake, our joint model suggests a single elongated normal faulting slip concentration in the northern segment 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 17 Nm/s) lasted for about 4 s, and then decreased rapidly. This event released a seismic moment of 1.5 × 10 18 Nm, corresponding to a M w 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 18 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 M w 6.7 event.

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 log 10 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.   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.

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 segment 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.

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 segment 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 suggest 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.

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 segment 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]: where, α is the fraction of transient slip that is seismic, b is the b-value in G-R law, and M 0 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 M 0 = GS L W.
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 bvalue, 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).
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 (Mw = 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). We can also directly estimate the earthquake recurrence interval T by T = S/V, as suggested by Shen et al. (2009) [38], where S is the mean coseismic slip on the fault segment 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 segment (~80and ~120-year recurrence intervals for the Amatrice and Visso fault segments). This is in general agreement 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.

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 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 segment 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 segment (~80-and~120-year recurrence intervals for the Amatrice and Visso fault segments). This is in general agreement 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 (M w 5-6) still exists, because they can rupture again with a much shorter return period.

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 18 Nm (M w 6.1). The 26 October and 30 October events both ruptured in the VBFS system, releasing seismic moments of 1.5 × 10 18 Nm (M w 6.1) and 1.1 × 10 19 Nm (M w 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 suggest 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).
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.               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 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 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 A9. Comparison of strong-motion velocity records (black) and the synthetic seismograms