Kinematic Finite Fault and 3D Seismic Wave Propagation of the 24 August, 2016, M w 6.0 Central Italy Earthquake

The magnitude M w 6.0 earthquake of 24 th August 2016 caused severe damages and nearly 300 fatalities in the central Italy region. Initial reports revealed an asymmetrical distribution of damage and coseismic effects, suggesting a major role of heterogeneities, both in the rupture history and in the geological structure of the region. Near realtime availability of seismological data afforded a timely determination of a ﬁnite fault model (Tinti et al., 2016). Here we test this source model by performing a 3D simulation of seismic wave propagation within a 3D structural model containing the major geological features of the region. Agreement between modeled seismograms and observed seismograms suggests that some complexities in the waveforms, such as high ampliﬁcation in the region of the Mt. Vettore fault system, can be accounted for by complexities in the fault rupture and 3D structural models. Finally, the consistency of the hypothesis of two distinct events has been analyzed.


I. Introduction
T he 24 August, 2016, M w 6.0 central Italy earthquake occurred near the village of Accumoli with a strong impact in an area of complex geological structure. Large ground motions and damages were the result of complex combination of seismic source slip history, path effects, and site response. The near realtime availability of high quality geodetic and seismological data provided the possibility to constrain a preliminary rupture history of seismic source in a timely schedule ( [Tinti et al.(2016)]). Nevertheless, one of the primary assumptions in the procedure implemented by [Scognamiglio et al.(2010)] is the adoption of 1D model that performs remarkably well in the central Italian region ([Magnoni et al.(2014)], [Scognamiglio et al.(2016)]). In this paper, we probe this initial reliable finite fault solution with a forward simulation adopting a 3D velocity model ([Di Stefano & Ciaccio.(2014)]) that takes into account major geological features of the region. As in similar tests ([Magnoni et al.(2014)], [Dreger et al.(2015)]), for a 3D forward simulation we expect an incompatibility in adopting a seismic source solution obtained by inversion with 1D velocity models. Nevertheless a 3D approach is required in order to obtain a better understanding of the complexity that surprisingly appears in this relatively small event.

II. Methods
In order to perform simulations of the seismic waves generated by the 24 August, 2016, central Italy earthquake, we exploit a numerical Spectral-Element Method (SEM) that allows to model seismic wavefields in complex heterogeneous structure with very high spatial accuracy and computational efficiency (e.g., [Komatitsch et al.(2005)]).  [Tinti et al.(2016)], Slip amplitudes are expressed in meters and rupture times (white contouring) in seconds. The inferred rupture velocity is 3100 m/s. b) Rupture model used in this paper, described as seismic moment release, considering the material properties on the fault plane. c) and d) are the imaging of Vp and Vp/Vs on the fault plane, derived from Di [Di Stefano & Ciaccio.(2014)]. The transition zone between the slower NW part and the faster SE part corresponds to the Olevano-Antrodoco discontinuity.
In particular, we use the SEM software package SPECFEM3D_Cartesian designed for local/regional waveform simulations ( [Peter et al.(2011)]). The software allows us to accommodate all complexities that affect seismic wave propagation, such as topography, lateral wavespeed variations, attenuation, anisotropy, absorbing conditions at model boundaries as well as a finite source description. Specifically, we use a kinematic source model, although the code also allows for dynamic rupture modeling (e.g., [Huang et al.(2013)]). The input elements for the code are models of the wavespeed structure and of the seismic source. To describe the material properties of the simulated region, we adopt the 3D travel time tomography produced by [Di Stefano & Ciaccio.(2014)], that provides independent models of v p and v s for the whole Italian region. The mass density has been derived as a function of v p using the empirical relationship in [Magnoni et al.(2014)]. Attenuation is also included in the simulations with a simple model of Q as a linear increasing function of v s ( [Magnoni et al.(2014)]). Concerning the kinematic source description, we use the model released by [Tinti et al.(2016)] 48h after the event. This rapid finite fault inversion has been performed using the method of [Dreger and Kaverina. (2000)], strong motion data from RAN and INGV stations, and the 1D CIA velocity model ( [Herrmann et al.(2011)]). The assumed fault has length of 26 km, width of 16 km, strike of 156 • and dip of 50 • to the SW, supported by contextual geodetic analyses. The fault was parameterized using 2x2 km 104 sub-faults. The slip inferred from the inversion is heterogeneously distributed featuring two main shallow slip patches: the first one located 3 km up-dip from the hypocenter and the second 8 km toward North along strike direction. The rupture history shows bilateral propagation and pronounced directivity effects both N-NW and SE of the hypocenter. In order to implement this model in SPECFEM3D, we consider a fault with the same dimension,  (Fig. 1).
To perform the SEM simualations, the model volume is discretized using a mesh of hexahedral elements and the values of the ma-terial properties are assigned to each nodal point. The constructed grid covers a volume of 90x120x60 km around the epicentral area and honors the topography with a resolution on the surface of 500 m.

III. Results
SPECFEM3D outputs involve complete (i.e., including all possible phase arrivals), threecomponent seismograms in displacement, velocity and acceleration for the whole duration of the simulation (40 s) (red seismograms in Fig. 2a). Moreover, snapshots of the seimic wavefield at fixed interval of time steps are produced, from which we can reconstruct a 3D animation of the wavefront propagation in the Earth (Fig. 2b). Finally, synthetic shakemaps of the ground motion behavior can be reconstructed from the SEM simulation ( Fig. 3a-b). We compare the synthetic seismograms with the accelerometer data integrated to obtain ground velocities, after a lowpass-highpass filter processing in the frequency band of 0.02-0.5 Hz (Fig. 2a). The original finite fault solution is determined using CIA 1D velocity model ( [Herrmann et al.(2011)]), but even if it follows that there is an intrinsic, epistemological incompatibility between the source and the structural model when this is described in 3D, the 3D seismograms show a good agreement for the first arrivals and a richer contribution in the coda part, mainly related to the role of topography (red seismograms in Fig.  2a). Nevertheless, observing station AMT (the closest to the epicenter) and considering the almost perfect fit of the first phases presented in [Tinti et al.(2016)] (green seismograms in Fig. 2a), the 3D synthetic seismograms show the contribution of faster material (Fig. 2a) with a first pulse anticipated of 1.2 s. Such difference highlights the ordinary unresolved and ambiguous map between the patches in the finite fault solution and the presence of heterogeneities in the geological volume. Significantly, this positive shift in time is less pronounced in the stations on the North (see NRC, CPS) and West (CSC, RM33) of the fault plane.

IV. Discussion
Although additional improvements of the tomography and a new finite fault inversion with 3D velocity model are needed to better capture the pulses in the data at these frequencies, the amplitude of the signals is well constrained by the 3D simulations. Fig 3a and 3b report respectively the peak ground displacement (PGD) and peak ground acceleration (PGA) maps. The two main slip patches described in the finite source solution are mirrored on the surface with prominent pulses in the region of Amatrice and toward the village of Castelluccio, close to the slopes of Mt. Vettore. The pattern of ground motion peak is strongly influenced by the inclusion of topography, with higher peaks in correspondece of higher reliefs. The field surveying of coseismic evidences by EMERGEO significantly reports 5 km long concentration of surface ruptures in the Mt.
Vettore fault system region, where the combination of the directivity of the second seismic pulse with the topographic amplification releases higher accelerations along the mountain edges.
As already remarked, the finite fault solution proposed by [Tinti et al.(2016)] features two large distinct slip pulses of high peak slip values and a relatively high rupture velocity (3.1 km/s). The second main patch radiated 3 s after the nucleation. This pattern is in agreement with aftershocks distribution and also with preliminary inversion of InSAR data ([Gruppo di lavoro IREA-CNR & INGV]), even if with a less pronounced separation in the latter case. Fig. 4 shows our simulation of the wavefront propagation for the vertical component velocity wavefield on the fault plane surface proposed by [Tinti et al.(2016)], superimposed to the slip distribution. The seismic waves impact the area of the second slip patch 3 s after the nucleation, when the proposed finite fault solution sees the activation of the second pulse. This timing coincidence and the clear separation between the two main patches suggest that the kinematic model of [Tinti et al.(2016)] is compatible with a description of two slightly smaller events (M w 5.8 and M w 5.9) separated by 3 s, with the second dynamically triggered by the first. Far from being able to demonstrate this hypothesis in this paper, we highlight here that this M w 6.0 event shows unsuspected heterogeneous history, that plays an important role in the high level of damages in part of the affected region. To obtain a better understanding of this complexity and improve the constrains on seismic source, future determination of kinematic fault history should adopt 3D velocity models and topography, in order to include three-dimensional propagation effects due to lateral wavespeed variations, alluvial basins, river valleys and topographical edges.