This electronic supplement contains earthquakes location and catalog, finite-fault analysis, figures of station location, velocities and VP/VS ratios, comparison of key parameters, time-dependent estimation of the completeness magnitude, and waveform fits.
We relocated 32,771 earthquakes detected in the period from 24 August 2016 to 29 November 2016 by the National Seismic Network (RSN) managed by the National Institute of Geophysics and Volcanology (INGV) by inverting P- and S-wave arrival times handpicked by the seismologists on duty in the seismic monitoring room of the INGV (Moretti et al., 2016). This analysis has been performed in quasi-real-time thanks to standard webservices giving prompt direct access to the real-time database. We added 12 stations to the RSN, installed in the epicentral area soon after the 24 August 2016 Mw 6.0 event by the INGV group handling the temporary and mobile network and connected in Universal Mobile Telecommunication System (UMTS) real time to the acquisition center, to improve both detection capacity and monitoring quality mainly in terms of aftershocks’ hypocentral depth (Fig. S1). Data from the Accelerometric National Network (RAN) stations have been added only for locating the major events (Fig. S1).
Hypocenter locations have been estimated with a 1D velocity model with a robust nonlinear inversion code (NonLinLoc [NLL]; Lomax et al., 2009), providing a comprehensive description of the location uncertainties. In addition, NLL separates the spatial and temporal location of earthquakes, reducing the trade-off between depth and origin time. The 1D P- and S-wave velocity model (Fig. S2) is a gradient version of the layered minimum 1D model estimated for the region by Carannante et al. (2013). In such a complex geological environment, we prefer a gradient velocity model being able to avoid possible clustering of hypocenters around sharp discontinuities. To further reduce systematic delay times due to the use of a 1D model, we used more than 300,000 P- and S-arrival times to calculate travel-time corrections for the (≈120) used stations. This correction substantially increased the location quality for the vast majority of the aftershocks but not for all of the higher-magnitude earthquakes (Mw ≥5.0) necessarily constrained also by very distant stations introducing instabilities due the strong 3D heterogeneities. Only for these events, we added an additional constraint by fixing the VP/VS ratio equal to 1.86, a value we derived for the area with the Wadati (1931) method (Fig. S3b).
The final catalog (available in Table S1), containing 25,620 earthquakes with 0.1 ≤ M ≤ 6.5, has been obtained by applying the quality selection criteria reported in Table S2. Because of the higher root mean square (rms) characterizing the majority of the highest magnitude earthquakes, a less severe selection criteria have been applied to the M >4.0 events, to keep them within the catalog for the sake of information completeness. The M >4.0 events have been reported in Table S3 with related errors and focal mechanism solutions when available.
To build the space–time chart of Figure 2 in the main article from the beginning of 2016, we also relocated with the same criteria the about 1500 events that occurred in an area of 50 km of radius from the first Mw 6.0 mainshock. Thus, the events in the figure are the ones detected by the RSN in the period from 1 January 2016 to 23 August 2016.
In Figure S4, we report the distribution of the estimated uncertainties related to longitude, latitude, depth, rms, and gap of the final dataset (red lines) versus the whole set of the relocated events. This comparison shows that the selected events are representative of the whole set of data. However, about 90% of the catalog has extremely high-quality locations with horizontal errors ≤ 0.5 km, vertical error ≤ 1.5 km, rms ≤ 0.3 s, and azimuthal gap ≤ 120° substantially unchanged during the seismic network deployment evolution. The highest azimuthal gaps, between 120° and 180°, obviously pertain to the edges of the network, especially to the western zone more poorly covered, that will benefit from the data of not real-time stand-alone stations.
Finally, we estimate the time-dependent completeness magnitude (Mc) of the aftershock catalog (Fig. S5), by assuming the Gutenberg–Richter law and by applying two different methods, all based on the maximum-likelihood estimation of b-value. The goodness-of-fit method (GFT; Woessner and Wiemer, 2005) measures the probability (prMc) that the Gutenberg–Richter law explains the data. We consider the value giving the largest value of prMc as the GFT Mc estimation, whereas the magnitude and b-value stability method (MBS; Woessner and Wiemer, 2005) measures the stability of the estimated b-values for consecutive magnitude bins (Nbin). In this study, we fix Nbin equal to 5. To balance the temporal evolution of seismicity with the number of events, we apply the GFT and MBS methods in overlapping time windows of three days, daily updated (at 00:00:00).
The two methods generally give consistent results for most of the interval times. Some differences are observed for the time windows starting 17, 18 September, 7, 8, 13, 14 October, and 23 November, due to some anomalies in the provisional data. Beside that, we observe the usual decrease in the completeness coinciding with the occurrence of the largest event and subsequent increase of seismic rate.
The inversion code is based on the method of Hartzell and Heaton (1983), as implemented by Dreger and Kaverina (2000) and Dreger et al. (2005). It consists of a nonnegative, least-squares inversion method of observed velocity seismograms to compute the finite-source kinematic models. Several constraining equations are applied to improve the stability of the inversion. These include the requirement that slip is positive (the fault cannot slip backward) and minimization of the spatial derivative of slip to smooth the results. The strength of the smoothing constraint is controlled by a weight, which is determined iteratively by searching for smoothed models that produce an acceptable level of fit to the waveform data, with simultaneous smoothing and damping. This approach also assumes constant rupture velocity and rise time but allows the use of multiple time windows to account for potential variations in rupture speed and local rise time. The slip velocity is modeled by imposing a simple boxcar source-time function. The best-fitting slip duration and rupture velocity in each time window is selected iteratively by performing inversions with different values of these parameters and quantitatively measuring the fit between synthetic and observed waveforms, which is quantified through the variance reduction parameter (VR)
in which i is the station index, xi(t) the synthetic waveform, di(t) the recorded waveform, and wi is the inverse epicentral distance weight, that is, the more distant stations have a larger weight, a weight proportional to the epicentral distance.
Rake parameter can be heterogeneous on the fault plane or assumed constant all over the fault and chosen iteratively as is the rupture velocity.
Synthetic seismograms are represented as a linear combination of eight fundamental Green’s functions that form the basis functions for any arbitrarily oriented double-couple (DC) mechanism. The Green’s functions were computed using a FORTRAN, frequency–wavenumber integration program (Saikia, 1994) on a regular grid sampling the focal volume every 1 km horizontally and 1 km vertically and filtered between 0.02 and 0.5 Hz, the same as for the recorded data. In this study, we adopt precalculated and stored Green’s functions obtained using the Central Italian Apennines (CIA) velocity model (Herrmann et al., 2011). Although the adopted methodology has very simplistic assumptions, it has been found generally effective to retrieve the gross slip distribution and to image the main features of the earthquake rupture history in terms of heterogeneous distribution of slip amplitudes, slip direction, and average rupture velocity.
We used waveforms recorded by the three-component digital accelerometers of the RAN (Rete Accelerometrica Nazionale), managed by the Italian Department of Civil Protection, and INGV (Michelini et al., 2016) networks (Fig. S1), with epicentral distances less than ~45 km. Prior to the inversion procedure, the recorded accelerograms were processed to remove the mean offset and instrument response, band-pass filtered between 0.02 and 0.5 Hz using a low-pass filter with three poles and a high-pass filter with three poles, and finally integrated in time to obtain ground velocities. The relatively low maximum frequency selected is required to reduce the possible contributions of local site effects (Bindi et al., 2011) to the source inversion modeling. The processed time histories were then resampled at 10 samples per second. Because the inverted stations have different elevations ranging between ~560 and ~1340 m.a.s.l. (meters above sea level), we determined the average elevation (~850 m) to which the finite-fault solutions are referred.
The 26 October Mw 5.9 earthquake is actually a double event, clearly visible at station CNE that is located 2 km from the epicenter (Fig. S6). This source complexity mainly affects the waveform fit of the closest stations CMI, CNE, and NRC, where we better fit the phase rather than the amplitude of the horizontal components that result underestimated. However, the retrieved source model shown in Figure 4b in the main article shows a relatively good variance reduction (47.5%), and it allows a good fit to all the other 33 stations, especially those located toward the northeast (Fig. S7).
Figure S8a,b shows the waveforms fit obtained for the 30 October Mw 6.5 earthquake. The synthetic seismograms match quite well the 39 recorded ground velocities, corroborating the goodness of the retrieved kinematic fault model (Fig. 4c in the main article). The obtained variance reduction (40%) is good, considering the high number of inverted stations. Some inaccuracies are present in the theoretical seismograms calculated for the southernmost stations. These features could suggest an uncalibrated slip model for the southernmost portion of the inverted fault that might be related to a geometrical source complexity. We will deeply explore this feature in the near future.
The details of the kinematic models for the Mw 5.9 and 6.5 events are available in Tables S4 and S5, respectively.
Table S1 [Plain Text Comma-separated Values; ~2.2 MB]. Location parameters for the 25,620 earthquakes included in this study. eqkID (earthquake id), OT (origin time, yyyy-mm-ddThh:mm:ss.sss), latitude (°N), longitude (°E), depth (km), errh (km) (error of the horizontal component), errv (km) (error of the vertical component), gap (°) (stations azimuthal gap), rms (s) (root mean square), nphs (number of phases used for the location) and M (magnitude; usually ML and MW for smaller and larger events, respectively).
Table S2. Quality selection criteria.
Table S3 [Plain Text Comma-separated Values; ~7 KB] Location parameters and focal mechanisms solutions for the larger (ML ≥4.0) earthquakes. eqkID (earthquake id), OT (origin time, yyyy-mm-ddThh:mm:ss.sss), latitude (°N), longitude (°E), depth (km), errh (km) (error of the horizontal component), errv (km) (error of the vertical component), gap (°) (stations azimuthal gap), rms (s) (root mean square), nphs (number of phases used for the location), ML (local magnitude), Mw (moment magnitude). S1 (°) strike, d1 (°) dip, and r1 (°) rake of the first nodal plane, s2 (°) strike, d2 (°) dip, and r2 (°) rake of the second nodal plane. The 26 October Mw 5.9 event is a double event (see the Finite-Fault Analysis section for explanation), thus we have two locations: 8669321(a) and 8669321(b).
Table S4. [Plain Text Space-delimited Table; ~7 KB]. Parameters of the kinematic model of the Mw 5.9 mainshock of the 26 October event.
Table S5. [Plain Text Space-delimited Table; ~11 KB]. Parameters of the kinematic model of the Mw 6.5 mainshock of the 30 October event.
Figure S1. Distribution of real-time seismic stations used in this work (triangles, green dotted triangles, and diamonds), and stand-alone, not-real-time stations deployed but not used in this work (gray squares and white dotted triangles); magenta stars are the 30 October Mw 6.5 (bigger) and 24 August Mw 6.0 events; and yellow stars are the 26 October Mw 5.9, 26 October 26 Mw 5.4, and 24 August Mw 5.4 events.
Figure S2. 1D gradient model modified after De Luca et al. (2009) for P- and S-wave velocities in the target area.
Figure S3. VP/VS reference value for the relocation of M <5.0 events; histogram representation of the VP/VS ratio distribution for the inverted dataset; values of VP/VS are calculated at each station for each of the selected events. (Inset) The mean value for the whole dataset, 1.8505 (σ: 0.1389; ε: 0.0005) has been calculated with the Wadati method by selecting only P and S readings with weight ≤ 0.
Figure S4. Distribution of key parameter values comparing the full dataset (solid black bars) and the final catalog (red bars): (a–c) longitude, latitude, and depth standard errors, (d) rms, (e) azimuthal gap, and (f) ML.
Figure S5. Time-dependent estimation of the completeness magnitude (Mc) retrieved by applying the GFT (red line) and the MBS (blue line) methods.
Figure S6. Three-component waveforms recorded at the CNE station. The acceleration is fitted up to 10 Hz. Blue and red lines indicate the P-onset arrival time of the two identified earthquakes.
Figure S7. Comparison between synthetic ground velocity (red lines) and recorded strong motions (blue lines) filtered between 0.02 and 0.5 Hz. Numbers in the left and right stripes represent the amplitude range in cm/s equal for all the inverted stations.
Figure S8. (a and b) Comparison between synthetic ground velocity (red lines) and recorded strong motions (blue lines) filtered between 0.02 and 0.5 Hz. Numbers in the left and right stripes represent the amplitude range in cm/s equal for all the inverted stations, with the exception of the three closest ones.
Bindi, D., L. Luzi, S. Parolai, D. Di Giacomo, and G. Monachesi (2011). Site effects observed in alluvial basins: The case of Norcia (central Italy), Bull. Earthq. Eng. 9, no. 6, 1941–1959, doi: 10.1007/s10518-011-9273-3.
Carannante, S., G. Monachesi, M. Cattaneo, A. Amato, and C. Chiarabba (2013). Deep structure and tectonics of the northern‐central Apennines as seen by regional‐scale tomography and 3‐D located earthquakes, J. Geophys. Res. 118, no. 10, 5391–5403.
De Luca, G., M. Cattaneo, G. Monachesi, and A. Amato (2009). Seismicity in central and northern Apennines integrating the Italian national and national network, Tectonophysics 476, 121–135, doi: 10.1016/j.tecto.2008.11.032.
Dreger, D. S., L. Gee, P. Lombard, M. H. Murray, and B. Romanowicz (2005). Rapid finite-source analysis and near-fault strong ground motions: Application to the 2003 Mw 6.5 San Simeon and 2004 Mw 6.0 Parkfield earthquakes, Seismol. Res. Lett. 76, no. 1, 40–48.
Dreger, D. S., and A. Kaverina (2000). Seismic remote sensing for the earthquake source process and near-source strong shaking: A case study of the October 16, 1999 Hector Mine earthquake, Geophys. Res. Lett. 27, no. 13, 1941–1944, doi: 10.1029/1999GL011245.
Hartzell, S. H., and T. H. Heaton (1983). Inversion of strong ground motion and teleseismic waveform data for the fault rupture history of the 1979 Imperial Valley, California, earthquake, Bull. Seismol. Soc. Am. 73, no. 6A, 1553–1583.
Herrmann, R. B., L. Malagnini, and I. Munafò (2011). Regional moment tensor of the 2009 L’Aquila earthquake sequence, Bull.Seismol. Soc. Am. 101, no. 3, 975–993.
Lomax, A., A. Michelini, and A. Curtis (2009). Earthquake location, direct, global-search methods, in Encyclopedia of Complexity and Systems Science, R. A. Meyers (Editor), Springer, New York, New York, 2449–2473.
Michelini, A., L. Margheriti, M. Cattaneo, G. Cecere, G. D’Anna, A. Delladio, M. Moretti, S. Pintore, A. Amato, A. Basili, A. Bono, et al. (2016). The Italian National Seismic Network and the earthquake and tsunami monitoring and surveillance systems, Adv. Geosci. 43, 31–38.
Moretti M., S. Pondrelli, L. Margheriti, L. Abruzzesse, M. Anselmi, P. Arrocoau, P. Baccheschi, B. Baptie, R. Bonadio, A. Bono, et al. (2016). SISMIKO: Emergency network deployment and data sharing for the 2016 central Italy seismic sequence (2017), Ann. Geophys. 59, no. 5, doi: 10.4401/ag-7212.
Saikia, C. (1994). Modified frequency–wave number algorithm for regional seismograms using Filon’s quadrature; modelling of Lg waves in eastern North America, Geophys. J. Int. 118, 142–158.
Wadati, K. (1931). Shallow and deep earthquakes, 3rd paper, Geophys. Mag. 4, 231–283.
Woessner, J., and S. Wiemer (2005). Assessing the quality of earthquake catalogues: Estimating the magnitude of completeness and its uncertainty, Bull. Seismol. Soc. Am. 95, no. 2, 684–698.
[ Back ]
