## 1. Introduction

On April 6, 2009 a strong earthquake (M_{L} 5.9, M_{W} 6.3), hereafter called main shock, struck the Aterno Valley in the Abruzzo region (central Italy) causing heavy damage in L’Aquila and in several nearby villages and killing more than 300 people. The event had a pure normal faulting mechanism, with a rupturing fault plane NW striking and 45°SW dipping; hypocentral location was at 9.5 km depth and epicenter at a distance of about 2 km WSW from L’Aquila center [1]. Few days later, the aftershock activity involved also the area NE of L’Aquila toward Arischia and Campotosto. The overall distribution of the aftershocks defined a complex, 40 km long and 10-12 km wide, NW trending extensional structure. The largest damage was mainly distributed in a NW-SE direction [2], according to the orientation of the Aterno river valley. The area has a high seismic hazard level in Italy [3] and has experienced in the past destructive earthquakes such as the 1349, I=IX–X; the 1461, l'Aquila, I=X and the 1703, I=X [4]. Many active faults are recognized in the area and several of them are indicated as potential sources for future moderate and large earthquakes by several authors (see for a review [5]).

The Aterno river basin is a complex geological structure with a carbonate basement outcropping along the valley flanks and elsewhere buried below alluvial and lacustrine deposits with variable thickness. The surface geology is even more complicated by the presence at L’Aquila of breccias consisting of limestone clasts in a marly matrix. Such complex geological scenario reflected in a large spatial variability of amplitude and frequency content of the ground motion (e.g. [6]). Among several studies performed on the recorded events, there is not a note dedicated to the modeling of the shear seismic velocities of the crust structures, yet it is useful for the geological reconstruction and fundamental for computing seismograms. Simulation of the ground motion has been performed at L’Aquila based on literature data [7] by using the Neo-Deterministic Seismic Hazard Analysis (NDSHA) [8,9], an innovative modeling technique that takes into account source, propagation and local site effects. In order to estimate realistic ground motion we need physical parameters of rocks from surface to depths greater than the earthquake hypocenter. At engineering scale, microzoning activities promoted by the Italian Civil Defense Department [10] have performed V_{S} measurements at depths around 25 m, in gravelly soils with different degree of cementation, alternating to thin layers of finer deposits (sands and/or silts) that often include carbonate boulders (www.cerfis.it). The investigated depths are too shallow to define the vertical and lateral passage from soft sediments to rock basement (V_{S} at least of 800 m/s) which was sporadically found. At regional scale, a physical model is available extending to depths of about 300 km [11].

Aim of this paper is to retrieve V_{S} models of the shallow crust in the Aterno river valley from the non-linear inversion of the group velocity dispersion curves of the fundamental mode extracted with the FTAN method (e.g. [12,13,14]) from recordings of earthquakes with M_{L} ≥ 2.9 (Table 1) between April 5 and November 10, 2009, in the selected coordinate window of 42.4 ± 0.2 N and 13.4 ± 0.2 E. In addition, V_{S} of the superficial 30 m of Aterno alluvial soils are defined by an active seismic experiment in the Coppito area, and compared with nearby cross-hole measurements. The V_{S} profiles vs. depth are then attributed to lithotypes along a geological cross section from the epicenter to a seismic station at L’Aquila. Simulation of the main shock is performed with the NDSHA approach and the computed response spectra and the H/V spectral ratios are compared with those recorded.

## 2. Geological and geophysical setting

The epicentral area of the L’Aquila seismic sequence mainly corresponds to the upper and middle Aterno river valley which is bounded by predominantly NW-SE-striking and SW-dipping active normal faults (e.g. [15] and references therein) and characterised by the high variability of the geologic and geomorphologic patterns (Fig. 1). The valley is superimposed on a Quaternary lacustrine basin of tectonic origin and surrounded by carbonates. The thickness of the Quaternary deposits is variable, from about 60 m in the upper Aterno river valley to more than 200 m in the middle valley.

The upper and middle crustal structure of the Abruzzo region has been reconstructed from V_{P} and V_{P}/V_{S} images on a 3D grid with a node vertical spacing of 4 km, by using local earthquake tomography and receiver function modeling [16]. Low V_{p} velocities (smaller than 4.8 km/s) were found in correspondence with the main Plio-Quaternary basins (Fucino, l'Aquila and Sulmona basins) and high V_{P} velocities (V_{p} larger than 5.0 km/s) were mostly correspondent with the outcropping Mesozoic carbonates. Low V_{P} and high V_{P}/V_{S} anomalies were found beneath the L'Aquila and Fucino Quaternary basins between 4 and 12 km depth, suggesting the existence of fluid-saturated rock volumes. Very high P-wave (6.7–7.0 km/s) and S-wave (3.6–3.8 km/s) velocity bodies were observed below 8–12 km depth, and interpreted as deep crustal or mantle rocks exhumed before the sedimentation of the Mesozoic cover. This layer has a regional character as it is found at ~15 km of depth below the 1°x1° cell containing L’Aquila epicentral area, and lying on the mantle detected at ~35 km [11].

The analysis of strong and weak motion recordings in 1996-98 put in evidence amplification effect at low frequencies (0.6 Hz) in the town of L'Aquila and 2D numerical modeling allowed to fit it along a SW-NE section [17]. A sedimentary basin was inferred, filled by lacustrine sediments, with a maximum thickness of about 250 m, below the breccias formation about 50 m thick.

## 3. V_{S} models

### 3.1. Data analysis

The 2009 seismic sequence was recorded by Rete Accelerometrica Nazionale (RAN) network, managed by the Italian Civil Defense Department, some of which located at L’Aquila (AQK station) or in the NW of it (AQG, AQA, AQM, AQV), and by the station AQU operating since 1988 as part of the Mediterranean Network (MedNet), managed by the Italian Istituto Nazionale di Geofisica e Vulcanologia (INGV). They are equipped with three-component accelerometers set to 1 or 2 g full-scale, coupled with high resolution digitizers, while AQU is also equipped with a very broadband Streckeisen STS-1.

In order to obtain V_{S} models for the L’Aquila basin shallow crust, we have analysed about fourty earthquakes (M_{L} ≥ 2.9) recorded at the RAN and AQU stations, and rotated to get the radial and transverse component of motion. Rayleigh wave group velocities of the fundamental mode have been measured from vertical and radial components of 17 events (Fig. 1, Table 1).

As regards the shallow 30 m subsoil, the same analysis has been applied to recordings of an active seismic experiment performed in the Coppito area, about 500 m far from the AQV station (Fig. 1).

The group velocity is measured as function of period by the Frequency Time Analysis (FTAN) on single waveforms (e.g. [12,13,14]). The FTAN method allows to isolate the different phases in a seismogram, in particular the fundamental mode of surface waves. A system of narrow-band Gaussian filters is employed, with varying central frequency, that do not introduce phase distortion and give the necessary resolution in the time-frequency domain. The source-receiver distance is commonly assumed to be the epicenter distance when it is much greater than the event depth. When this assumption is not valid, in order to extract the correct dispersion curve of Rayleigh waves we have to add a time delay to seismograms as Δt=h/V_{P}, h being the source depth, V_{P} the average P-wave velocity from surface to the hypocenter, and then analyze the seismograms by FTAN, considering the hypocenter distance (e.g. [19] and references therein).

The dispersion curves obtained in such a way can be inverted to determine S-wave velocity profiles versus depth. A non-linear inversion is made with the Hedgehog method ([20,14] and references therein) that is an optimized Monte Carlo non-linear search of velocity-depth distributions. In the inversion, the unknown Earth model is replaced by a set of parameters and the definition of the structure is reduced to the determination of the numerical values of these parameters. In the elastic approximation, the structure is modeled as a stack of N homogeneous isotropic layers, each one defined by four parameters: V_{P} (dependent parameter), density (fixed parameter), V_{S} and thickness (independent parameters).

Given the error of the experimental phase and/or group velocity data, it is possible to compute the resolution of the parameters (parameter step), computing partial derivatives of the dispersion curve with respect to the parameters to be inverted ([14] and references therein). The theoretical phase and/or group velocities computed during the inversion with normal-mode summation are then compared with the corresponding experimental ones and the models are accepted as solutions if their difference, at each period, is less than the measurement errors and if the r.m.s. (root mean square) of the differences, at all considered periods, is less than a chosen quantity (usually 60–70% of the average of the measurement errors). Being the parameter step indicative of the parameter resolution, all the solutions of the Hedgehog inversion differ by no more than ±1 step from each other. A good rule of thumb is that the number of solutions is comparable with the number of the inverted parameters. From the set of solutions, we accept as representative solution the one with r.m.s error closest to the average r.m.s error of the solution set, and hence reduce, at the cost of loosing in resolution, the projection of possible systematic errors [20] into the structural model. Other selection criteria of the representative solution are discussed in detail by [21].

### 3.2. Results

Dispersion curves of Rayleigh wave fundamental mode have been extracted from the vertical and/or radial components of recordings of 17 events (Table 1, Fig. 1). Average dispersion curves have been computed along 12 paths and inverted with Hedgehog method [20,14] to get V_{S} models. A V_{P}/V_{S} ratio equal to 1.8 turned out to be a suitable value after a set of tests made varying it between 1.8 and 2.1. In other words, keeping all other values of the parameterization unchanged, the number of solutions maximizes for V_{P}/V_{S}=1.8. Moreover, the analysis of group velocity derivatives with respect to elastic parameters versus depth has allowed to define the sensitivity of the investigated periods on the S-wave velocity structure. Time corrections have been applied to the recordings by assuming the regional average V_{P} computed from the V_{S} model relative to the cell 1°x1° containing L’Aquila [11], assuming a V_{P}/V_{S} ratio of 1.8. Such value is in very good agreement with those used by INGV (5 km/s) and [1] (V_{P}=5.5 km/s) for earthquake location and attributed to the shallow 11.1 km of crust.

In the following, the results (dispersion data and Hedgehog V_{S} solutions) are presented by grouping the paths as northern paths (1-2-10-11-12), middle Aterno river valley paths (6-7-8), middle paths (3-5), and southern paths (4-9). The events are listed in Table 1 and the paths are located in Fig. 1. Moreover, the results obtained with the same methods from an active seismic experiment in the Coppito area, about 500 m far from the AQV station, are presented.

#### 3.2.1. Northern paths

All the results obtained for the northern paths (1-2-10-11-12) are shown in Fig. 2. The dispersion curves relative to path 1 have been extracted from the recordings of the 4/13/2009 event, located nearby Campotosto, at AQG, AQV and AQM stations. Rayleigh data, sampled at periods of 1-3.8 s, have been inverted in V_{S} models of the shallow 6 km of crust. The representative solution is characterized by velocities increasing from 0.9 to 3.0 km/s at 3.9 km of depth.

The path 2 is relative to the recording of the 6/22/2009 event at AQK station. Dispersion curves have been extracted at periods of 0.5-1.6 s and V_{S} models of the shallow 2 km have been retrieved. The representative solution is characterized by velocities increasing from 0.9 to 2.7 km/s at 1.3 km of depth.

The dispersion curves relative to path 10 have been extracted from the recordings of 2 events at AQU station. Rayleigh data have been sampled in the 0.4-1.2 s period range and the V_{S} models are relative to the shallow 2 km. The representative solution is characterized by shear velocities increasing from 0.9 km/s to 2.2 km/s at 0.8 km of depth.

The dispersion curves along the path 11 have been extracted from the radial and vertical components of the 04/13/09 event recording at AQU station. The average dispersion curve is defined in the 1.0-2.5 s period range and V_{S} models of the shallow 3 km have been retrieved from it. The representative solution shows velocities increasing from 0.9 km/s to 2.9 km/s at 1.8 km of depth.

The average dispersion curve along the path 12 has been computed from the dispersion curves extracted from the vertical components of 2 earthquakes on 04/09/09 and is defined at periods between 0.7 and 1.2 s. The retrieved V_{S} models are relative to the shallow 2 km and the representative solution is characterized by velocities increasing from 1.1 to 2.8 km/s at 1 km of depth.

The interpretation of the obtained V_{S} models, attributed beneath the middle points of the paths (bottom of Fig. 2), is performed by taking into account available geological data [18]. Outcropping rocks along the cross section A through the paths 12-1-11-10 consist of limestone and marls, which are found with a thickness of ~1 km in the Campotosto drilling (located in Fig. 1), below a 1.2 km thick layer of clays and sandstones. The Mesozoic carbonate horizon is found at the well bottom (2.45 km). V_{S} range between 0.9 and 1.5 km/s in the shallowest 100-200 m of weathered rocks, and reach the value of 2.5 km/s at depths between 0.6 and 1.7 km. The velocity of 2.5 km/s can be reasonably attributed to fractured limestone rocks and, based on the Campotosto drilling stratigraphy, to the top of the Mesozoic limestone horizon. Such horizon has V_{S} of ~2.8 km/s (V_{P} ~5 km/s) at depths of 1-2.5 km and of ~3 km/s (V_{P} ~5.5 km/s) at 4 km of depth (detected only below the path 1).

#### 3.2.2. Middle Aterno river valley paths

All the results obtained for the middle Aterno river valley paths (6-7-8) are shown in Fig. 3. The paths are relative to the recording of the 4/6/2009 mainshock event at AQK station (path 6), and of the 4/7/2009 events at the AQG station (paths 7 and 8).

As regards the path 6, dispersion curves have been extracted at periods of 0.9-2 s and V_{S} models of the shallow 1 km have been retrieved. The representative solution is characterized by velocities decreasing from 0.85 to 0.7 km/s in the shallow 0.1 km, and deeper increasing from 0.9 to 1.45 km/s at 0.4 km of depth. Dispersion curves relative to path 7 are defined in the period range of 0.9-1.8 s and V_{S} models of the shallow 1.3 km have been retrieved from their average curve. The representative solution presents velocities increasing from 0.6 to 2.1 km/s at 0.8 km of depth. The average dispersion curve relative to path 8, is defined in the period range of 0.4-1.4 s and V_{S} models of the shallow 1.3 km have been obtained. The representative solution presents velocities increasing from 0.5 to 1.7 km/s at 0.5 km of depth. Taking into account the geological data relative to the shallowest 0.2-0.3 km [10], stratigraphies may be attributed to the V_{S} profiles.

Lacustrine soils have a thickness increasing from about 0.2 km in the center of the valley (path 7) to about 0.4 km towards L’Aquila (path 6) with V_{S} increasing from 0.5 km/s to ~0.9 km/s. Maiolica and flysch layers have an average V_{S} of ~1.2 km/s while breccia has a V_{S} of ~0.9 km/s. The calcarenites with V_{S} of ~1.4 km/s deepen from 0.15 km to 0.45 km towards S (path 6).

#### 3.2.3. Middle paths

All the results obtained for the middle paths (3-5) are shown in Fig. 4. The path 3 is relative to the recording of the 4/11/2009 event at AQU station. It crosses the Aterno valley delimited on the eastern side by the Paganica fault. Dispersion curves have been extracted at periods of 0.7-1.9 s and V_{S} models of the shallow 2 km have been retrieved. The representative solution is characterized by velocities increasing from 0.7 to 2.0 km/s at 0.8 km of depth.

The path 5 is relative to the recording of the 4/9/2009 event at AQU station and crosses the lower border of the middle Aterno valley, on the west of L’Aquila. Dispersion curves have been extracted at periods of 0.8-1.4 s and the retrieved V_{S} models are relative to the shallow 2 km. The representative solution is characterized by velocities increasing from 0.6 to 2.5 km/s at 0.9 km of depth.

From the comparison of the representative V_{S} profiles (chosen Hedgehog solutions) below paths 3, 5 and 6, it turns out that the western sector is characterized by high velocities (~1.8 km/s) at very shallow depth (~0.1 km) which are detected at 0.8 km in the eastern sector and are not found in the investigated shallow 1 km in the epicentral area of the L’Aquila main shock.

#### 3.2.4. Southern paths

The paths 4 and 9 cross the south of L’Aquila consisting of outcropping alluvial deposits and carbonate rocks (Fig. 1). All the results are shown in Fig. 5. Both the paths are averaged on 3 events recorded at AQU station. Dispersion curves have been extracted at periods of 0.5-0.9 s, for the path 4, and of 0.6-1.2 s for the path 9. V_{S} models of the shallow 1 km have been retrieved which show homogeneous velocities increasing from 0.9 to 2 km/s at ~0.5 km of depth.

#### 3.2.5. Active seismic experiment

V_{S} models have been also retrieved from an active seismic experiment performed in the Coppito area (Fig. 1) with geophone offsets of 64 m and by using FTAN and Hedgehog methods. The models are characterized by an average velocity (V_{S30}) of 190 m/s in the shallow 30 m of alluvial soils (Fig. 6b). Instead a V_{S30} of 473 m/s is obtained from a cross-hole test, 500 m distant, performed close to the AQV station, in the same alluvial soils [22]. Such discrepancy has important consequences in the respect of the national building code as the soil classification changes from C to B [23].

The comparison of the frequency of the maximum peak of the H/V spectral ratio, relative to the main shock recorded at the AQV station, with the 1D spectral amplifications, computed with SHAKE program [24], by assuming the two different V_{S} data sets has evidenced the agreement with the V_{S} profiles relative to FTAN-Hedgehog measurements (Fig. 6c). Once more, such comparison evidences: 1) the strong lateral and vertical heterogeneities of such alluvial soils; 2) the cross-hole (and down-hole) point-like measurements, even though quite precise, may not be representative of the average seismic path (e.g. [25,26]). We remind that the surface measurements for the FTAN analysis do not need boreholes or conventional arrays, hence are particularly suitable in the urban areas.

## 4. Ground motion modeling

Simulations of the 2009 L’Aquila main shock have been performed with the Neo-Deterministic Seismic Hazard Analysis (NDSHA), an innovative modeling technique that takes into account source, propagation and local site effects [8,9].

This approach uses a hybrid method consisting of modal summation and finite difference methods. The path from the source up to the region containing the 2-D heterogeneities is represented by a 1-D layered anelastic structure. The resulting wavefield for both SH- and P-SV- waves is then used to define the boundary conditions to be applied to the 2-D anelastic region where the finite difference technique is used. Synthetic seismograms of the vertical, transverse and radial components of ground motion are computed at a predefined set of points at the surface. Spectral amplifications are computed as response spectra ratios, RSR, i.e. the response spectra computed from the signals synthesized along the laterally varying section (2D) normalized by the response spectra computed from the corresponding signals, synthesized for the bedrock (1D). A scaled point-source approximation ([27] as reported in [28]) has been considered to scale the seismogram to the desired scalar seismic moment.

Modeling of the main shock has been done along a geological cross section at L’Aquila, from the epicentre through the AQK station (Fig. 7).

Along the cross section, the outcropping units are represented by megabreccias except in the Aterno river, where recent fluvial sediments are present. At engineering scale, after the 2009 seismic sequence, geological and geophysical studies, beside several drillings with down-hole tests generally reaching depths around 25-30 m, have been performed at L’Aquila to reconstruct the shallow 200-300 m of subsoil [10]. The new interpretation of available gravity data indicated a meso-cenozoic carbonate unit (average density=2.6 g/cm^{3}) lying at maximum depth of 200-300 m in the Aterno valley, below flysch or breccia unit (average density=2.4 g/cm^{3}) and Quaternary products including alluvial and lacustrine deposits and fan alluvial material (average density=1.9 g/cm^{3}). Strong lateral and vertical geological heterogeneities have been evidenced which, in the center of L’Aquila, are mainly due to a discontinuous stiff top layer of breccias, called megabreccias, overlying soft lacustrine sediments. The geometry of the vertical and lateral passage of the lacustrine to megabreccia deposits is still poorly known due to the shallow drillings. Taking into account the available geological cross sections [10], a computing cross section has been prepared (Fig. 8). We have attributed V_{S} of 0.2 km/s to the thin silt deposits according to surface measurements at Coppito (Fig. 6). Taking into account the V_{S} profile obtained for the path 6 (Fig. 3), we have assigned velocities of 0.9 km/s both to megabreccias and sandstones, and of 1.45 km/s both to marls and calcarenites. We are in agreement with [17] as regards velocities of megabreccias, instead a strong discrepancy results for velocities of calcarenites (our 1.45 km/s against 2.5 km/s). Literature V_{S} [17] have been attributed to Aterno river recent deposits (colluvium and alluvial deposits) and to upper and lower lacustrine soils.

The 1-D reference model has been chosen according to the regional model [11]. A parametric study has been done for the best dip angle between 43° and 60°. A dip of 56° turned out to be able to fit the observed response spectra and is in agreement with the geometry of the seismogenic fault [29]. A small valley of colluvium had to be hypothesized beneath the AQK station, realistically assumed on the basis of geological considerations, in order to fit the observed response spectra and the frequency of the main peak of the H/V spectral ratio obtained from the main shock recording at the AQK station (Fig. 9). A good fit results between observed and synthetic response spectra despite the fictitious greater distance (6 km) of the cross section in the computation of the ground motion for a point source.

Acceleration time series for P-SV and SH-waves, including surface waves, have been computed along the cross section (Fig. 10), using a grid spacing of 3 m since the lowest seismic velocity is 200 m/s and 10 points are requested for a good sampling of the minimum wavelength corresponding to 7 Hz frequency. They are shown at an array of sites along the cross section, with 50 m spacing. The alluvial and colluvial soil cover is responsible for the amplification of the peak values and duration of accelerations along the radial and transverse components. This amplification is higher along the vertical component.

Spectral amplifications computed for the vertical component of ground motion show maximum values of 10 at frequencies lower than 1 Hz, in correspondence of the layer of megabreccias and of about 5 at 4 Hz in correspondence of the Aterno river alluvial sediments (Fig. 11). As regards the radial and transverse components, spectral amplifications of 2-3 are computed for a wide frequency range (1-7 Hz), along the cross section. Taking into account that the majority of the buildings, generally 2-5 floor, at the historical center of L’Aquila suffered serious damage, we can argue that structures lying on soils suffered amplifications of 2-3 along the horizontal components and up to 5 along the vertical component of the ground motion.

## 5. Conclusions

A realistic estimation of the ground motion at L’Aquila for the M_{W} 6.3 earthquake is obtained by the NDSHA approach, an innovative modeling technique that takes into account source, propagation and local site effects [8,9]. A key point is the definition of V_{S} models representative of the seismic path, like those obtained from the non-linear inversion of Rayleigh group velocities of the fundamental mode extracted with the FTAN method from earthquake recordings and active seismic surveys. Very fractured carbonatic rocks with V_{S} of ~1.4 km/s, covered by alluvial soils with a maximum thickness of ~0.4 km in the center of the Aterno valley are retrieved. The top of the carbonates with the average velocities of 5 and 2.8 km/s, compression and shear respectively, lays at 2-2.5 km of depth and rises to 1 km in the NW part of the valley. This result contradicts available seismic [17] and gravity [10] modeling of the carbonate horizon with V_{S}=2.5 km/s and density ρ =2.6 g/cm^{3} at some hundred meters of depth. The carbonate horizon with V_{S} of ~3 km/s is found at 4 km of depth. Moreover, the shallowest 30 m of alluvial soils have average V_{S} of ~0.2 km/s against ~0.5 km/s as obtained from cross-hole measurements at the AQV station, about 500 m distant. Such velocity difference evidences that strong lateral and vertical geological heterogeneities are present and that the cross-hole (and down-hole) point-like measurements, even though quite precise, may not be representative of the average seismic path.

The soundness of the synthetics is in the good fitting of the recorded H/V spectral ratio and response spectra, despite the point-source approximation. The lateral and vertical geological variability mainly due to the covering of megabreccias on soft soils are responsible of spectral amplifications, mostly for the vertical component, for a wide frequency range (0.5-7 Hz). Taking into account that the majority of the buildings, generally 2-5 floor, at the historical center of L’Aquila suffered serious damage, we can argue that spectral amplifications might have been responsible for damage, beside the near-field conditions. This study shows that realistic ground motion can be computed in advance for the several active faults of the L’Aquila district and that a sounded building code can be formulated for the restoration of the existing damaged buidings and for new building design.