ER3D: a structural and geophysical 3D model of central Emilia- Romagna (Northern Italy) for numerical simulation of earthquake ground motion

ER3D: a structural and geophysical 3D model of central EmiliaRomagna (Northern Italy) for numerical simulation of earthquake ground motion Peter Klin, Giovanna Laurenzano, M. Adelaide Romano, Enrico Priolo, Luca Martelli 5 Centro Ricerche Sismologiche (CRS), Istituto Nazionale di Oceanografia e di Geofisica Sperimentale (OGS), Sgonico (TS), 34010, Italy Servizio Geologico Sismico e dei Suoli, Regione Emilia-Romagna, Bologna, 40127, Italy Correspondence to: Peter Klin (pklin@inogs.it) 10 Abstract. During the 2012 seismic sequence of Emilia region (Northern Italy), the earthquake ground motion in the epicentral area featured longer duration and higher velocity than those estimated by empirical-based prediction equations typically adopted in Italy. In order to explain these anomalies, we (1) build up a structural and geophysical 3D digital model of the crustal sector involved in the sequence, (2) reproduce the earthquake ground motion at some seismological stations through physics-based numerical simulations and (3) compare the observed recordings with the simulated ones. In this way 15 we investigate how the earthquake ground motion in the epicentral area is influenced by local stratigraphy and geological structure buried under the Po Plain alluvium. Our study area covers approximately 5000 km and extends from the Po river right bank to the Northern Apennines morphological margin in N-S direction, and between the two chief towns of Reggio Emilia and Ferrara in W-E direction, involving a crustal volume with 20 km of thickness. We set up the 3D model by using already published geological and geophysical data, with a detail corresponding to a map at scale 1:250,000. The model 20 depicts the stratigraphic and tectonic relationships of the main geological formations, the known faults and the spatial pattern of the seismic properties. Being a digital vector structure, the 3D model can be easily modified or refined locally for future improvements or applications. We exploited high performance computing to perform numerical simulations of the seismic wave propagation in the frequency range up to 2 Hz. In order to get rid of the finite source effects and validate the model response, we choose to reproduce the ground motion related to two moderate-size aftershocks of the 2012 Emilia sequence 25 that were recorded by a large number of stations. The obtained solutions compare very well to the recordings available at about 30 stations, in terms of peak ground velocity and signal duration. Snapshots of the simulated wavefield allow us to explain the exceptional length of the observed ground motion as due to surface waves overtones that are excited in the alluvial basin by the buried ridge of the Mirandola anticline. Physics-based simulations using realistic 3D geo-models show eventually to be effective for assessing the local seismic response and the seismic hazard in geologically complex areas. 30

Correspondence to: Peter Klin (pklin@inogs.it) 10 Abstract. During the 2012 seismic sequence of Emilia region (Northern Italy), the earthquake ground motion in the epicentral area featured longer duration and higher velocity than those estimated by empirical-based prediction equations typically adopted in Italy. In order to explain these anomalies, we (1) build up a structural and geophysical 3D digital model of the crustal sector involved in the sequence, (2) reproduce the earthquake ground motion at some seismological stations through physics-based numerical simulations and (3) compare the observed recordings with the simulated ones. In this way 15 we investigate how the earthquake ground motion in the epicentral area is influenced by local stratigraphy and geological structure buried under the Po Plain alluvium. Our study area covers approximately 5000 km 2 and extends from the Po river right bank to the Northern Apennines morphological margin in N-S direction, and between the two chief towns of Reggio Emilia and Ferrara in W-E direction, involving a crustal volume with 20 km of thickness. We set up the 3D model by using already published geological and geophysical data, with a detail corresponding to a map at scale 1:250,000. The model 20 depicts the stratigraphic and tectonic relationships of the main geological formations, the known faults and the spatial pattern of the seismic properties. Being a digital vector structure, the 3D model can be easily modified or refined locally for future improvements or applications. We exploited high performance computing to perform numerical simulations of the seismic wave propagation in the frequency range up to 2 Hz. In order to get rid of the finite source effects and validate the model response, we choose to reproduce the ground motion related to two moderate-size aftershocks of the 2012 Emilia sequence 25 that were recorded by a large number of stations. The obtained solutions compare very well to the recordings available at about 30 stations, in terms of peak ground velocity and signal duration. Snapshots of the simulated wavefield allow us to explain the exceptional length of the observed ground motion as due to surface waves overtones that are excited in the alluvial basin by the buried ridge of the Mirandola anticline. Physics-based simulations using realistic 3D geo-models show eventually to be effective for assessing the local seismic response and the seismic hazard in geologically complex areas.

1 Introduction
Computer aided three-dimensional (3D) geological modeling (e.g. Mallet, 2002) is becoming an increasingly important tool in geo-science studies for both the management of natural resources and the prevention of natural disasters. 3D geological modeling allows the combination of multi-disciplinary data in the shaping and visualization of the current knowledge of the geological structures and allows the integration with new data or interpretations, as they become available (Calcagno, 2015). 35 Moreover, 3D geological models represent the basis for the execution of physics-based numerical simulations, provided that a reliable scientific procedure is defined to convert the different types and levels of the available complex geological information into that needed by the proposed numerical simulation at the predefined scale level (e.g. Fischer et al., 2015).

2
The present study concerns the set-up of a 3D structural model starting from geological data and the development of the corresponding geophysical model by assigning visco-elastic properties to each structural unit. The scope of the final 3D geophysical model is to allow physics-based forward modeling of seismic wave propagation aimed at 1) explaining the ground motion peculiarities observed in past earthquakes and 2) increasing the reliability of ground motion predictions for possible future events. (e.g. Moczo et al., 2014;Taborda and Roten, 2015;Cruz-Atienza et al., 2016). Our study focuses on 5 the Emilia region (Northern Italy), where in 2012 a relevant seismic sequence featuring the two main shocks MW 6.1 on 20/5/2012 at 02:03:53 UTC and MW 5.9 on 29/5/2012 at 07:00:03 UTC (Rovida et al. 2016), occurred (Fig. 1).
Seismic-hazard studies are usually based on the empirical-statistical method, which makes use of ground motion prediction equations (GMPE) (e.g.: Barani et al., 2017a and2017b), with possible corrections deduced from local geological conditions (Grelle et al., 2016). However, occasionally the observed ground motion characteristics deviate considerably from the 10 empirical-statistical predictions. Those deviations imply the presence of case-specific features in wave generation or propagation (e.g., complex fault ruptures, complex geological structures, such as deep basins), which are not adequately considered in the derivation of the GMPE. In order to predict the effects of these features we may apply numericaldeterministic methods.
An emblematic case of such deviations occurred during the 2012 Emilia seismic crisis, when unexpectedly long duration and 15 large peak ground velocity (PGV) characterized the earthquake ground motion at some sites in the epicentral area (Castro et al., 2013;Luzi et al., 2013;Barnaba et al., 2014;De Nardis et al., 2014). Those deviations have been attributed to the complexity of the geological structure beneath the Po Plain, which features a very large and deep alluvial basin bounded by two largely buried thrust-and-fold systems, the Northern Apennine chain at South and the Southern Alpine ridge at North, respectively (Boccaletti et al., 1985). In order to explain quantitatively the observed ground motion characteristics, we have 20 built a 3D model that describes the morphology of the buried geological structure and assigns visco-elastic properties (mass density, elastic modula and elastic quality factors), to each formation, so that it can be used for physics-based numerical modeling of the seismic wave propagation in the studied volume.
Our model is not the first 3D model that was developed for the Po Plain area. At least three research groups have carried out 3D numerical computations of the earthquake ground motion in the Po Plain so far. A first study was performed by Vuan et 25 al. (2011), who simulated long-period (T > 5 s) surface waves generated in the basin by strong (MW > 6) earthquakes. They used a 3D model featuring realistic, irregular basin edges and a simplified depth-dependent velocity profile for the sedimentary filling of the basin. A more complex 3D geological model was set-up by Molinari et al. (2015) for simulating the earthquake ground motion in the long-period band (T > 3 s). The simulated waveforms were compared only qualitatively with the recorded waveforms at some far-source stations in order to demonstrate the effectiveness of the 3D geological 30 model. A third model is that one developed by Paolucci et al. (2015), who simulated the near-source strong ground motion for the MW 6.1 20 May 2012 earthquake in the frequency range 0.1-1.5 Hz. The overall satisfactory agreement of their simulated waveforms with the empirical records is due to two key-elements: the extended source model (i. e. slip distribution and rupture propagation) and the 3D structural model, which contains only two main geologic interfaces, (i. e. the base of the Pliocene formation and that of the Quaternary deposits). In particular, the satisfactory simulation of the surface waves trains 35 stems mainly from the shape of the interface of the base of the Quaternary deposits. We have to mention also Turrini et al. (2014), who defined the whole structure of the Po basin from its deep roots, at the Moho level, through an exhaustive analysis of all the existing structural-geological and geophysical studies. They summarize the current knowledge of the Po basin structural-geology into a digital, editable model that can be used to improve the geodynamic interpretation of the area. However, their model does not contain any geophysical parameterization, neither it reaches the level of detail that is required 40 in our study.
3 Among the cited works only Paolucci et al. (2015) provided the elements for understanding the peculiar features of the nearsource strong-motion observed during the 2012 events (such as the propagation of prominent trains of surface waves in the Northern direction), by adopting a reasonably simple 3D model of an area centered on the 2012 MW 6.1 mainshock epicenter. In the present work we instead focus on the southern sector of the 2012 epicentral area, characterized by a very deep basin with sediment thickness exceeding 8000 m in some points. In order to investigate the effects of this complex 5 geological setting, we set up a 3D geological model with unprecedented detail of a limited area of the Po Plain, bounded by the Po river right bank at North, by the Northern Apennines morphological margin at South, and located between the two chief towns of Reggio Emilia at West, and Ferrara at East (Fig. 1). The area includes the epicenters of the 2012 seismic sequence as well as some other potential seismogenic structures (DISS Working Group, 2015). In order to set up the 3D geological model we considered only data available in scientific literature. The physical properties assigned to the geological 10 units were deduced from literature as well. We relied on the commercial software GeoModeller released by Intrepid Geophysics for merging and interpolating the geological data in a 3D digital model, which constitutes the input for numerical simulations of the earthquake ground motion and represents a basis for further improvements when new data will be available. In this first version, denoted ER3D, the model is based on the elaboration of a Digital Terrain Model, a seismotectonic map, three deep geological sections crossing the study area as well as the isobaths of two interfaces between 15 some relevant geological formations. As discussed in section 3.2, the detail level included in this model is consistent with numerical computations of the ground motion in the frequency range up to 2 Hz, therefore comparable with the frequency range of the computations performed by Paolucci et al. (2015).
We performed the computations of the ground motion by applying the high-performance computing (HPC) code FPSM3D (Klin et al., 2010), dedicated to the numerical modeling of the propagation of viscoelastic waves in heterogeneous media.

20
The code FPSM3D is based on the Fourier pseudo-spectral method for the solution of hyperbolic equations; its accuracy performance compares well with other computer codes used in the scientific community for the 3D simulation of the earthquake ground motion in alluvial basins, as emerged during recent verification exercises (Maufroy et al., 2015;Chaljub et al., 2015). The validation of the constructed 3D geological model consisted in a comparison in terms of PGV and duration, between the numerical predictions and the empirical recordings of two 2012 events at the several stations that were deployed 25 in the area during the seismic sequence ( Fig. 1). In order to put in evidence the peculiarities of the ground motion that are due only to propagation effects, we considered two weak events (MW 4.0 and 4.1), that can be modeled using a point source, and not the mainshocks, which would require a finite source model. The computations were run using the HPC resources of the CINECA consortium in Bologna.

30
The fundamental step for physics-based numerical prediction of the earthquake ground motion consists in the set-up of a 3D model of the geological structure. In order to set up a reliable geological model we need a sound geological interpretation of well-constrained geophysical data. Thanks to oil exploration and research widely undertaken since 1960, a comprehensive synthesis of the structural features of the Po Plain subsurface was possible in the past decades (e.g. Pieri and Groppi, 1981;Fantoni and Franciosi, 2010;Boccaletti et al., 2011;Martelli et al., 2017). In the following we give an overview of the 35 known geological features of the study area and describe how we synthetized these data in a digital 3D structural model.
Finally, we discuss how we assigned the physical properties to each geological formation for characterizing the 3D model also from a geophysical point of view. 4

Geological and seismotectonic setting of the study area
The study area is in the Emilia-Romagna region (Northern Italy), and specifically it occupies the sector of the Po Plain between Reggio Emilia (West) and Ferrara (East), as shown in Fig. 1. The Po Plain is a foredeep-foreland zone interposed between two chains with opposite vergence: Northern Apennines to the South and Southern Alps to the North. Terrigenous sediments originating from the erosion of the two growing chains accumulated in the basin (Dondi et al., 1982), first those of 5 alpine origin (Miocene-Quaternary) then those of Apennine provenance (Pliocene-Quaternary). From the Middle Pleistocene the sedimentation is mainly continental and results from the depositional activity of the Po River and its tributaries. The substrate of the terrigenous sediments is made up by a carbonatic succession of mainly Mesozoic age, whose top consists of marly sediments of Paleogene age. This carbonate succession is separated from the metamorphic basement by a thick evaporite succession of Triassic age (Fig. 2). From the tectonic point of view, the area is affected by numerous compressive 10 structures, with northern vergence (Fig. 1). The southern zone, coinciding with the Apennines hills between Albinea, Sassuolo, Vignola and Casalecchio di Reno municipalities, is characterized by the reverse faults of the Pedeapenninic Thrust (Boccaletti et al., 1985), which is responsible for the morphological transition between the Northern Apennines and the Po Plain. Subsoil investigations for oil exploration (Pieri and Groppi, 1981) showed that the Apennine outer front does not coincide with the Apennine -Po Plain morphological margin and that in the Po Plain subsoil many blind faults and folds are 15 present. Actually, the Apennine outer front is located in the subsoil around the present course of the Po River, coinciding with the reverse faults of the Ferrara Folds overthrusting the Lombardy-Veneto monocline ( The relationships between tectonic structures, sedimentary bodies and the surface morphology indicate that Pedeapenninic Thrust and Ferrara Folds were active also in recent times, as demonstrated by the Quaternary deposits which are deformed and uplifted. Conversely, the Emilia and Romagna Folds were active mainly in the Pliocene, being the Quaternary deposits 30 not deformed by these structures but included in the syncline between the Pedeapenninic Thrust and the Ferrara Folds (Pieri and Groppi, 1981;Burrato et al., 2003;Boccaletti et al., 2004 andMartelli et al., 2017).
The entire area is seismically active, and the distribution of historical and instrumental earthquakes seems to confirm the major actual activity of the Pedeapenninic Thrust and the Ferrara Folds. In fact, the main historical earthquakes of the area have been located along the Apennines-Po Plain margin (Rovida et al., 2016), while Ferrara Folds are responsible for 35 earthquakes of 20 and 29 May 2012, respectively MW 6.1 and MW 5.9 (Fig. 1). For these reasons, these structures are included in the database of the seismogenic structures capable of generating strong earthquakes (DISS Working Group, 2015). The instrumental data indicate that in this area the most part of earthquakes has a compressional source mechanism (Pondrelli et al., 2006), and that the hypocentral depth (http://cnt.rm.ingv.it, last accessed May 2018) in the northern zone

Integration of geological data in the 3D digital model
Geological 3D modeling consists in the representation of a solid Earth sector by using surface and subsurface data in a computer-aided process (Mallet, 2002), which allows to shape and to visualize the current knowledge and/or to update it with new data. Numerous methodologies were implemented in several packages dedicated to the geological 3D modeling.
The package we adopt for the present work is GeoModeller (Calcagno et al., 2008;Guillen et al., 2004), a commercial 5 software originally developed by the French Bureau de Recherches Géologiques et Minières (BRGM) and more recently by the Intrepid Geophysics (http://www.geomodeller.com, last accessed May 2018). GeoModeller is a software tool for the integration of different geometrical, geological, and geophysical data in a geometrically coherent 3D geological model. The procedure is based on the potential field interpolation (Lajaunie et al., 1997) and is particularly well suited when the available geological data consist only in some geological maps, sparse cross-sections or boreholes. The method requires in 10 input the location of the geology interfaces and orientation data at some points. The theory behind the method describes the 3D geologic surfaces as iso-potential surfaces of a scalar potential-field with orientation vectors playing the role of the field's gradient. The stratigraphic pile is defined by the chronological order of the strata and the relationships between the formations in terms of either 'onlap' or 'erode'. The complex geology is described by different domains, each characterized by a geological serie, separated by stratigraphic or tectonic discontinuities. For each domain, the geology is modeled by a set of 15 sub-parallel, smoothly curving surfaces using the potential-field functions. Cokriging is used to obtain a solution that honors the input data (McInerney et al., 2005). Faults are taken into account as discontinuous drift functions into the cokriging equations (Chilès et al., 2004). Refer to Calcagno et al. (2008) for a more comprehensive description of the methods implemented in GeoModeller..
To build the 3D geological model of the central Emilia we considered a crustal volume with 70 km x 70 km area and depth 20 of 20 km, in order to include most of the 2012 seismic sequence hypocenters, associated to the deepest segments of the active thrusts. We defined the stratigraphic pile according to that one reported on the seismotectonic map of the Emilia-Romagna region (Boccaletti et al., 2004 andMartelli et al., 2017). We imported in Geomodeller the following data: The geological cross-sections of Boccaletti et al. (2004) and Boccaletti et al. (2011) are based on more recent seismic profiles than those used by Pieri and Groppi (1981) and take into account also stratigraphic data derived from RER and ENI-Agip (1998), for the definition of the superficial part (down to a depth of approximately 300-400 m).
To constrain the geometry of the geological bodies, we manually digitized the interfaces separating the oldest geological formations, both on the excerpt of the seismotectonic map and on three mentioned geological cross-sections, and merged 40 6 them with the digital data outlining the isobaths of the two youngest formations bottoms. Similarly, we digitized the fault traces on the cross sections and attributed them their extension, their relationship with the geology series (in order to take in consideration the faults when interpolating the geology series) and their relationship with other faults (to define the termination of one fault on another). The building process consists in several steps, with a progressive integration of the available data, generally starting from the top surface. At each step we perform a computation of the implicit surfaces of the 5 formation boundaries and review the partial result before adding new elements. Whether new data were available for the project, we can obtain a revised model with little effort.
The 3D model obtained for the Emilia region is displayed in Figs. 4, 5 and 6. Two different views of the model sampled on the three input geological cross-sections are shown in Fig. 4. Figure 5 shows four parallel equally spaced North-South 2D vertical sections across the investigated volume,, while Fig. 6 evidences the surfaces corresponding to the fault system.

10
In order to use the model for numerical simulations, we exported it into a 'voxet' format by sampling the geological formation volumes with a regular 3D grid.

Physical properties of geological formations
In order to perform the physics-based numerical simulations of the seismic waves propagation we have to assign the values of the physical properties to each 3D geological volume. Considering valid the assumption of an isotropic and visco-elastic 15 medium, we assigned to each formation the values of the following parameters: the velocities VP and VS of the compressional and the shear seismic wave, respectively the mass density ρ the elastic quality factors QK and Qµ for bulk and shear deformations, respectively.
We assumed that each geological formation belonging to the stratigraphic pile of Fig. 4 is characterized by different values 20 of the above-mentioned parameters. In order to simplify the assignment of the physical properties values to each formation, we decided to characterize each unit only by VP and to evaluate from the latter the other four properties using some well-established empirical relations.
Considering the velocities expressed in km/s we adopted the relation: VS =0.7858 − 1.2344·VP +0.7949·VP 2 − 0.1238·VP 3 +0.0054·VP 4 (1) 25 that was found by Brocher (2005) from a large number of measurements made in a variety of lithologies including Quaternary alluvium and Miocene sedimentary rocks, which constitute a fundamental part of our model. We also adopted the well established relation: found by Gardner et al. (1974) for the mass density  in g/cm 3 and VP in km/s.

30
The intrinsic attenuation is described with the shear quality factor, which is evaluated from VS with the widely used rule of thumb (e.g. Paolucci et al., 2015) Qµ =100·VS (km/s) and with the bulk quality factor, which value is set as

35
in accordance with the theory exposed by Morozov (2015).

7
We assumed that the value of VP assigned to each geological formation might depend on the depth through a linear gradient The values of the coefficients VP(0) and ∂zVP for each formation are given in  Montone and Mariucci (2015).
We tested the validity of eq. (1) by analyzing the consistency of the predicted Vs with some measures of Vs resulting from geophysical surveys performed in the Po plain. According to eq. (1), the value Vp=1.5 km/s assigned to the uppermost formation A (table 1)

The FPSM3D code
In order to compute the seismic waves propagation in the constructed 3D geological model we adopted the code FPSM3D (Klin et al., 2010), which is based on the Fourier pseudo-spectral method (FPSM) for the integration of hyperbolic equations.
The peculiarity of FPSM (first introduced by Kreiss and Oliger, 1972) consists in the evaluation of the spatial derivatives by means of a multiplication in the wavenumber domain. The transition from the spatial domain to the wavenumber domain, 30 and back, is performed by means of the Fast Fourier Transform. FPSM combines the optimal accuracy of the global spectral differential operators and the simplicity of the spatial discretization with a structured rectangular grid. According to the Nyquist's sampling theorem, FPSM works with a relatively coarse spatial sampling (Fornberg, 1987), which represents a valuable advantage when solving 3D problems. The code FPSM3D performs the time integration by means of the 2 nd -order explicit Finite-Difference (FD) scheme and adopts the Convolutional Perfectly Matching Layer (C-PML) approach 35 (Komatitsch and Martin, 2007) to prevent the effects of the spatial domain boundaries on the computed wavefield. The effects of the staircase approximation of the material interfaces in the regular grid are avoided using the volume harmonic averaging of the elastic moduli and volume arithmetic averaging of the mass density, as proposed by Moczo et al. (2002).
The adequateness of the FPSM3D code in this kind of applications is demonstrated in the works by Chaljub et al. (2015) and 8 Maufroy et al. (2015), aimed to estimate the accuracy of a number of numerical methods currently used for physics-based predictions of earthquake ground motion in 3D models of sedimentary basins.

Setup for the computations
A critical step in the setup for the numerical simulations consists in the choice of the frequency range. In order to reproduce accurately the wave propagation at high frequencies it is required a fine spatial and temporal sampling and therefore a larger The numerical computations were performed using the spatial and temporal sampling as exposed in Table (

Comparisons between numerical predictions and data
In order to investigate whether the ER3D model is capable to reproduce the peculiar features of the observed earthquake ground motion, we performed a comparison between the ground motion recorded by 29 seismological stations deployed in 25 the study area during the 2012 seismic sequence (see Fig. 1 and Table 4) and the numerically predicted ground motion at the same locations. The considered seismic stations belong to the Italian Strong Motion Network (IT) managed by DPC, and to the Italian National Seismic Network (IV) managed by INGV. For reference, we considered also the physics-based numerical predictions resulting from the simplified model PADANIA (Malagnini et al., 2012), which is composed of horizontal homogeneous layers (therefore a 1D model in contrast to our 3D). The numerical simulations regarding the 30 PADANIA model were performed using the Wavenumber Integration Method (WIM) (Herrmann and Wang, 1985;Herrmann, 1996a, b), which solves the wave equation in a horizontally layered medium. The synthetic seismograms contain all the phases and are accurate in both the near-and far-field.

The simulated earthquake ground motion
In order to focus the study on the propagation of seismic waves rather than on their generation, we decided not to consider 35 the main shocks of the 2012 Emilia sequence, which were simulated in previous works concerning the 3D modeling of the Po Plain (e.g. Molinari et al., 2015;Paolucci et al., 2015). For those events, of magnitude MW 6.1 (on 20/5/2012) and MW 5.9 (on 29/5/2012) respectively, the effects of the peculiarities of the seismic source in the recorded waveforms are not 9 negligible. Since the main topic of this work is the estimation of the wave propagation effects on the earthquake ground motion (in particular the influence the Po Plain underground geological structure has on the wave propagation), we decided to simulate events of lower magnitude. With an upper frequency limit of 2 Hz (see section 3.2), we can roughly assume that the complexities (i.e. unpredictable irregularities in the spatial extension and time evolution) in the seismic sources are negligible for earthquakes up to MW = 4.0. Nevertheless, such events are strong enough to be well recorded in almost all the 5 considered stations. We therefore computed the seismic wavefield for the two MW ≈ 4.0 events listed in Table 3 with sources located at the NE and NW ends of the studied area (c and d labeled events in Fig. 1). The hypocenters and the magnitude of these events were taken from the latest relocation study . The generation of the wavefield was modelled as a double-couple point source with a time function corresponding to the low-pass minimum phase Butterworth filter plotted in Fig. 7 and with an inverse focal mechanism, in accordance to the fault plane solutions of the 2012 sequence 10 found by Saraò and Peruzza (2012).

Comparison with the empirical earthquake ground motion
The permanent and temporary seismological stations deployed in the study area during 2012 are mapped in Fig. 1 and listed in Table 4. The time series recorded at these stations during the two events listed in (Table 3) Fig. 1. Since the empirical time-series have a much wider frequency content than the simulated ones, they were low-pass filtered using the same minimum phase Butterworth filter plotted in Fig. 7 that was used as the source time function in the numerical simulations.
We compared the simulated ground motion with the empirical one in terms of horizontal peak ground velocity (PGV) defined as the peak modulus of the vector sum of the two horizontal components and in the duration defined as the time 20 interval length between 5% and 95% of the Arias intensity (Arias, 1970). The vertical component was excluded from this comparison since it was systematically lower than the horizontal ones.
In Fig. 8 we plotseparately for the two eventsthe logarithm of the measured and computed PGV at each station against the epicentral distance. We represent there also the linear fit for the three series of data (empirical, 3D model ER3D, 1D model PADANIA). The plot shows that, in both cases, the 3D model numerical predictions fit better the observations, 25 whereas the 1D model prediction underestimate the observed PGV at most stations (by a factor of almost 2). The high variability shown by stations at similar epicentral distance is probably due to the different source-station azimuth and focal mechanism-radiation pattern. As observed in Maufroy et al (2015), the uncertainty in source characteristics may impact the numerical predictions especially at short distances. The remarkable underestimation of PGV for the event 2 at station T800, located just above the hypocenter is therefore not too surprising and could be attributed to the combined effect of inaccurate 30 hypocentral location, focal mechanism, and near-source heterogeneities. In fact, considering that source 2 has a dip of 33° (Table 2), T800 is near to the P-wave radiation maximum and at the margin of the S-wave lobe. Similarly, in Fig. 9 we plot the duration of the measured and computed ground motions against the epicentral distance.
Again, the 3D model numerical predictions fit better the observations than the 1D model predictions, which underestimate the duration at almost all the stations. In particular, it can be observed how the 3D model is able to predict quite well the very long duration values observed at some stations located in the southern part of the model (for example the MDN, MODE and 40 10 ZPP stations for the event labelled as "d" in Fig. 1). In order to analyze the reasons of the exceptional length of the observed ground motion, we analyze, in the following section, the snapshots of the wave field propagation across a North-South vertical section that encompass both the source of the event "d" and the neighborhoods of the station MODE.

Wave propagation across a vertical profile
A remarkable advantage of 3D numerical modeling consists in the possibility to visualize the wave phenomena, which 5 causes unexpected features in the observed ground motion. The most apparent anomalies in the observed ground motion that were reasonably well predicted with the 3D model are the PGV and the long duration at the stations south of the epicenter during the MW 4.1 12/6/2012 event (labeled with d in Fig. 1). In Fig. 10, we compare the horizontal components of the empirical and the computed time series at the stations T828, T824 and MODE. Even though we could not reach a match between the time series in a strict sense, the results obtained with the 3D model represent a significant improvement if 10 compared to the 1D model results.
In order to investigate the cause of the particular features of the ground motion in these stations we can follow the modelled propagation of the seismic waves on a vertical profile extracted from the 3D spatial domain (Figs. 11, 12 and 13), whose trace on the surface is shown by the red dashed line in Fig. 1 The grey shadow on the profiles represents the wave amplitude whereas the yellow lines represent the interfaces between the structural units. For the sake of clearness, the structural units are labeled in Fig. 11a only. The profile samples three different areas of the geological structure, the northern, the central and the southern. The northern area is characterized by the Ferrara Folds ( Fig. 1), where high velocity layers (VS =1.7 km/s) are lifted up to few tens of meters below the surface. The central 20 area is characterized by a deep syncline with thick, low velocity (VS < 1 km/s) superficial layers of sediments and alluvial deposits. In the southern part, we find the Emilia folds, which again reduce the thickness of the soil cover. In the first snapshot, taken after 2 s of propagation (Fig. 11a), we can see the initially concentric wavefronts propagating from the source located in the Ferrara folds. After 4 s (Fig. 11b) the wavefronts propagating towards south assume an almost plane shape, after having been deformed in the slower formations of the basin. We can clearly discern the compressional waves 25 (denoted by the letter P), and the shear waves (denoted by the letter S), being the latter stronger and slower, with a shorter wavelength. After 8 seconds (Fig. 12a) the direct S has reached T828 and we can observe how at that time the S wave is reflected from the soil surface above the ridge and channeled in the dipping layers south of it. Because of the layers dip, the reflected S wave hits the layers at a post-critical angle and generates a number of diffracted waves, which correspond to surface waves overtones, if we adopt a mathematically more elegant formalism. After 16 seconds (Fig. 12b), the 30 aforementioned diffracted waves can be well recognized in the profile across the wave-field and we can associate them with the strong phases following the direct S arrival at the 3 considered stations.
For example, the strong wave-train predicted at T824 between 16 and 20 s of propagation (Fig. 10b)

Code availability
The digital 3D geological model was set up with the commercial software GeoModeller 3.