ANNALS OF GEOPHYSICS, 59, Fast Track 5, 2016; DOI: 10.4401/ag-7213 1 Investigating the effectiveness of rup- ture directivity during the August 24, 2016 Mw 6.0 central Italy earthquake. SPAGNUOLO ELENA*, CIRELLA ANTONELLA, AKINCI AYBIGE Istituto Nazionale di Geofisica e Vulcanologia, via di Vigma Murata 605, 00143 Roma (Italy) *elena.spagnuolo@ingv.it Abstract In this study we investigate directivity effects associated to the August 24, 2016 Mw 6, central Italy earth- quake taking into account the source rupture heterogeneities. We use the directivity predictor proposed by Spudich et al. (2004) which is derived from the isochrones theory. The directivity is computed using a source to site geometry and a focal mechanism. For its simplicity it can be computed once that a moment tensor solution is available. We use this technique to validate the real time solutions. Moreover, because the directivity predictor depends on the rupture velocity it can be used as a proxy to validate the possible rupture history. For the aforementioned reasons our method revealed fruitful for real time applications and helpful to constrain a few main rupture features for further analysis. I. INTRODUCTION he 2016 Amatrice earthquake (Mw=6.0) occurred in the Central Apennines (Italy) on August 24th at 01:36 UTC. The hypo- center has been located at 13.238E, 42.704N, at a depth of 8.1 km (http://cnt.rm.ingv.it). The earthquake caused casualties and heavy dam- ages in the Amatrice town and in several vil- lages nearby. Preliminary analyses have been conducted to retrieve fault and rupture pa- rameters in nearly real time. Cirella and Piata- nesi (2016) investigated the rupture process of the main shock, by inverting strong motion da- ta of 22 stations operated by RAN (Rete Accel- erometrica Nazionale,) and by INGV. Their preliminary source model is characterized by two principal patches (Figure1); a largest slip concentration, located at nearly 5 km NW from the nucleation, and a second smaller slip patch, located above the hypocenter. The retrieved front propagation shows a clearly bilateral rupture with heterogeneous rupture velocity. Strong motion data registered by the Italian strong motion network (RAN Rete Accelero- metrica Nazionale) which are equipped with Kinemetrics Episensor (FBA-3 200 Hz) and with ETNA 18 bits or K2-Makalu 24 bits digi- tizers, report amplitudes stronger than ex- pected with presumed azimuthal variability (ran.protezionecivile.it). Acceleration time se- ries recorded at nearby stations (epicentral dis- tances <50 km) show high variability in both amplitude and duration with PGA values ranging from 0.45g to 0.02g. Stations located in the N-W region of the source show relatively larger ground accelerations respect to those T ANNALS OF GEOPHYSICS, 59, Fast Track 5, 2016; DOI: 10.4401/ag-7213 2 located to the S-W. It is also observed that four ground motions recorded at Amatrice (AMT), Norcia (NRC), Montereale (RM33) and Fiastra (MNF) present impulsive characteristics over a multitude of orientations (INGV-RELUIS Re- port, 2016). Directivity could be implicated. Directivity has several derivations. Here we adopt the concept of directivity due to interfer- ing S waves. Waves radiated from points that rupture progressively on the fault plane inter- fere to modulate the energy in the direction of the site. The directivity effect can have im- portant implications for hazard. However, the validations of directivity model on the real da- ta can also trace-back important information regarding the rupture process on the fault. Here we apply a method based on the Iso- chrones Directivity Predictor by Spudich at al. (2004) that shows remarkable similarity with the observed azimuthal variations of the ground motion. Three models are proposed to evaluate possible contributions of mechanical heterogeneities of the fault to the observed ground motion amplifications. This study also suggests that a quasi-real time evaluation of the directivity effect can be used to put con- straints on a particular fault geometry or pro- vide hints on the rupture process. II. METHODS Empirical ground motion predictive equations (GMPEs) are commonly used to evaluate a ground motion intensity measure (i.e. the spec- tral acceleration Sa) from the source to site dis- tance and magnitude, at first-order. The ground motion predictive equation can be modified to incorporate directivity effects (Sa- DIR) using the following equation [Spudich and Chiou, 2008]: (1) ln (SaDIR) = ln (Sa) +fD where fD is the directivity effect defined as: (2) fD = Tapers(Rrup, M)*(a + b*IDP) a and b are empirical constant generally cali- brated respect to the available GMPEs, Ta- pers(Rrup, M) tapers the fD outside a region which provides the best correlation with the residuals (Rrup > 70 km and M < 5.6). The IDP is calculated using a close-form analytical solu- tion for the ground displacement for a finite fault in a whole space. It uses a parametric ray form derived from the isochrones theory [Spudich et al. 2004]. The IDP is calculated considering the source to site geometry, the hypocenter position, a fixed rupture velocity, rupture direction, and focal mechanism: (3) IDP = S Rri C where Rri is the radiation pattern that depends on the focal mechanism and assumes a water level of 0.2, S is a functional form which de- pends on source to site fault geometry, and C is a functional form of the isochrones velocity c’ defined as: (4) c’ = [1/Mach – (Rhyp – Rrup)/D]-1 Mach in Eq. (4) is the seismic Mach number de- fined as Mach=Vr/β where Vr is the rupture ve- locity and β is the S-waves speed. Rhyp is the hypocentral distance, Rrup is the distance to the fault rupture and D is the distance between the hypocenter and closest point in the source- to-site direction. For a complete description of the parameters in Eq. (4) we refer the reader to the original paper by Spudich and Chiou, [2008]. We also recall that a new model is re- cently available and developed for the Next Generation Attenuation (NGA) 2 project, ANNALS OF GEOPHYSICS, 59, Fast Track 5, 2016; DOI: 10.4401/ag-7213 3 which unfortunately was not implemented yet for application to the present case. In the present paper first we computed a fast and straightforward directivity evaluation by using the parameters listed in Table 1. We adopt the focal mechanism resulting from the nearly real time domain moment tensor solu- tion (TDMT) available at the INGV website (http://cnt.rm.ingv.it) within 6-9 minutes after the occurrence of an earthquake. Then we al- low for a variable rupture velocity and a depth dependent Vs profile to investigate the role of mechanical properties of both fault and propa- gating media. We avail of the preliminary kin- ematic inversion model resulting from Cirella and Piatanesi [2016] (Figure 1), which was resampled and processed as explained below. The advantage of this analytical approach is that it is fast and can be computed nearly real time as soon as few information about the seismic source become available. Table 1: Fault parameters and focal mechanism Parameter TDMT Mw 6.0 Width 17.5 Length 30 Dip 50 Rake -85 Rupture Velocity Variable Strike 156 Depth to top 0.5 In the present paper we perform a non- conventional calibration of the IDP on the real observations of the mainshock. The advantage of this technique is that it is independent from a particular choice of a GMPE. Equation (1) models how the observed ground motion de- viates from the mean value along a path of constant Rrup (the racetrack, Spudich at al. 2013). The mean value along the racetrack can be computed using the parametric equation in Spudich et al. [2013, appendix E]: (5) IDPM = a1 + - b1/ (cosh(c1 max(Rrup-z1,0)) Where a1, b1, c1, and z1 are depth dependent parameters provided in the author’s paper. The IDPM represents a general distance de- pendent model which can be used to remove the non-directive part of the motions from the data. We calculate the residuals t of the observations respect to the IDPM using a linear additive model with a constant term (@Matlab function regstats). The IDPM from Eq. (5) was then iter- atively shifted in amplitude to minimize the residuals. We thus calibrate the IDP on the ob- servations using a linear regression procedure: (6) t= a + b*IDP In the regression procedure we use stations with Rrup < 140 km in order to include a sig- nificant number of observations by doubling the distance tapering in Eq. (2). Parameters a and b are listed in Table 3. We also list in Table 3, parameters a and b resulting from the con- ventional calibration of Spudich and Chiou [2008] on the Abrahamson and Silva [2008] (hereafter AS08) equation for comparison. We calculated the fD (Eq. 2) for two cases, us- ing the fault geometry and focal mechanism listed in Table 1. In Model 1 we assumed con- ANNALS OF GEOPHYSICS, 59, Fast Track 5, 2016; DOI: 10.4401/ag-7213 4 stant rupture velocity and constant Vs value (Mach number = 0.8). In Model 2 we assume both a variable rupture velocity and depth de- pendent Vs profile (Table 2). Figure 1: Inverted rupture model [Cirella& Pi- atanesi, 2016] projected on the Earth surface with villages up to 20 km distant from the nucleation. Colors on the fault plane indicate the slip distribu- tion; black contours represent the position of the propagating rupture at 1 s interval. Violet line displays the location of the observed surface breakages (provided by S. Gori & E. Falcucci- EMERGEO WG). The use of a variable rupture velocity in the IDP model is not a standard practice. In this study we propose this implementation to veri- fy possible implications of mechanical hetero- geneities along the fault. However, the model resulting from this novel approach is still pre- liminary and possibly it introduces mathemat- ical singularities that should be considered with care. The rupture velocity model originates from the preliminary inversion results from Cirella and Piatanesi [2016, Figure 1]. We simplify the Mach number resulting from the depth dependent velocity profile and the rupture avoiding velocity inversions and su- pershear, considering only the high Mach numbers correlated to a substantial slip, set- ting maximum and minimum values in a typi- cal range between a water level of 0.7 and a maximum of 0.9. The resulting Mach number distribution along fault strike and down dip is shown in Figure 2. The grid spacing for slip distribution is 5 km. Ground motion results are plotted for a grid of sites uniformly distributed (2 km) around the fault. Table 2: Crustal Vs structure. Name Depth (km) Vs (km/s) 0.0 2.1436 1.5 2.1436 1.5 2.8210 4.5 2.8210 4. 3.4336 7.5 3.4336 7.5 3.1475 14.5 3.1475 14.5 3.3583 29.5 3.3583 29.5 4.0081 35.5 4.0080 35.5 3.9864 ANNALS OF GEOPHYSICS, 59, Fast Track 5, 2016; DOI: 10.4401/ag-7213 5 Figure 2: Along strike and down dip distribution of the Mach number=Vr/β. The hypocenter posi- tion is indicated with a star. The chosen ground motion intensity measure is the pseudo spectral acceleration PSA at T=3.0 s. This particular choice of T=3 s is due to: (1) the fact that the effectiveness of a di- rectivity model is larger at period greater than 1s (Spagnuolo et al. 2016); (2) fault geometry and details of the rupture process are more important at periods larger than 1s whereas site effects are less important; (3) to allow a di- rect comparison with the Shakemaps (resampled over the same grid of virtual sites), which are readily available after the revised localization of the event. Table 3: Parameters resulting from the residual regres- sion Eq.(6) Regression Parameters with 95% confidence bounds T=3 seconds a -0.063 (-0.241, 0.114) b 0.189 (0.07623, 0.562) a (AS08) -0.1979 b (AS08) 0.1319 III. RESULTS Figures 3 shows the directivity amplification for the two models proposed. The effect of di- rectivity (in color-bar scale) is expressed by the term [exp(fD) – 1] which represents a relative difference. If the directivity correction for a model ranges from -5% to 45% of the non- directive motion, then the color-bar axis runs from -0.05 to 0.45 [Spudich et al. 2013, PEER report]. Contour lines represent the PSA (3s) %g resulting from the Shakemap procedure [Faenza and Michelini, 2010, 2011]. The super- position is here meant to compare the spatial distribution of the real data (absolute values) with the spatial distribution of the directivity effect (relative values). The fD correlates fairly well with the expected azimuthal variations of the ground motion. The constant rupture velocity Model 1 (Figure 3, top panel) essentially underestimates the ground motion on the NE region and overes- timates the ground motion in the area along the surface fault projection. However, it cap- tures peculiar features in ground motion spa- tial distribution in the south-western fault edge and in the north-eastern region (i.e. BVG recorded PSA (3s) = 4.5 %g on the EW compo- nent). The variable rupture velocity fD Model 2 (Fig- ure 3, bottom panel) re-modulates the ground motion which results stronger than Model1 towards the SW edge of the fault with broader amplitudes lobes that amplify the ground mo- tion towards Amatrice (AMT). Importantly, the use of a variable rupture model re- modulates the up-dip directivity amplifica- tions along the 3 %g contouring of the Shake- ANNALS OF GEOPHYSICS, 59, Fast Track 5, 2016; DOI: 10.4401/ag-7213 6 map and better represents the ground shaking observed in the eastern area. Figure 3: fD resulting from Eq. (2). Model 1 (top): Constant rupture velocity and constant Vs result- ing in Mach = 0.8; Model 2 (bottom): variable rup- ture velocity and depth variable Vs profile (Table 2). Shakemaps results (from the peak horizontal value of 5% critically damped pseudo-acceleration, %g) are presented in white contouring. It is worth noting that the fD is not symmetric respect to the average IDPM. A stronger posi- tive contribution results from the source to site geometry (the hypocentral position in particu- lar) modulated by the radiation pattern (re- spectively S and C, Rri in Eq. 3). IV. DISCUSSION Our results show that the analytic model for directivity can be effectively used in the par- ticular case of the Amatrice earthquake to pre- dict the azimuthal variability of the ground motion once that a the fault geometry, the hy- pocentral position and the focal mechanism are known from preliminary reports [Tinti and Scognamiglio, 2016; Cirella and Piatanesi, 2016]. In this study we demonstrated that the variable Mach number (the ratio between the rupture velocity and S wave velocity), re - modulates the directivity effect along the top surface projection of the fault and introduces features that correlate with a few characteristic aspects of the azimuthal ground motion distri- bution. Moreover, the intervention of structur- al and geological control is not totally unex- pected since variable rupture velocity can orig- inate from the existing lithological contrasts between the synorogenic flysch of the NE por- tion of the fault (lower seismic velocities) and the carbonatic sequence of the SW portion of the fault [Gruppo di Lavoro INGV sul terre- moto di Amatrice, 2016]. Though still preliminary, the opportunity of introducing a variable rupture velocity is intri- guing as it allows complex features in a rather simple analytical formulation. The relative amplification values in Figure 3 are still not able to predict the strong ground motion recordings of the main event especially far from the rupture (Rrup > 70 km). This in- consistency is likely due: i) to the fD (Eq. 2) it- self which best performs for fault dip > 65° and in the near-field [Spudich and Chiou, 2008]; ii) to the comparison with the ShakeMap which interpolates the real data using a procedure that differs from the analytical formulation ANNALS OF GEOPHYSICS, 59, Fast Track 5, 2016; DOI: 10.4401/ag-7213 7 proposed here (Spagnuolo et al. 2013); iii) to the calibration of Eq. (2) that here results from a fast evaluation of the non-directive motion (Eq. 5) and could not exhaustively explain re- gional features. The latter Point iii) is related to the fact that first order effects superimpose their features on the directivity pattern (e.g. anelastic attenu- ation, site effects, topographic effects, rupture features such as the position of the patch of slip) and these effects are not taken into ac- count in the average non-directive model (Eq. 5) proposed here. In order to account for such effects, regionally derived GMPEs should be used to calibrate Eq. (2) for a and b. This strat- egy was originally proposed for the NGA da- tabase by Spudich and Chiu, (2008, see Table 3 for comparison) with the purpose of reducing the residuals between the observed (and not explained) azimuthal variability of the ground motion and the regional ground motion pre- dictive equation. The a and b parameters (Ta- ble 3) are markedly different in the case of AS08 (from Spudich and Chiou [2008]) and the average model proposed here. This difference highlights the relevance of selecting a proper ground motion prediction over which operat- ing the correction for directivity. In the case of the Italian database, the ground motion varia- bility shows a better agreement with Bindi et al. [2011], and Malagnini et al. [2011], than oth- er NGA models at almost all distances [Pischi- utta et al. 2016]. Unfortunately, GMEPs de- rived from Italian strong motion database do not provide corrective factors for the SC08 and do not include themselves a directivity term. The application to other Italian and European GMPEs need ad-hoc calibrations to enhance the performance of the IDP predictor. In this respect it is also worth stressing out that the directivity model developed by Spu- dich and Chiou [2008] provides the average IDP over a rich data of both synthetics and empirical recordings. Thus, it is likely that the use if the IDP for forward modeling of a single event is not able to predict the observations. However, the use of the IDP could be still ef- fective in a probabilistic framework when all the possible occurrences of rupture parameters are adequately accounted for [Spagnuolo et al. 2012]. Although the directivity correction is still not representative of the real contribution to the ground shaking, the azimuthal variations are well explained. Given the spatial agreement between the prediction and the observations for this particular earthquake, our results seem to indicate a strong contribution of the di- rectivity effect resulting from radiated waves interfering from the source to the site. The agreement also suggest the predictor can be effectively used in real time to seek the best geometry of the fault, relative position of the hypocenter, focal mechanism and rupture ve- locity, over a statistical consistent number of simulated directivity scenarios. REFERENCES Abrahamson N., and W. Silva, (2008). Summary of the Abrahamson & Silva NGA ground-motion relations, Earthquake Spectra, 24(1), 67–97. Bindi D., F. Pacor, L. Luzi , R. Puglia, M. Massa, G. Ameri, and R. Paolucci (2011). Ground motion prediction equations derived from the Italian strong motion database, Bulle- ANNALS OF GEOPHYSICS, 59, Fast Track 5, 2016; DOI: 10.4401/ag-7213 8 tin of Earthquake Engineering. doi 10.1007/s10518-011-9313-z Boore, D.M. and G.M., Atkinson, (2008). Ground-motion prediction equations for the average horizontal component of PGA, PGV, and 5 per cent damped PSA at spectral periods between 0.01 s and 10.0 s, Earthq. Spectra, 24, 99–138. Chiou, B.S-J., and R. R. Youngs, (2008). An NGA Model for the Average Horizontal Component of Peak Ground Motion and Re- sponse Spectra, Earthquake Spectra 24, no. 1, 173-215,doi:http:/dx.doi.org/10.1193 /1.2894832 Cirella A., and A. Piatanesi (2016). Source Complexity of the 2016 Amatrice Earthquake from Non Linear Inversion of Strong Motion Data: Preliminary Results; DOI:10.5281/zenodo.153821 Faenza, L., and Michelini, A., (2010). Regres- sion analysis of MCS intensity and ground mo- tion parameters in Italy and its application in ShakeMap. Geophys. J. Int.,180, 1138–1152. Faenza L. and A. Michelini, (2011). Regression analysis of MCS intensity and ground motion spectral accelerations (SAs) in Italy. Geophysical Journal International, vol. 186, pp. 1415–1439. Gruppo di Lavoro INGV sul terremoto di Amatrice (2016). Primo rapporto di sintesi sul Terremoto di Amatrice Ml 6.0 del 24 Agosto 2016 (Italia Centrale), doi: 10.5281/zenodo.61121 Malagnini, L., Akinci, A., Mayeda, K., Munafo, I., Herrmann, R.B.,` Mercuri, A., 2011. Charac- terization of earthquake-induced ground mo- tion from the L’Aquila seismic sequence of 2009, Italy. Geophys. J. Int., 184: 325–337. Pischiutta, M., Akinci, A., Malagnini L., and Herrero A., (2016). Characteristics of the Strong Ground Motion from the 24th August 2016 Amatrice Earthquake. Annals of Geophys- ics (This Issue). ReLUIS-INGV Workgroup, (2016). Preliminary study of Rieti earthquake ground motion rec- ords V5, available at http://www.reluis.it. Spagnuolo E., Herrero A., and G. Cultrera (2012). The effect of directivity in a PSHA framework, Geophyisical Journal Internation- al, 0.1111/j.1365-246X.2012.05630.x Spagnuolo E., L. Faenza, G. Cultrera, A. Herrero, and A. Michelini (2013). Accounting for source effects in the ShakeMap procedure: the 2000 Tottori and the 2008 Miyagi earth- quakesGeophys. J. Int. 194 (3), 1836-1848 2013doi:10.1093/gji/ggt195 Spagnuolo E., A. Akinci, A. Herrero, and S. Pucci (2016). Implementing the Effect of the Rupture Directivity on PSHA for the city of Istanbul, Turkey, Bulletin of Seismological So- ciety of America, Vol. 106, No.6, doi:10.1785/0120160020). Spudich P., B. S. J. Chiou, R. Graves, N. Col- lins, and P. Sommerville P. (2004). A Formula- tion of Directivity for Earthquake Sources Us- ing Isochrone Theory. USGS Open-File Report 2004-1268] http://dx.doi.org/10.1193/1.2894832 ANNALS OF GEOPHYSICS, 59, Fast Track 5, 2016; DOI: 10.4401/ag-7213 9 Spudich P. and B.S. Chiou, (2008). Directivity in NGA earthquake ground motions: analysis using isochrone theory, Earthquake Spectra, 24 (1), 279–298. Spudich P., et al. (2013) Final Report of the NGA-West2 Directivity Working Group, USGS Peer 2013/09. Tinti E. and L. Scognamiglio (2016). Finite fault rupture history of the 2016/08/24 Mw6.0 Ama- trice Earthquake by inverting strong motion data. Preliminary results. doi.org/10.5281/zenodo.61370