Investigating the effectiveness of rupture directivity during the August 24 , 2016 Mw 6 . 0 central Italy earthquake .

In this study we investigate directivity effects associated to the Mw6.0 Amatrice earthquake 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 hypocenter 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 damages in the Amatrice town and in several villages nearby.Preliminary analyses have been conducted to retrieve fault and rupture parameters in nearly real time.Cirella and Piatanesi (2016) investigated the rupture process of the main shock, by inverting strong motion data of 22 stations operated by RAN (Rete Accelerometrica 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 Accelerometrica Nazionale) which are equipped with Kinemetrics Episensor (FBA-3 200 Hz) and with ETNA 18 bits or K2-Makalu 24 bits digitizers, report amplitudes stronger than expected with presumed azimuthal variability (ran.protezionecivile.it).Acceleration time series recorded at nearby stations (epicentral distances <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 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 Report, 2016).Directivity could be implicated.Directivity has several derivations.Here we adopt the concept of directivity due to interfering S waves.Waves radiated from points that rupture progressively on the fault plane interfere to modulate the energy in the direction of the site.The directivity effect can have important implications for hazard.However, the validations of directivity model on the real data can also trace-back important information regarding the rupture process on the fault.Here we apply a method based on the Isochrones 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 constraints on a particular fault geometry or provide 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 spectral acceleration Sa) from the source to site distance 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 calibrated respect to the available GMPEs, Tapers(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 solution 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: 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 depends on source to site fault geometry, and C is a functional form of the isochrones velocity c' defined as: Mach in Eq. ( 4) is the seismic Mach number defined as Mach=Vr/β where Vr is the rupture velocity 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 sourceto-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 recently available and developed for the Next Generation Attenuation (NGA) 2 project, 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 solution (TDMT) available at the INGV website (http://cnt.rm.ingv.it)within 6-9 minutes after the occurrence of an earthquake.Then we allow for a variable rupture velocity and a depth dependent Vs profile to investigate the role of mechanical properties of both fault and propagating media.We avail of the preliminary kinematic 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.Depth to top 0.5 In the present paper we perform a nonconventional 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 deviates 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]: Where a1, b1, c1, and z1 are depth dependent parameters provided in the author's paper.
The IDPM represents a general distance dependent 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 iteratively shifted in amplitude to minimize the residuals.We thus calibrate the IDP on the observations 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 significant 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 conventional 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, using the fault geometry and focal mechanism listed in Table 1.In Model 1 we assumed con-stant rupture velocity and constant Vs value (Mach number = 0.8).In Model 2 we assume both a variable rupture velocity and depth dependent Vs profile (Table 2).The use of a variable rupture velocity in the IDP model is not a standard practice.In this study we propose this implementation to verify possible implications of mechanical heterogeneities along the fault.However, the model resulting from this novel approach is still preliminary and possibly it introduces mathematical 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 supershear, considering only the high Mach numbers to a substantial slip, setting maximum and minimum values in a typical 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.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 directivity 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 direct 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 regression Eq.( 6) If the directivity correction for a model ranges from -5% to 45% of the nondirective 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 andMichelini, 2010, 2011].The superposition 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 overestimates the ground motion in the area along the surface fault projection.However, it captures peculiar features in ground motion spatial 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 component).
The variable rupture velocity fD Model 2 (Figure 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 motion towards Amatrice (AMT).Importantly, the use of a variable rupture model remodulates the up-dip directivity amplifications along the 3 %g contouring of the Shake-map and better represents the ground shaking observed in the eastern area.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 positive contribution results from the source to site geometry (the hypocentral position in particular) modulated by the radiation pattern (respectively S and C, Rri in Eq. 3).

IV. DISCUSSION
Our results show that the analytic model for directivity can be effectively used in the particular case of the Amatrice earthquake to predict the azimuthal variability of the ground motion once that a the fault geometry, the hypocentral 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), remodulates 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 distribution.Moreover, the intervention of structural and geological control is not totally unexpected since variable rupture velocity can originate from the existing lithological contrasts between the synorogenic flysch of the NE portion of the fault (lower seismic velocities) and the carbonatic sequence of the SW portion of the fault [Gruppo di Lavoro INGV sul terremoto di Amatrice, 2016].
Though still preliminary, the opportunity of introducing a variable rupture velocity is intriguing 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 inconsistency is likely due: i) to the fD (Eq.2) itself 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 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 regional features.
The latter Point iii) is related to the fact that first order effects superimpose their features on the directivity pattern (e.g.anelastic attenuation, site effects, topographic effects, rupture features such as the position of the patch of slip) and these effects are not taken into account 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 strategy was originally proposed for the NGA database 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 predictive equation.The a and b parameters (Table 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 operating the correction for directivity.In the case of the Italian database, the ground motion variability shows a better agreement with Bindi et al. [2011], and Malagnini et al. [2011], than other NGA models at almost all distances [Pischiutta et al. 2016].Unfortunately, GMEPs derived 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 Spudich 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 effective 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 directivity 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 velocity, over a statistical consistent number of simulated directivity scenarios.

Figure 1 :
Figure 1: Inverted rupture model [Cirella& Piatanesi, 2016] projected on the Earth surface with villages up to 20 km distant from the nucleation.Colors on the fault plane indicate the slip distribution; 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).

Figure 2 :
Figure 2: Along strike and down dip distribution of the Mach number=Vr/β.The hypocenter position is indicated with a star.
Figures 3 shows the directivity amplification for the two models proposed.The effect of directivity (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 nondirective 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 superposition 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.

Table 1 :
Fault parameters and focal mechanism