Rupture imaging for the 2016 Au­ gust 24, Mw=6.0 central Italy earthquake, from back­projection of strong­motion array data.

By extending the conventional Beam­Forming frequency­wavenumber power spectral estimate to the case of arbitrarily­shaped wavefronts, we obtained images of rupture propagation during the 2016 August 24, Mw=6.0 Central Italy earthquake. Using a set of strong­motion accelerometers, we evaluate the beam power along the travel time curves associated with synthetic sources spanning a model fault surface. This allows deriving time­dependent images of the distribution of energy radiation throughout the fault plane. Results indicate bi­lateral rupture propagation toward SE and NW, in rough agreement with surface co­ seismic displacement and surface damage pattern. To a first order, our results are also consistent with those obtained from full­waveform inversion of strong­motion data.


I. INTRODUCTION
n August 24, 2016, at 01:36:32 UTC a Mw=6.0 struck the central sector of the Apennines chain (Italy), causing almost 300 casualties and extensive destruc tion. According to the TimeDomain Moment Tensor (TDMT) solution dispatched by INGV soon after the event, the earthquake was caused by normal faulting with planes striking along the Apennines direction, i.e. SSENNW (GdL INGV, 2016;Fig. 1). Subsequent analyses of the GPS and Synthetic Aperture Radar (SAR) deformation patterns [GdL IREACNR & INGV, 2016] and aftershocks distribution al low constraining the causative fault to the O southwesterndipping nodal plane of the TDMT solution. This mechanism is consistent with the structural features of this sector of the Apennines, characterized by NNWSSEtrend ing, westdipping extensional Quaternary faults which are responsible for most of the de structive earthquakes that struck Italy over the last decades. In recent years, thanks to the growing availability of data from dense digital networks, backprojection of seismic array data allowed to track rupture propagation for sev eral major earthquakes worldwide (e.g., Kiser and Ishii, 2012;Meng et al., 2012, and refer ences therein). This approach is particularly at tractive, as it constrains the spatiotemporal evolution of the rupture solely on the base of the phase of coherent array signals. Thus, it does not require any detailed knowledge of Green's functions and fault geometry, or re strictive parameterizations of the rupture kine matics. Following these premises, in this note we present timedependent images of rupture propagation during the 2016, August 24 main shock by backprojection of Pwave recordings from two set of strongmotion accelerometers. We adopt a standard beamforming procedure (e.g., Abrahamson and Bolt, 1987), which seeks the maximum power of the recorded seismo grams along the travel time curves associated with trial source locations. The method suffers from low resolution, especially in the case of multiple, simultaneously emitting sources (e.g., Goldstein and Archuleta, 1987). Nonethe less, the retrieved pictures are consistent with both the coseismic displacement imaged by satellite interferometry and damage distribu tion at the surface.

II. DATA
For this study we used data from nine stations pertaining to the National Accelerometric Net work (hereinafter: RAN), owned and managed by the Italian Department for Civil Protection (see section on Data and Sharing Resources). These stations are all located west of the fault, at epicentral distances ranging between 26 and 61 km (Fig. 1). Preliminary inspection of the threecomponent recordings evidenced that the Swaves suffered important loss of coher ence even for small interstation distances. Moreover, those waves are characterized by a low frequency content (~0.3 Hz on velocity seismograms), which make them not appropri ate for investigating the dynamical evolution of the source over short time intervals. Due to these reasons, we conducted our analyses only on the P wavetrain as recorded at the vertical components of ground motion. Accelerograms were preprocessed to remove the offset and linear trend, highpass filtered with a two pole, singlepass 0.1 Hz Butterworth filter, and finally integrated in time to obtain ground ve locity. As for any other array processing scheme, beamforming requires that the wave field maintains significant coherence among the different array channels. Therefore, we conducted the analysis separately for two sub arrays, composed by stations [LNS, ANT, CTD, TRL, RTI] and [SPM, TRN, NRN, SNI], hereinafter referred to as subarrays AR1 and AR2, respectively (see Fig. 1). The vertical component seismograms for the mainshock as recorded by the two cluster of stations are il lustrated in Figure 2.

III. METHODS
Following hypocentral data and focal mecha nisms dispatched by INGV (http://cnt.rm.in gv.it/event/7073641; last accessed October 20, 2016), the mainshock's fault is modeled as a 25x16km plane striking 156° and dipping 50°, centered on the catalog hypocenter. This model fault is then discretized using a regular grid of nodes spaced by 0.5km along both the strike and dip directions. For each grid node, we calculate the theoretical travel times to all the array elements using a smoothed version of the reference 1D model reported by Caran nante et al. (2013) (Fig. 3). Station residuals are then calculated as the dif ference between the observed Pwave travel times, and those calculated for the grid node coinciding with the hypocenter. These residu als are then used to derive the corrected travel times D' which are stored for the subsequent utilization in beam forming (BF) estimates. The analysis proceeds by aligning the velocity seis mograms using the interstation time differ ences estimated via crosscorrelation to adjust the manuallypicked Pwave arrival times. The main advantage of this preprocessing is that it allows for the use of time windows which are shorter than the average propaga tion time across the array, thus providing more precise estimates of phase differences among array channels (e.g., Goldstein and Archuleta, 1991a). The array seismograms are then used for deriving timevarying estimates of source location on the fault plane using a conven tional, frequencydomain beam forming esti mator (e.g., Abrahamson and Bolt, 1987;Rost and Thomas, 2002). Separate tests with syn thetics (not shown here) demonstrated that this method, though less resolutive, provides results which are more reliable than those ob tained from much sophisticated approaches such as Capon's HighResolution (Capon, 1969) or MUSIC (Goldstein and Archuleta, 1987) estimates. For a given time frame, the evaluation of BF power spectra is conducted through the following steps: (a) Fourier transform of the array signal; (b) Given a reference frequency w 0 , calculation of the spatial covariance matrix R(w 0 ), by smoothing the crossspectral estimates over the three frequency bins centered at w 0 (eq. 21 in Abrahamson and Bolt, 1987); (c) For each fault grid node located at X, calcu lation of the array steering vector A whose ele ments are given by: where D' and T are the theoretical and ob served travel times, respectively, x is the posi tion vector of the generic array element, and N is the number of stations. Once applied to the crossspectral matrix, these terms have the ef fect of shifting the phase of the estimated spa tial crossspectra according to what predicted for a generic source located at X.
(d) Estimate of the BF power spectral estimator which, in matrix notation, is written as: where the H superscript indicates the complex conjugate transpose. When X coincides with the true source location, all the array cross spectra in R are brought in phase by the opera tor A, and their squared sum (i.e., the beam power P) will thus take a maximum. Steps (b d) are repeated for all the discrete frequencies spanning the 0.54 Hz frequency band, and a final broadband power spectrum is obtained by linearly stacking the narrowband BF spec tra (eq. 2) obtained at individual frequency bins (Gal et al., 2014). The entire procedure is iterated over succes sive time frames, thus obtaining timedepen dent images of the distribution of beam power over the modeled fault surface spanned by the synthetic sources X. The resolving capabilities of the method are tested using a set of five syn thetic sources distributed at the center and cor ners of the model fault plane. The source time function is given by a 2Hz Ricker wavelet, which is propagated across the two subarrays using the traveltime tables described earlier.
Seismograms at individual stations are then contaminated by white, Gaussian noise with a signaltonoiseratio equal to 50, and processed using the same time window and frequency band as for the real case (see next Section). The different sources are analyzed separately, so that the obtained results do not account for in terference phenomena associated with waves simultaneously issued by separate sources. Peaks of the synthetic BF spectra (Fig. 4) cor rectly recover the location of the different sources; in order to assess resolution, we fol low Goldstein and Archuleta (1991b) and con sider as significant only those P(X) values which are 1dB above the background spectral level. This latter quantity is defined as two standard deviations above the mean of the BF spectrum. Resolution is on the order of ~7km and ~5 km, respectively, along the downdip and alongstrike directions. In principle, the technique should thus be able to recover spa tiallydistinct point sources distributed throughout a fault surface whose inferred length and width are on the order of 2025km and 515km, respectively (GdL INGV, 2016).

IV. RESULTS
BF power spectra are calculated over subse quent 2slong time windows shifting along the seismograms with 1s increment. We lim ited the analysis to the first 5 seconds of signal, in order to avoid contamination by the Swave arrivals. For each time frame, the power spec tra obtained at the two subarrays are multi plied, so that the final backprojection image only contains the power contributions which are common to the two cluster of stations. Over the first 1.5s of rupturing, spectra are dominated by radiation from the hypocentral area, with a slight updip propagation toward NNW. Within the 1.75s2.25s time interval, bi lateral rupture propagation becomes evident, with rupture fronts propagating updip to ward SSE, and both up and downdip toward NNW. From the timing and location of those spectral peaks one gets rough estimates of rup ture velocities in the 2.53.5 km/s range, which is consistent with the 3.1km/s inferred by Tinti et al. (2016) from full waveform inversion. Time frames in between 2.75s and 3.75s are dominated by spectral peaks located NNW of the hypocenter, but at closer distances with re spect to those observed at previous times. This apparent paradox, that would imply a retreat of the rupture front, can be interpreted in terms of the interference between waves simul taneously radiated by the previouslyidenti fied, separate rupture fronts expanding toward opposite directions.

V. DISCUSSION AND CONCLUSION
In this work we used backprojection of strong motion records to obtain a first estimate of rupture behavior during the early few seconds of the 2016, Mw=6.0, Central Italy earthquake. The reduced number of available stations, and their large relative distances hindered applica tion of sophisticated, highresolution methods. White, dashed circles are isolines of equal distance from the hypocenter.
As a consequence, our images are likely blurred by the interference of waves simulta neously radiated by distinct portions of the ex panding rupture front. Nonetheless, the main features of the rupture history thus far pre sented are consistent with what reported by separate studies and observations. Our images indicate bilateral rupture, with peak energy ra diations located SE updip and NW of the hypocenter. This is consistent with (a) the ob served directivity focusing seismic energy mainly toward the NNW as well as to the SE, as indicated by marked pulses in strong mo tion records (Tinti et al., 2016); (b) the spatial distribution of surface damages [GdL INGV, 2016], and (c) the pattern of coseismic dis placement. The overall picture is also in agree ment with the results from waveform inver sion of strong motion data by Tinti et al. (2016), who inferred heterogeneous slip distribution characterized by two shallow slip patches lo cated updip and NW from the hypocenter, and bilateral rupture propagation with rupture velocity on the order of 3 km/s. The simplicity of the processing steps pre sented in this note suggests the possibility of deploying multiple seismic arrays at close dis tances from seismogenic faults, aimed at the realtime tracking of rupture evolution in case of large earthquakes. This would provide a quick estimate of fault size and hence earth quake magnitude, contributing to the robust ness of earlywarning systems and rapid dam age assessment.