# Broad-band strong motion simulations coupling k-square kinematic source models with empirical Green's functions: the 2009 L'Aquila earthquake

Sergio del Gaudio, Mathieu Causse, Gaetano Festa

► To cite this version:

Sergio del Gaudio, Mathieu Causse, Gaetano Festa. Broad-band strong motion simulations coupling k-square kinematic source models with empirical Green's functions: the 2009 L'Aquila earthquake. Geophysical Journal International, 2015, 203 (1), pp.720 - 736. 10.1093/gji/ggv325. hal-01692245

HAL Id: hal-01692245

https://hal.science/hal-01692245

Submitted on 8 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Broad-band strong motion simulations coupling k-square kinematic source models with empirical Green's functions: the 2009 L'Aquila earthquake

Sergio Del Gaudio,<fnref n="1" />,<fnref n="2" />,* Mathieu Causse<fnref n="3" /> and Gaetano Festa<fnref n="2" />

Accepted 2015 July 30. Received 2015 July 29; in original form 2014 October 1

## SUMMARY

The use of simulated accelerograms may improve the evaluation of the seismic hazard when an accurate modelling of both source and propagation is performed. In this paper, we performed broad-band simulations of the 2009, M 6.3 L'Aquila earthquake, coupling a $k^{-2}$ kinematic model for the seismic source with empirical Green's functions (EGFs) as propagators. We extracted 10 EGFs candidates from a database of aftershocks satisfying quality criteria based on signal-to-noise ratio, fault proximity, small magnitude, similar focal mechanism and stress drop. For comparison with real observations, we also derived a low-frequency kinematic model, based on inversion of ground displacement as integrated from strong motion data. Kinematic properties of the inverted model (rupture velocity, position of the rupture nucleation, low-frequency slip and roughness degree of slip heterogeneity) were used as constraints in the $k^{-2}$ model, to test the use of a single specific EGF against the use of the whole set of EGFs. Comparison to real observations based on spectral and peak ground acceleration shows that the use of all available EGFs improves the fit of simulations to real data. Moreover the epistemic variability related to the selection of a specific EGF is significantly larger (two to three times) than recent observations of between event variability, that is the variability associated with the randomness of the rupture process. We finally performed 'blind' simulations releasing all the information on source kinematics and only considering the fault geometry and the magnitude of the target event as known features. We computed peak ground acceleration, acceleration Fourier and response spectra. Simulations follow the same trend with distance as real observations. In most cases these latter fall within one sigma from predictions. Predictions with source parameters constrained at low frequency do not perform better than 'blind' simulations, showing that extrapolation of the low-frequency description of the rupture front as inferred by the kinematic inversion may introduce some bias in the final simulations.

**Key words:** Earthquake ground motions; Earthquake source observations; Wave propagation.

## INTRODUCTION

Seismic hazard is commonly assessed by the use of ground motion prediction equations, appropriately calibrated to tackle specific scenarios of future earthquakes. As well, the dynamic response of the structures is frequently computed using real accelerograms as input, eventually scaled to be compatible with the design spectrum or with the target magnitude of the event. Nevertheless ground motion prediction equations do not represent the variety of source, propaga-

tion and soil conditions and this can lead to incorrect estimates for some specific case studies, such as in the near-fault range (Hall et al. 1995; Somerville 2003), for large crustal events (Chiou et al. 2008), for deep basins, where long period waves are trapped (Miyake & Koketsu 2005; Koketsu & Miyake 2008) and for stations located on the hanging wall of the fault (Donahue & Abrahamson 2014). In these conditions, possible biases may arise from the scarcity of data and/or an inadequate physical modelling of the source and wave propagation.

In such cases ground motion simulations can be a valid tool to complement prediction equations for scenario studies, provided that both source and propagation are accurately described and uncertainties properly addressed. In the low to intermediate frequency range

<footnote n="1">¹AMRA s. c. a r. l., Napoli, Italy. E-mail: sergio.delgaudio@unina.it</footnote>

<footnote n="2">²Department of Physics, Università di Napoli Federico II, Napoli, Italy</footnote>

<footnote n="3">³ISTerre, Université Joseph Fourier, IFSTTAR, Grenoble, France</footnote>

<footnote>*Now at: Institut de Radioprotection et de Sûreté Nucléaire (IRSN), Fontenay aux Roses, France.</footnote>
($<$1–2 Hz) numerical simulations are commonly used for gener-
ating ground motion synthetics. Depending on the knowledge of
the velocity structure and the required computational costs, they
can be based either on a 1-D (e.g. Gallovič & Brokešová 2008;
Ameri et al. 2011b) or a 3-D description of the geological struc-
ture. In the latter case, finite differences (Graves 1996; Olsen 2000;
Graves et al. 2008), spectral element (Komatitsch & Vilotte 1998;
Stupazzini et al. 2009; Cuppillard et al. 2012) or finite elements
(Xu et al. 2003; Oglesby & Mai 2012) accurately model the wave
propagation in heterogeneous media, with the inclusion of topog-
raphy, basins and steep geological interfaces inside the Earth. Such
simulations can rely on either a kinematic or a dynamic description
of the earthquake source.

Alternatively, the impulse response of the medium can be re-
placed with the record of a small earthquake, referred to as em-
pirical Green's function (EGF), under the assumption of similarity
with the main event in terms of location, focal mechanism and stress
drop. Since the pioneering work of Hartzell (1978), many applica-
tions have made use of the EGFs for ground motion simulation
(Hutchings 1994; Irikura & Kamae 1994; Hartzell et al. 1999;
Causse et al. 2008) and kinematic source inversion (Miyake et al.
2003; Nozu & Irikura 2008). An EGF behaves as an effective
Green's function up to its own corner frequency, beyond which the
finite source effects of the small event start to appear (Hutchings
1994). Such a frequency limit could be in principle also extended
beyond the small event corner frequency, correcting the EGF by a
specific source model (Hutchings et al. 2007). Hence, it properly
describes the viscoelastic propagation of seismic waves from the
source to the receiver. The main advantage of using EGFs is that
the maximum frequency at which the Green's function is accurate
is larger than the one from numerical simulations, provided that
EGF magnitude is small enough to assimilate it to the medium
Green's function up to higher frequencies (Hutchings 1994). Al-
though the availability of high-quality EGFs is significantly in-
creasing for several regions, thanks to the installation of seismic
networks for monitoring the microseismicity and the aftershock ac-
tivity, this technique suffers from several limitations. In many cases
location of the small events does not completely cover the fault plane
and interpolation and distance correction are required to build up
a complete dataset of EGFs to be used for ground motion com-
putation, reducing the accuracy of the method (Hutchings & Wu
1990). Moreover, aftershocks are generally located in the regions
of the fault plane that provided no or little coseismic slip (Das &
Henry 2003). Also, estimation of stress drops through corner fre-
quency measurements is affected by large uncertainty (Cotton et al.
2013), introducing epistemic uncertainties in the selection of spe-
cific EGFs for simulations. Moreover, the magnitude of the EGFs
cannot be chosen too small, so as to maintain a good signal-to-noise
ratio in a broad frequency range. To improve the low frequency con-
tent of synthetics, it is possible to complete EGFs with numerical
waveforms computed in 1-D or 3-D models (e.g. Causse et al.
2009). Finally, EGFs do not account for possible nonlinear site ef-
fects in the shallow layers, although a first order 1-D correction
can be derived by elastic deconvolution to obtain the signal at the
bedrock convolved subsequently to include site effects (Heuze et al.
2004).

In order to reproduce ground motion observations, low fre-
quency kinematic models are commonly obtained from inversion
of strong motion data. Thus, rupture models for several moderate
to large magnitude events have already been inferred and collected
in a database (SRCMOD, Mai & Thingbaijam 2014). The upper
limit of the modelled frequency is however limited due to partial

knowledge of the velocity structure (<0.5 Hz), and hence smaller
than the frequency range of interest for engineering applications.
Moreover, kinematic inversion is well known to be an ill-posed
problem and differences in solutions may arise from data selec-
tion and processing, inversion strategies, regularization and Green's
functions computation (Mai et al. 2007; Festa & Zollo 2012; Shao &
Ji 2012). Array back-projection methods may track the rupture evo-
lution at higher frequencies (up to 1–2 Hz), but interpretation of
imaged energy on the fault in terms of physical quantities is still
debated (Ishii et al. 2005; Festa & Zollo 2006; Maercklin et al.
2012). To complement kinematic models at high frequency, the
self-similar rupture model of Herrero & Bernard (1996) is widely
used. With this description of the rupture, the slip amplitude de-
cays in the wavenumber domain **k** as k⁻². This model correctly
reproduces the high frequency f⁻² decay of the far-field waveform
displacement spectra, providing that the rise time decreases as k⁻¹
in the wavenumber domain and that the rupture velocity is constant
(Gallovič & Brokešová 2004). Additionally, to prevent overestima-
tion of ground motion amplitudes for directive stations, incoher-
ence in the rupture evolution has been introduced by a small scale
random propagation (Gallovič & Brokešová 2007; Causse et al.
2009).

In this paper, we perform simulations for the 2009, M<sub>w</sub> 6.3,
L'Aquila event, which occurred on April 6 (01:32 UTC), along a
normal fault. It struck the central area of the Abruzzo region in
Central Italy, heavily damaging the city of L'Aquila and its sur-
roundings, with several causalities and huge losses to the artistic
and cultural heritage. The hypocentre of the event was located at
around 9 km depth beneath the city of L'Aquila, and the rupture
activated a 50° dipping and 135° NW trending normal fault with a
length of about 20 km. This fault has the same geometrical charac-
teristics as the typical extensional structures of central Apennines
(Chiarabba et al. 2009; Chiaraluce et al. 2011). Different studies
pointed out the complexity of the earthquake nucleation and propa-
gation (e.g. Cirella et al. 2009, 2012). The initial stage of the rupture
history showed a predominant updip propagation which produced
large velocity peaks clearly observable on the three components
of the stations located in L'Aquila city. The subsequent southward
propagation along strike generated a second, larger scale directivity
effect that can be observed in the intensity map of the area and on
the data at stations located in the southern part of the region (Ameri
et al. 2011a).

The strong motion data recorded in the near-fault region exhibit
large amplitudes, localized peaks and large space variability, this
latter being also observed in the pattern of the damage distribution.
Peak ground accelerations reach values of 30–40 per cent of g in or
nearby L’Aquila city.

Beyond the main event, almost 3000 earthquakes with local mag-
nitude M_L > 1.9 occurred in the area during 2009 that were recorded
and located (Chiaraluce et al. 2011). These events consist of a fore-
shock sequence that started in January 2009 and culminated with a
M_w 4 event that occurred on March 30 (13:38 UTC) and an after-
shock activity that occurred with a 15-km-thick upper crust. This
earthquake hence represents an interesting case study from which
we can investigate near fault effects and understand the role of the
source and of the large scale wave propagation on the observed
ground motion and its spatial variability. Moreover, records of af-
tershocks and foreshocks provide a database of EGFs that can be
used as the impulse response of the medium in the ground motion
simulations. Specifically, records at the accelerometric network al-
low for the calibration of source parameters by comparison with
real observations, whilst records at the temporary network that was
installed after the main event, allow for an estimation of the ground motion in sites where no seismic data from the main event are available.

In this paper we compute synthetic waveforms for the L'Aquila earthquake assuming different *a priori* levels of knowledge for source properties (fault, magnitude, hypocentre location, low-frequency rupture kinematics, etc.), to be compared with real observations. In the first section, we summarize the methodology used to compute broad-band ground motion synthetic waveforms. We describe the hypotheses under which the displacement can be written as a summation of EGFs and we recall the main ingredients of the $k^{-2}$ source model that are used in this work. Then, we present the analysis of the aftershock database, from which we extracted a subset of EGFs and we discuss their variability. Later, we derive a low-frequency kinematic model from inversion of strong motion data. We next use this kinematic model to constrain the rupture process and perform broad-band ground motion predictions of the L'Aquila earthquake to analyse differences and advantages in using one single EGF as compared to the whole dataset of EGFs. Finally, we perform simulations releasing all of the information about the rupture description and assuming instead *a priori* distributions of source parameters (referred to as 'blind' predictions), and we compare the results to real observations both in terms of PGA and spectral amplitudes.

# BROAD-BAND SIMULATION METHOD

## Empirical Green's functions (EGFs) summation

Simulations are performed by discretely solving the representation integral in the form (Aki & Richards 2002)

$$u_n(\mathbf{r}, t) = \int_{\Sigma} s_i(\vec{\xi}, t) * [c_{ijkl}(\vec{\xi})v_j(\vec{\xi}) \times G_{nk,l}(\mathbf{r}, \vec{\xi}; t)] d\Sigma(\vec{\xi}), \quad (1)$$

where $u_n$ is the nth component of the displacement at the receiver located at position $\mathbf{r}$ and at time $t$, $\mathbf{s}$ is the slip vector, $\mathbf{c}$ is the elastic coefficient matrix, $\vec{v}$ is the normal to the fault, $\mathbf{G}$ is the Green's function and $\Sigma$ is the fault surface. Summation over repeated indices is assumed and the symbol $*$ denotes time convolution. The contribution $c_{ijkl}v_j G_{nk,l}$ is the nth component of the displacement produced by a unit slip in the ith direction, localized in space and time. Assuming a uniform stress drop and focal mechanism, we replaced this contribution with the displacement record of a small earthquake $\tilde{\mathbf{G}}(\mathbf{r}, \vec{\xi}; t)$ having the same location as the Green's function, sharing the same focal mechanism and stress drop as the target event, after multiplication by the shear modulus $\mu$ and normalization by the seismic moment of the small event $M_{\tilde{\mathbf{G}}}$:

$$u_n(\mathbf{r}, t) = \int_{\Sigma} \frac{\mu(\vec{\xi})}{M_{\tilde{\mathbf{G}}}^*(\vec{\xi})} S(\vec{\xi}) * \overline{G}_n(\mathbf{r}, \vec{\xi}; t) d\Sigma(\vec{\xi}) \quad (2)$$

$\tilde{\mathbf{G}}(\mathbf{r}, \vec{\xi}; t)$ is generally referred to as EGF. We also removed the index for the slip vector because we assumed a fixed rake direction along the fault plane, so that $S$ only represents the amplitude of the slip. The main advantage of using EGFs is that they naturally contain information about propagation and site effects between the source and the receivers, under the assumption of linearity of the soil response.

The size of the EGF controls the space discretization of the fault plane. Assuming self-similarity between the rupture of small and large events (Brune 1970), the value $N$ is

$$N = \frac{L}{l} = \frac{W}{w} = \left(\frac{M_0}{M_{\tilde{\mathbf{G}}}}\right)^{1/3}, \quad (3)$$

$N$ represents the scale factor in space between the main shock and any associated EGF. In the above formula $L$, $W$ and $M_0$ are the length, width, and seismic moment of the main event, respectively, while $l$, $w$ and $M_{\tilde{\mathbf{G}}}$ represent the same quantities for the small event. In our study, the main event rupture has a larger extent along strike than along dip. We then used squared subfaults with size $\Delta x = w$. After discretization, the displacement can be written as a summation of EGFs

$$u_n(\mathbf{r}, t) = \sum_{ij} \frac{\mu^{(ij)} \Delta x^2}{M_{\tilde{\mathbf{G}}}^{(ij)}} S^{(ij)} * \tilde{G}_n^{(ij)}(\mathbf{r}, t). \quad (4)$$

Here the indices $i$ and $j$ are related to subfault $(i,j)$ and run over strike and dip directions, respectively. Note that eq. (4) holds if the various $\tilde{G}_n^{(ij)}(\mathbf{r}, t)$ have the same seismic moment and stress drop (and hence the same corner frequency), and is valid at frequencies smaller than the EGF corner frequency $f_c$.

## Coupling EGFs with kinematic source models

To generate high frequency source energy, we adopted a kinematic description of the rupture based on the $k^{-2}$ source model, which correctly reproduces the high frequency $f^{-2}$ decay of the far-field waveform displacement spectra, using a rise time which decreases as $k^{-1}$ in the wavenumber domain and a constant rupture velocity (Gallovič & Brokešová 2004; Herrero & Bernard 1996). The amplitude of the static slip distribution is described in the wavenumber domain by an asymptotic $k^{-2}$ decay beyond the corner wavenumber $k_c = K/L_c$, where $L_c$ is the characteristic fault length and $K$ a non-dimensional parameter. The quantity $L_c/K$ represents the correlation length of the static slip distribution. Large $K$ values result then in high amplitude slip patches of small size, while small $K$ values give rise to a single smooth slip patch. As such, $K$ controls the roughness degree of slip heterogeneity on the fault. In order to simulate an earthquake with the same stress drop as the EGF, Causse et al. (2009) showed that $K \approx 0.74$. While the slip amplitude follows a $k^{-2}$ decay, the phase spectrum is fixed at small wavenumbers ($k \sim 1/L_c$) to localize the maximum slip at the centre of the fault plane or to match the slip distribution with a specific low frequency source model. At larger wavenumbers ($k \sim 1/L_c$), phases are randomly chosen, giving rise to a variety of static slip distributions. Furthermore, the rupture evolution on the fault is modelled as a slip pulse of width $L_0$ propagating at constant rupture velocity $v_r$ from the hypocentre. Following Bernard et al. (1996), the rise time is also assumed to have a $k$ dependency, being fixed at $\tau = L_0/v_r$ for $k < 1/2L_0$ and decaying as $k^{-1}$ for larger wavenumbers. Gallovič & Brokešová (2004) showed that in the case of a line fault embedded in a homogeneous medium, such source model generates far-field ground motion with a corner frequency

$$F_c = v_s C_d K / L_c, \quad (5)$$

where $C_d = \frac{v_s}{v_s - v_r \cos \theta}$ is the directivity coefficient (Ben-Menahem 1961), $\theta$ denotes the angle between the rupture front propagation and the source-receiver direction (the directivity angle), and $v_s$ the shear wave velocity nearby the source. Finally, small scale slip heterogeneities (with $k > 1/2L_0$) are assumed to have random rupture
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Location of seismic events (black dots) recorded in the first month after the main shock, whose epicentre is represented by a red star. Stations that properly recorded the main event and some of the EGFs are represented with a blue triangle, those that only recorded some of the EGFs are represented by a cyan triangle, and those that only recorded the main event are represented by a blue square. Green stars point to the selected EGF epicentres while the red rectangle represents the surface projection of the considered fault plane.</figcaption></figure>

directions (Causse *et al.* 2009). Rupture becomes then incoherent at frequencies larger than $v_r/L_0$. This choice strongly reduces the overestimation of directivity effect at high frequency (Herrero & Bernard 1994). In the following, source parameters (slip roughness parameter *K*, rise time $\tau$, rupture velocity $v_r$, position of the hypocentre) will be either deterministically defined (i.e. fixed from a reference source model obtained from low frequency kinematic inversion), or randomly selected according to the strategy described in Causse *et al.* (2009) (i.e. chosen from a priori distributions).

As abovementioned, the EGF summation scheme of eq. (4) is valid up to the EGF corner frequency $f_c$. Beyond $f_c$, EGFs cannot be considered as point sources and summation becomes incoherent. This leads to several numerical artefacts such as an overestimation of the high frequency source energy (e.g. Pavic *et al.* 2000; Kohrs-Sansorny *et al.* 2005) or appearance of artificial peaks due to the regular space-time slip discretization (e.g. Bour & Cara 1997). In order to perform broad-band simulations beyond $f_c$, we then adapted the EGF summation scheme according to the procedure proposed by Causse *et al.* (2009): (1) the number *N* of EGFs to sum up on each subfault is increased by a factor of $\gamma = (\alpha(N)/\beta)^{N/2}$, where $\alpha(N) = 2\sqrt{\ln(N/4)}$ and $\beta \approx 3.5$; (2) the acceleration amplitude spectrum is finally divided by $\gamma$ to preserve the seismic moment of the simulated event. Note that the spatial discretization of the fault plane remains unchanged. This procedure guarantees that the $k^{-2}$ source energy level is properly modelled above $f_c$.

# EGFs SELECTION

## Selection criteria

The 2009 L'Aquila earthquake provided a very large dataset of aftershocks and foreshocks that can be used as the impulse response

of the medium in the ground motion simulations. We started from a database of more than 400 events recorded in the first month after the main shock, whose location is shown in Fig. 1, together with the station positions. We then extracted a subset of events that satisfy the following quality criteria: (1) a large signal-to-noise ratio, (2) proximity to the fault plane of the main event, (3) focal mechanism similar to the one of the main shock and (4) similar magnitude range. Specifically, we screened out events with signal-to-noise ratio smaller than 3 in the frequency range 0.2–20 Hz. To compute the signal-to-noise ratio, we selected 5 s of record before the first *P* arrival time to be used as noise and 5 s around the maximum peak to peak amplitude to be used as signal. The ratio was measured on the amplitude spectra of the signal and the noise, frequency by frequency. We considered only events with distance from the fault plane smaller than 500 m, the locations being provided by Chiaraluce *et al.* (2011). The fault plane definition is described in the section ‘Deterministic source model using kinematic inversion’. We ruled out events with focal mechanisms significantly different from the one of the main event (differences of strike larger than 30°, differences in dip larger than 15°), based on focal mechanisms computed by Scognamiglio *et al.* (2010). In addition, we selected only events with moment magnitude ranging from 3.3 to 3.7. We finally ended up with a suite of 10 aftershocks, whose characteristics are reported in Table 1.

## Computation of seismic moment and stress drop

To compute the moment magnitude and the corner frequency of the selected EGFs and the main shock, a two-step inversion of displacement spectra was performed. We manually picked the *S* wave on the records at all stations and we selected a time window containing the whole direct *S*-wave train. The duration of the time window also depends on the magnitude and was fixed to 5 s for the
<table><caption>**Table 1.** List of the earthquakes used in this work. The first event corresponds to the main shock. Parameters in the columns represent local magnitude $M_L$ and moment magnitude $M_w$, corner frequency $f_c$ and associated uncertainty ($\sigma_{fc}$), stress drop $\Delta\sigma$, strike, dip and rake from focal mechanisms. Local magnitude is reported from INGV bulletin while moment magnitude and corner frequency were computed from inversion of displacement spectra.</caption><thead><tr><th>#EGF</th><th>Date-Time</th><th>M<sub>L</sub></th><th>M<sub>w</sub></th><th>f<sub>c</sub> (Hz)</th><th>Δσ<sub>fc</sub></th><th>Δσ (MPa)</th><th>Strike</th><th>Dip</th><th>Rake</th></tr></thead><tbody><tr><td>1</td><td>2009.04.06.01.32.40</td><td>5.9</td><td>6.0</td><td>0.29</td><td>0.26</td><td>8.9</td><td>133</td><td>54</td><td>-102</td></tr><tr><td>2</td><td>2009.04.06.04.47.54</td><td>3.6</td><td>3.9</td><td>1.7</td><td>1.6</td><td>1.0</td><td>129</td><td>59</td><td>-123</td></tr><tr><td>3</td><td>2009.04.06.10.36.18</td><td>3.5</td><td>3.5</td><td>2.9</td><td>1.7</td><td>1.3</td><td></td><td></td><td></td></tr><tr><td>4</td><td>2009.04.06.22.47.14</td><td>3.7</td><td>3.6</td><td>3.8</td><td>3.3</td><td>4.9</td><td>104</td><td>58</td><td>-120</td></tr><tr><td>5</td><td>2009.04.08.03.00.35</td><td>3.6</td><td>3.8</td><td>2.6</td><td>3.3</td><td>2.6</td><td>131</td><td>73</td><td>-155</td></tr><tr><td>6</td><td>2009.04.08.11.35.57</td><td>3.5</td><td>3.3</td><td>4.5</td><td>4.5</td><td>3.3</td><td>107</td><td>79</td><td>-152</td></tr><tr><td>7</td><td>2009.04.10.15.46.18</td><td>3.5</td><td>3.4</td><td>8.1</td><td>4.6</td><td>23.3</td><td>120</td><td>60</td><td>-44</td></tr><tr><td>8</td><td>2009.04.12.09.48.59</td><td>3.5</td><td>3.4</td><td>5.1</td><td>3.9</td><td>5.8</td><td>122</td><td>74</td><td>-150</td></tr><tr><td>9</td><td>2009.04.13.19.17.58</td><td>3.5</td><td>3.5</td><td>2.9</td><td>2.3</td><td>1.2</td><td>130</td><td>61</td><td>-120</td></tr><tr><td>10</td><td>2009.04.21.15.44.36</td><td>3.7</td><td>3.5</td><td>2.8</td><td>2.5</td><td>1.2</td><td>139</td><td>54</td><td>-107</td></tr><tr><td>11</td><td>2009.04.30.13.01.01</td><td>3.5</td><td>3.4</td><td>6.7</td><td>5.2</td><td>11.0</td><td>131</td><td>60</td><td>-128</td></tr></tbody></table>

EGFs and 20s for the main event; in both cases it starts 1s before the S wave pick. The displacement spectra were computed for each component in the frequency domain and the analysed spectrum was obtained as the modulus of the spectral amplitude computed on the horizontal components. For the inversion, we assumed in the angular frequency domain $\omega$ an $\omega^{-2}$ decay of the displacement spectrum $D(\omega)$ beyond its corner frequency $f_c$:

$$D(\omega) = \frac{M_0}{1 + (f/f_c)^2} \Gamma(\omega, R), \quad (6)$$

where $M_0$ is the seismic moment and $\Gamma$ an attenuation function that depends on $\omega$ and the hypocentral distance $R$:

$$\Gamma(\omega, R) = \frac{R_p F_s}{\beta_s^{5/2} \beta_r^{1/2} \rho_s^{1/2} \rho_r^{1/2}} \frac{1}{R} e^{-\frac{\omega T(R)}{2Q}}. \quad (7)$$

In the above formula, $R_p$ is the radiation pattern coefficient, $F_s$ is the free surface coefficient, $\beta$ is the S-wave velocity, $\rho$ is the density, $Q$ is the anelastic attenuation quality factor and $T(R)$ is the traveltime between the source and the receiver. Finally, subscripts s and r refer to the source and the receiver, respectively. We used $R_p = 0.63$ and $F_s = 2$ as suggested in Zollo et al. (2014), while other quantities are computed according to the velocity model of Bianchi et al. (2010).

We found a strong correlation between the corner frequency $f_c$ and the attenuation factor $Q$ for the selected dataset. For retrieving the source parameters, we then fixed the quality factor to a constant value $Q = 100$, this value representing a good compromise for the attenuation in the analysed frequency and distances ranges (Bindi et al. 2009). With this condition, the inverse procedure determines the plateau level and the corner frequency of the events, using the Levenberg-Marquardt algorithm (Zollo et al. 2014). Assuming a log-normal distribution for both the seismic moment and the corner frequency, we estimated the source parameters as the geometrical mean of the single estimates provided by the different stations (Zollo et al. 2014). Uncertainties were computed using the confidence limits at 68.8 per cent, being directly extracted from the distribution of the estimated values. The threshold in the confidence limit was selected in order to compare this choice for the uncertainty with the classical error estimation based on standard deviation. The retrieved values of moment magnitude and associated corner frequency are shown in Table 1. The associated stress drops were computed using the Brune source model (Brune 1970) and are also listed in Table 1. Fig. 2, represents the corner frequency as a function of the seismic moment in a logarithmic scale for the analysed events. Stress drop for the main event is 8.9 MPa, which is compatible with the value of 9.6 MPa retrieved by Bindi et al. (2009).

Corner frequencies computed for small events at single stations show a large dispersion, indicating that the associated estimation of the stress drop is affected by large uncertainties. For the selected events, we found that the stress drop ranges between 1 and 20 MPa, but uncertainties are so large that the stress drop value obtained for the main event falls within one sigma from the mean value of the stress drops for all the small earthquakes. Different values in the stress drop may lead to significant differences in the ground motion amplitude, increasing the epistemic variability associated with the selection of a specific EGF. In order to generate a catalogue of EGFs with the same stress drop and seismic moment, we applied a correction to each EGF according to Hutchings et al. (2007). The EGF amplitude spectra were first divided by their corresponding Brune spectra, considering the seismic moment and stress drop values of Table 1. The spectra were next multiplied by a Brune spectrum associated to a reference seismic moment $M_0 = 2.23 10^{14}$ N m ($M_w$ 3.5) and stress drop $\Delta\sigma = 3$ MPa, this latter value being the geometric mean of the EGF stress drops. This new set of EGFs has the same amplitude spectrum, while keeping their original phase spectrum. An example of this correction is shown in Fig. S1 of the supplemental material for the EGF1 and the EGF6, which have the smallest and the largest stress drops, respectively.

## EGF variability

The selected EGFs are available at 16 stations located at different azimuths and distances from the fault plane. These stations are indicated by triangles in Fig. 1. Specifically, stations indicated by a blue triangle recorded both the main shock and some of the EGFs, while stations represented by a cyan triangle only recorded some of the EGFs. Nine of these stations belong to the National Accelerometric Network RAN (Rete Accelerometrica Nazionale owned by Civil Protection), four to INGV (Istituto Nazionale di Geofisica e Vulcanologia), one to MEDNET (Mediterranean Network) and two to ISTerre (Institut des Sciences de la Terre de Grenoble).

To investigate the variability of the corrected EGFs at each station, we computed frequency by frequency the standard deviation of the natural logarithm of the velocity spectra. These latter were scaled to get a common reference distance equal to the main shock hypocentral distance, considering a geometrical spreading of $1/r$, consistent with findings of Bindi et al. (2009). Fig. 3 represents the variability in the velocity spectra, averaged in three frequency bands (0.2–1 Hz; 1–3 Hz, 3–20 Hz) as a function of the distance from the main shock hypocentre. The variability does not show a
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Corner frequency as a function of the seismic moment in a logarithmic scale with associated error bars. On the same plot, dashed lines join points at the constant stress drop values of 1 and 10 MPa. Single estimations of the corner frequency are affected by a large uncertainty.</figcaption></figure>

clear dependence on the distance for the analysed frequency bands. Changes are more pronounced in the low and intermediate frequency ranges than in the high frequency interval. In the former case the variability ranges between 0.5 and 1.5. Besides, the com-

parison of the variability at the nearby stations located in L'Aquila city (AQK, AQM, AQU, AQA, AQG and AQV) indicates a difference in the low and intermediate frequency range that does not appear in the high frequency range, where the variability is almost

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Standard deviation of the velocity spectra associated to the selected EGFs in three frequency bands (0.2–1 Hz; 1–3 Hz, 3–20 Hz) for different stations as a function of the hypocentral distance from the main shock. $\sigma_{\ln |\dot{u}(f)|}$ is calculated as the average standard deviation measured frequency by frequency on the natural logarithm of the velocity amplitude spectra. We found that the variability is almost independent of the distance, with larger variations for the low and intermediate frequency ranges.</figcaption></figure>
constant. This variation may be related to changes in the radiation pattern contribution. Since such stations are located just above the fault that generated the main event, the radiation pattern contribution may change as the rupture moves on the fault plane.

DETERMINISTIC SOURCE MODEL
USING KINEMATIC INVERSION

In order to obtain a reference source model to be used in the EGF ground motion simulations, we derived a low-frequency kinematic model for the L'Aquila main shock from inversion of strong motion data. The geometry of the fault was defined according to the focal mechanism given by Chiaraluce et al. (2011), and the dimension of the fault along the strike direction was assumed to be 28 km and the width 20 km (Cirella et al. 2009).

We inverted 13 strong motion data recorded at the RAN network and at station AQU, from MEDNET. The analysed stations are indicated in Fig. 1. Three component acceleration data were twice integrated to obtain displacement, and then band-pass filtered between 0.05 and 0.5 Hz. The lower frequency was selected to ensure a large signal-to-noise ratio on the displacement time histories, while the high-frequency was chosen as the largest frequency for which a 1-D model can be assumed as a good approximation of the structure crossed by the seismic waves (Cirella et al. 2009; Gallovič et al. 2015). Green's functions were computed in the 1-D model of Bianchi et al. (2010) modified by Ameri et al. (2012), using the discrete wavenumber method (AXITRA; Coutant 1994). The inversion was performed by fixing the shape of the slip rate function (we selected a boxcar) and the rise time duration to 1s, this value being compatible with the expected rise time of M 6.3 earthquake (Somerville et al. 1999). It is worth to note that the spectrum of the slip rate function is almost flat up to the inverse of the rise time, independently of the shape of the function. Since this frequency (1 Hz for this case) is significantly above the upper frequency limit of the inverted data (0.5 Hz), the particular choice of the source time function does not significantly affect our results. We should note, however, that the rise time is poorly constrained by the inversion, due to the band limitation of the data. Moreover, we assumed that the slip direction on the fault plane is constant, equal to the rake angle of the focal mechanism.

For the slip parametrization, we used a linear superposition of overlapping Gaussian functions (Lucca et al. 2012). The width of the Gaussian functions was set to 4 km, this value being compatible with the minimum data wavelength. The rupture speed may change in 16 control nodes located on a regular grid, and its value is interpolated elsewhere in all subfaults, by the use of a Lagrangian bi-linear interpolation. Rupture times were then obtained by integration of the rupture velocity, using the eikonal equation (Festa & Zollo 2006). Inversion was again based on the minimization of the $L^2$ norm between observed data and synthetics built on the resulting kinematic model. We used a two-step nested strategy (Lucca et al. 2012), for which the rupture velocity is explored by a global algorithm, while for each fixed rupture velocity model, the slip is retrieved by a non-negative linear inversion, based on the least squares technique.

The kinematic model obtained from the inversion is shown in Fig. 4, while the misfit between real data (black records) and synthetics (red records) is shown in Fig. 5. We found a major slip patch between 5 and 10 km southwards of the hypocentre, which is responsible for the directivity effect observed south of the fault at stations CLN and SUL. We also found a smaller asperity in the up-

per part of the fault with slip as large as 50 cm and a third one nearby the surface, with a 30 cm of slip, both of them being responsible for the up-dip directivity observed at L'Aquila and GSA stations (Fig. 5). Rupture times indicate that the small asperity above the hypocentre was activated before the large asperity. In the upper part of the fault the average rupture velocity is about 2.7 km s⁻¹. The antiplane rupture southward of the hypocentre is delayed in time, with a resulting average rupture velocity of 2.3 km s⁻¹. Location of slip asperities and their time activation are consistent with kinematic models obtained by Cirella et al. (2009) and Gallovič & Zahradník (2011), although absolute amplitudes of slip do not coincide, likely because of differences in model parametrization. Specifically, the slip amplitude in the latter inversion (Gallovič & Zahradník 2011) is larger above the hypocentre than southwards.

The final variance reduction, computed according to Lucca et al. (2012) is 0.60, and for each component of each station are superimposed to the plots of Fig. 5. We observe that in the selected frequency range the synthetics fit pretty well the real observations, at least for two of the three components, and for almost all the stations. Specifically, the amplitudes of one of the horizontal components of the L'Aquila stations AQK and AQU are not reproduced, while the synthetic waveform shape fairly agrees with the real one. Short duration wiggles are not reproduced by the model, indicating that the velocity model is depleting higher-frequency signals, that could be related to the basin effects and resonating monochromatic waves. This effect is added to the simplified description of the source model (smooth parametrization for slip amplitude, constant rake and rise time, bi-linear interpolated rupture velocity) that by itself is not able to catch small scale features on the fault plane.

EFFECT OF THE EGF SELECTION
ON GROUND MOTION PREDICTIONS

Ideally, the set of EGFs $\tilde{G}_n^{(ij)}$ in eq. (4) should cover the whole fault plane, with each $\tilde{G}_n^{(ij)}$ having its hypocentre at the subfault $(i,j)$. In practice, a single EGF is often available and used to perform ground motion simulations. The availability of several good quality EGFs for the L'Aquila earthquake (Table 1) represents then an opportunity to quantify the uncertainty due to the use of various EGFs.

For each station, we performed a set of ground motion simulations, using each of the available EGFs. Hereafter, these sets of simulations will be referred to as 'single EGF' simulations. To better highlight the uncertainty that comes from the EGF choice, we fixed the rupture velocity, the position of the rupture nucleation, the low frequency static slip distribution and the roughness parameter K according to our reference source model obtained by kinematic inversion (Fig. 4). Those simulations are referred to as 'best level of source knowledge' predictions. The roughness parameter K was estimated by fitting the spectrum of the reference static slip model with a $k^{-2}$ model (Somerville et al. 1999). We found a K value of 0.74, which is consistent with a Brune (1970) source model (Causse et al. 2009). The only source of variability remains then the randomness of the static slip phases that controls the location of the small scale slip heterogeneities on the fault plane. This resulted in a set of 25 static slip distributions (Fig. S3). Note that for each subfault $(i,j)$, the selected EGF was scaled by a factor $R_{i,j}/R_{hypo}$, where $R_{i,j}$ is the distance between subfault $(i,j)$ and the receiver, and $R_{hypo}$ is the EGF hypocentral distance, so as to account for geometrical spreading that changes from one subfault to the other. Each EGF was also time shifted by a factor of $(R_{i,j} - R_{hypo})/\beta$, where $\beta$ is the average S-wave velocity along the path between source and
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Final slip distribution and contours of rupture times as obtained from the kinematic inversion of strong motion data. We found a large slip asperity southward the hypocentre, with maximum amplitude of 90 cm and two smaller patches above the hypocentre and close to the surface, that were responsible for the up-dip directivity observed at L'Aquila stations. Orange stars represent the positions of the hypocentres of the EGF with respect to the fault plane.</figcaption></figure>

receiver, to account for differences in source-station distances from one subfault to the other, assuming that the EGF is dominated by *S* waves (e.g. Irikura & Kamae 1994). We next performed a new set of simulations, incorporating the whole set of EGFs. In this approach, referred to as ‘multiple EGFs’, we associated to each subfault (*i*,*j*) the closest EGF. Also in this case, EGFs related to each subfault are scaled to account for the geometrical spreading. The position of the EGFs, superimposed with the kinematic source model, is shown in Fig. 4. It is worth to note that none of the EGFs falls within the large slip patches. Most of them are located nearby the main shock hypocentre while one EGF (E4) is located at the boundary of the largest slip asperity. The number of times that a single EGF is used in the multiple EGF approach is shown in the histograms of Figs 6 and S2 of Supplemental Material. Note that the histograms are station-dependent because not all the EGFs satisfy the quality criteria at some stations.

For each approach ('single EGF' and 'multiple EGFs'), we finally generated 25 time histories by randomly choosing the slip phase spectrum. Fig. 7 displays an example of ground motion simulation at station AQU, located within the projection of the fault plane on the surface, and at station FMG, located about 25 km away from the fault plane, using the same source model. Despite the careful selection of the EGFs and the applied correction, which results in a homogeneous set in terms of moment magnitude, focal mechanism and stress drop, the figure shows that ground motion predictions are strongly sensitive to the EGF choice.

In order to compare the ‘single EGF’ and the ‘multiple EGF’ approaches, we represented the misfit between synthetics and real data for all the stations at which the waveform of the main shock was recorded (Fig. 8). The misfit was computed as the $L^2$ difference on the acceleration Fourier spectra between the average of the 25 simulations and the real data. The comparison indicates that in most cases, the ‘multiple EGF’ approach improves the fit with real data. In Fig. 9 we also compare the average simulated PGA with

the real one (represented by red stars). The PGA was computed as the maximum of the geometrical mean of the two horizontal components of the acceleration. Although for many stations one of the ‘single EGF’ simulations fits observations better than the ‘multiple EGF’ approach, again ‘multiple EGF’ simulations in overall better fit real observations. Figs 8 and 9 show that the use of several EGFs distributed over the fault plane may significantly improve ground motion predictions. A first obvious explanation is that the use of several EGFs may better model the propagation and radiation pattern effects, especially for near field stations. Second, potential azimuthal dependence of site effects may be better accounted for, especially in the context of L’Aquila earthquake where there is a strong coupling between the basin and the fault. Third, the use of several EGFs may reduce epistemic uncertainty associated with the determination of EGF source parameters (moment magnitude and stress drop), and average EGF particular source features, such as stress drop and potential directivity effects.

Since the use of a specific EGF may lead to large differences in the PGA estimation, we investigated the variability associated to the ‘single EGF’ simulations. For each of the 25 source models considered above (slip distributions of Fig. S3 and nucleation, rise time and rupture velocity fixed from the kinematic inversion), we performed ‘single EGF’ simulations at all stations using all the available EGFs per station; we then estimated the PGA and computed the standard deviation of the PGA natural logarithm. Fig. 10 represents the variability averaged over the selected source models (blue dots). This quantity is compared to the variability introduced by the source model (red dots). The latter is computed by estimating the variability in the source models for fixed EGF and averaging the estimations over all the EGFs. In literature this variability is generally referred to as ‘between event variability’, that is the variability one would expect at a given station for repeating earthquakes of the same magnitude. A review of between event variability associated to ground motion prediction equations provided by various authors
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Comparison of real data (solid black curves) and synthetics (dashed red curves) obtained from the kinematic model of Fig. 4. We also report the variance reduction for each component.</figcaption></figure>

is available in Cotton *et al.* (2013), who reported values between 0.2 and 0.5. While the variability in our source models is slightly smaller, the epistemic variability associated to the EGF selection is up to three times larger for the analysed stations (~0.5–1). We do not

observe any distance dependency, although the explored distances are limited to 50 km.

Finally, the PGA observed at the nearby sites AQU and AQK is the same, but simulations provide different results in terms of
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Number of times that a single EGF is used for simulations for the stations which recorded all 10 EGFs (AQU, AQK, CAMP and FAGN), and for the stations GSA, FMG, LG12 and LG13.</figcaption></figure>

PGA. When using the multiple EGF approach, there is a factor of almost 2 in the synthetic estimates between the two stations and the simulations better predict the observation at AQU than at AQK. (Figs 9, 11 and S6). To investigate this discrepancy, we plotted the acceleration records of EGF E4 at both stations (Fig. S4) and the associated amplitude spectra (Fig. S5). This EGF was selected because it is predominantly used in the ground motion simulations at these two stations (Fig. 6). Although the signals are very similar in terms of waveform and duration, the amplitude at AQK is almost twice larger than at AQU. This discrepancy can also be observed in the Fourier spectra in a wide frequency range. Hence we believe that specific soil conditions, eventually characterized by a non-linear behaviour at AQK may be responsible for this difference between EGF and main shock records.

## BLIND GROUND-MOTION PREDICTIONS

A fundamental stage to reduce seismic hazard is the prediction of realistic ground motion for potential future earthquakes. Such predictions, usually referred to as ‘blind’, require the generation of a population of ground motion records that represent the natural variability of the source process for the target earthquake scenario (e.g. Causse *et al.* 2008; Cultrera *et al.* 2010).

We performed blind predictions of the L'Aquila event, assuming that the only available information about the source is the magnitude, the focal mechanism and the overall rupture size, the latter being the same as used in the kinematic inversion. We then defined a priori distributions for the source parameters, based on the analysis of published source models of past earthquakes, and release the source information from the kinematic inversion model used in

the previous section. The distribution of the hypocentre was chosen according to Mai *et al.* (2005), as derived from inspection of a database of kinematic source models obtained from inversion of seismic/geodetic data. We distinguished the position of the nucleation point along the fault strike ($x_n$) and the dip ($y_n$). To avoid a nucleation close to the fault edges, we chose a normal distribution for $x_n$, centred in the middle of the fault plane. For $y_n$, we chose the non-symmetric Weibull distribution to increase the probability of having a nucleation at depth. The distribution of the slip roughness was assumed to be log-normal with a mean value of 0.74, as expected from earthquake self-similarity (Causse *et al.* 2009), and a standard deviation $\sigma_{ln} = 0.28$, estimated by Causse *et al.* (2010). At low wavenumbers, the phase of the static slip spectrum was chosen in order to concentrate the maximum slip in the middle of the fault plane, while at large wavenumbers the phase was randomly chosen (Causse *et al.* 2009). Finally, the rupture velocity (here constant along the fault) was assumed uniformly distributed between 70 per cent and 90 per cent of the shear wave velocity at the source. The rise time is chosen equal to 1s for a $M_w$ 6.3 event, according to Somerville *et al.* (1999). For large wavenumbers, it is conversely proportional to the wavenumber $k$ (see section ‘kinematic source modelling’). Details on the selected distributions and associated parameters are given in Table 2. The Latine hypercube sampling (LHS) method (McKay 1988) was then applied to extract 100 values of each source parameter according to its probability density function. For each station, we then performed 100 simulations combining the selected values of the different source parameters, using the multiple EGFs approach. For comparison with real observations, we computed the median and the standard deviation of both the Fourier and response acceleration spectra. Fig. 11 shows the mean value and the mean value plus/minus one standard deviation of the simulated spectra for the stations AQU and FMG, compared with
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Simulations of the vertical component of the ground motion at the stations AQU and FMG. Real accelerograms (black curve) are compared with several 'single EGF' simulations (red curves) and with the 'multiple EGF' simulation (blue curve). All the records are computed for the same rupture model.</figcaption></figure>

real observations. Fig. S6 shows the same quantities for the remaining stations. The results indicate that the variability of the spectra, measured in natural logarithmic scale, is about 0.5–0.6 in the whole frequency range, but the mean spectrum follows the trend of the real

one almost in the whole frequency range, with several local maxima retrieved. For the horizontal components of near fault stations located at L'Aquila (AQU and AQK), real data falls within about one sigma from the average spectrum in the whole frequency range
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Misfit between synthetics and real data for all the stations at which the main shock was recorded. The misfit is represented as a function of the epicentral distance and was computed for each component (rows) as the $L^2$ norm on the acceleration Fourier spectra between the mean spectrum of the simulations and the spectrum of the real data. We found that in most cases the multiple EGF approach minimizes the misfit.</figcaption></figure>

(0.2–20 Hz). Vertical components, instead, are underestimated at low frequencies (below 1 Hz). Stations on the transect (AQG, AQA and AQV) show a different behaviour: in the low-intermediate frequency range, the simulations overestimate the real data. This anomalous behaviour in the spectrum can be ascribed to a source effect, since velocity data show a narrow peak in this frequency band related to the initial up-dip rupture propagation, eventually coupled with the shallow resonance in the basin (Ameri *et al.* 2012). At

distant stations, real data are quite well reproduced with the exception of station ORC, for which only three EGFs meet the chosen selection criteria.

To investigate the influence of the assumptions underlying the choice of source parameters, in Fig. 12 we compare the ‘blind’ simulations with those calibrated from our low-frequency source inversion model (section ‘Effect of the EGF selection on ground motion predictions’), in which the variability only arises from the

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Comparison between the mean PGA from single EGF (dots) and multiple EGF (triangle) simulations with the PGA measured on real data (star).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Variability of the natural logarithm of the PGA in the single EGF approach as a function of the distance. Blue points correspond to the average value of the variability obtained with source parameters fixed according to our low-frequency deterministic source model and various EGFs. Red points indicate the average variability measured for fixed EGF, and variable phase spectrum of the slip distribution.</figcaption></figure>

slip phase randomness. Fig. 12 also compares the simulated PGA with the real ones. Although we observe a slight improvement in the PGA estimate with simulations obtained fixing the source parameters, we do not observe significant differences between the blind and calibrated simulations, neither in the average value nor in the variability. This is even more evident for the response spectrum at 1 Hz, which is plotted in Fig. 13 as a function of distance. This could be due to the fact that the values of rupture velocity inferred from the kinematic inversion at large wavenumbers are not representative of the evolution of the rupture at smaller scales. Additionally, we assumed a specific constant value of rupture velocity which may introduce some bias in the final estimation of ground motion parameters.

We have also computed the aleatoric variability of the PGA for blind predictions as we did for the variability associated to the source model showed in Fig. 10. Values ranges from ~0.25 to ~0.4 and are of the same order as the ones obtained by Cotton *et al.* (2013).

## DISCUSSION AND CONCLUSIONS

In this study, we performed ground motion simulations coupling a k⁻² slip distribution to represent the high frequency source process and EGFs to model the elastic response of the propagation medium. We analysed as a case study the M6.3, 2009 L'Aquila earthquake, for which a large number of foreshocks and aftershocks was recorded at sites that also recorded the main shock. From a database of more than 400 events, we extracted a subset of only 10 events, which satisfy strict quality criteria based on signal-to-noise ratio, vicinity to the fault plane, similar focal mechanism and stress drop, magnitude

in the narrow range 3.5–3.7. For accurate computation of moment magnitude and stress drop, we inverted the displacement spectra to retrieve the plateau level and the corner frequency of both the EGFs and the main shock. We found that the uncertainty in the corner frequency is as large as the corner frequency value, leading to a large error in the stress drop estimation. We found that the main shock has a corner frequency of 0.25 Hz, yielding a Brune stress drop of 8.9 MPa. Before using EGFs in the simulations, amplitude spectra were corrected to account for the differences in seismic moment and stress drop.

To improve the low-wavenumber description of the source, we inverted strong motion data recorded in the source vicinity to get a low-frequency kinematic description of the rupture. Inversion was performed by a two steps nested procedure, where the slip function was obtained from a non-negative linear inversion for a fixed rupture velocity map. The use of a Gaussian parametrization for the slip allowed to rule out all the wavelengths smaller than the shorter wavelength represented by the data. We found a major slip patch between 5 and 10 km southwards of the hypocentre and two smaller asperities in the upper part of the fault. The up-dip patches were activated before the main asperity and because isochrones for near-fault sites are elongated in the up-dip direction, the signal from the main asperity was delayed in time at L'Aquila stations. Location, moment release and time activation of the asperities are compatible with other models in literature (Cirella *et al.* 2009; Gallovič & Zahradník 2011) and were related to fluid filled regions, characterized by large values of the ration between P- and S-wave velocities in tomographic models (Di Stefano *et al.* 2011).

<table><caption>**Table 2.** Source parameters, associated probability density functions used for ground motion simulations, and selected parameters for the distributions.</caption><tr><td>Rupture velocity v<sub>r</sub></td><td>Uniform distribution</td><td>Lower value: 2.1 km s<sup>-1</sup></td><td>Upper value: 2.7 km s<sup>-1</sup></td></tr><tr><td>Nucleation point along strike, x<sub>n</sub> [0;1]</td><td>Normal distribution</td><td>Mean value: 0.5</td><td>Standard deviation: 0.23</td></tr><tr><td>Nucleation point along dip, y<sub>n</sub> [0;1]</td><td>Weibull distribution</td><td>a = 0.63</td><td>b = 3.92</td></tr><tr><td>Roughness parameter, K</td><td>Lognormal distribution</td><td>Mean value: K = 0.74</td><td>Standard deviation: σ<sub>ln</sub> = 0.28</td></tr></table>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Comparison between blind simulations (grey curves) and real observations (red curves) for the stations AQU and FMG in terms of Fourier acceleration spectra (left-hand panels) and response spectra (right-hand panels). Black solid lines represent the mean spectra while black dashed lines represent the mean spectrum ±1 standard deviation.</figcaption></figure>

Since we were able to select a database of several EGFs, we investigated the variability in the ground motion simulations related to the use of a single EGF and compared these simulations to those obtained using the whole set of EGFs. For this comparison, ground motion simulations were performed using source parameters calibrated from our low-frequency kinematic inversion model. Variability in the source process was then introduced by choosing random phases of the slip spectrum at high wavenumbers. We compared the results of the simulations to the real observations in terms of the $L^2$ difference of the acceleration spectra. We found that the use of multiple EGFs improves the fit with data as compared to the single EGF approach. This is likely due to a better representation of path effects and changes in the focal mechanism. In addition,

position of the station may change on the focal sphere for the different EGFs, at least for near source stations. Finally, we found that the epistemic variability related to the EGF choice in ‘single EGF’ predictions (~0.5–1) is up to three times larger than the between event variability empirically retrieved from ground motion prediction equations by Cotton *et al.* (2013) (~0.4).

We finally released all the hypothesis on the kinematic model (low-frequency slip, rupture velocity and slip roughness parameter) and performed blind simulations for the L’Aquila event, assuming a *priori* distributions of source parameters (Table 2). We found that both Fourier and response acceleration spectra reproduce the main data trend. In the intermediate frequency range (0.5–1 Hz), real observations are about one sigma distant from the median simulations,
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Comparison between the average PGA from blind simulations (black) and simulations obtained with source parameters fixed according to our low frequency deterministic source model (blue). Observed PGA are represented by red triangles.</figcaption></figure>

with equal occurrence of over and underestimations. This might be due to a coherent slip patch above the hypocentre for the real event that generated an impulsive signal at L'Aquila stations (Ameri et al. 2012). We also compared 'blind' simulations with the ones calibrated from our low-frequency kinematic inversion model in terms of PGA and response spectral ordinate at 1 s. Except for a few cases, we observe similar mean values for PGA and $S_a(1s)$. This is not surprising because the deterministic source model mainly affects the low frequency ground motion (up to ~0.5 Hz),

while $S_a(1s)$, and particularly PGA, are controlled by higher frequencies. Hence calibrated predictions do not provide any significant improvement as compared to the blind ones. Moreover for some stations blind simulations approach real observations better than simulations with low frequency calibrated source model. There are two main possibilities to explain why constraining source parameters do not always lead to a better prediction. First, the extrapolation of the low-frequency description of the rupture, as inferred by the kinematic inversion, may introduce some bias in the final

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Spectral ordinate at 1 s for blind simulations (black) and for fixed low frequency source model (blue). Red triangles represent the observed values.</figcaption></figure>
simulations. Second, the low frequency source model has been obtained using Green's functions computed in a 1-D-layered velocity model, while our predictions use EGFs as an approximation of the medium response.

Finally we have computed the aleatoric variability of the PGA associated to blind predictions. Values ranges from ~0.25 to ~0.4 and are of the same order as the ones obtained by Cotton *et al.* (2013). They are also significantly smaller than the epistemic variability associated to the choice of a specific EGF. This is probably due to the averaging effect of the multiple EGFs configuration.

## ACKNOWLEDGEMENTS

We are grateful to Martin Mai and an anonymous reviewer, who helped in improving the manuscript. This research activity was supported by the EU project NERA, Network of European Research Infrastructures for Earthquake Risk Assessment and Mitigation. NERA has received funding from the European Union's Seventh Programme for research, Combination of Collaborative Project and Coordination and Support action, under grant agreement No. 262330.

## REFERENCES

* Aki, K. & Richards, P.G., 2002. *Quantitative Seismology*, 2nd edn, University Science Books.

* Ameri, G., Bindi, D., Pacor, F. & Galadini, F., 2011a. The 2009 April 6, $M_w$ 6.3, L'Aquila (central Italy) earthquake: finite-fault effects on intensity data, *Geophys. J. Int.*, **186**, 837–851.

* Ameri, G., Emolo, A., Pacor, F. & Gallovič, F., 2011b. Ground-motion simulations for the M 6.9 Irpinia 1980 Earthquake (Southern Italy) and scenario events, *Bull. seism. Soc. Am.*, **101**, 1136–1151.

* Ameri, G., Gallovič, F. & Pacor, F., 2012. Complexity of the $M_w$ 6.3 2009 L'Aquila (central Italy) earthquake: 2. Broadband strong motion modeling, *J. geophys. Res.*, **117**, B04308, doi:10.1029/2011JB008729.

* Ben-Menahem, A., 1961. Radiation of seismic surface-waves from finite moving sources, *Bull. seism. Soc. Am.*, **51**, 401–435.

* Bernard, P., Herrero, A. & Berge, C., 1996. Modeling directivity of heterogeneous earthquake ruptures, *Bull. seism. Soc. Am.*, **86**, 1149–1160.

* Bianchi, I., Chiarabba, C. & Piana Agostinetti, N., 2010. Control of the 2009 L'Aquila earthquake, central Italy, by a high-velocity structure: a receiver function study, *J. geophys. Res.*, **115**, B12326, doi:10.1029/2009JB007087.

* Bindi, D., Pacor, F., Luzi, L., Massa, M. & Ameri, G., 2009. The $M_w$ 6.3, 2009 L'Aquila earthquake: source, path and site effects from spectral analysis of strong motion data, *Geophys. J. Int.*, **179**, 1573–1579.

* Bour, M. & Cara, M., 1997. Test of a simple empirical Green's function method on moderate-sized earthquakes, *Bull. seism. Soc. Am.*, **87**, 668–683.

* Brune, J.N., 1970. Tectonic stress and the spectra of seismic shear waves from earthquakes, *J. geophys. Res.*, **75**, 4997–5009.

* Causse, M., Cotton, F., Cornou, C. & Bard, P.Y., 2008. Calibrating median and uncertainty estimates for a practical use of empirical Green's functions technique, *Bull. seism. Soc. Am.*, **98**, 344–353.

* Causse, M., Chaljub, E., Cotton, F., Cornou, C. & Bard, P.Y., 2009. New approach for coupling $k^{-2}$ and empirical Green's functions: application to the blind prediction of broadband ground-motion in the Grenoble basin, *Geophys. J. Int.*, **179**(3), 1627–1644.

* Causse, M., Cotton, F. & Mai, P.M., 2010. Constraining the roughness degree of slip heterogeneity, *J. geophys. Res.*, **115**, B05304, doi:10.1029/2009JB006747.

* Chiarabba, C. *et al.*, 2009. The 2009 L'Aquila (central Italy) MW6.3 earthquake: main shock and aftershocks, *Geophys. Res. Lett.*, **36**, L18308–L18314.

* Chiaraluce, L., Valoroso, L., Piccinini, D., Di Stefano, R. & De Gori, P., 2011. The anatomy of the 2009 L'Aquila normal fault system (central Italy) imaged by high resolution foreshock and aftershock locations, *J. geophys. Res.*, **116**, B12311, doi:10.1029/2011JB008352.

* Chiou, B., Darragh, R., Gregor, N. & Silva, W., 2008. NGA project strong-motion database, *Earthq. Spectra*, **24**(1), 23–44.

* Cirella, A., Piatanesi, A., Cocco, M., Tinti, E., Scognamiglio, L., Michelini, A., Lomax, A. & Boschi, E., 2009. Rupture history of the 2009 L'Aquila (Italy) earthquake from non-linear joint inversion of strong motion and GPS data, *Geophys. Res. Lett.*, **36**, L19304–L19309.

* Cirella, A., Piatanesi, A., Tinti, E., Chini, M. & Cocco, M., 2012. Complexity of the rupture process during the 2009 L'Aquila, Italy, earthquake, *Geophys. J. Int.*, **190**, 607–621.

* Cotton, F., Archuleta, R. & Causse, M., 2013. What is sigma of the stress drop?, *Seismol. Res. Lett.*, **84**, 42–48.

* Coutant, O., 1994. Expression of the Green's functions in cylindrical coordinates used with a reeffictivity method. Axitra program, ftp://ftp-lgit obs.ujf-grenoble.fr/pub/archive/axitra/axitra.ps

* Cultrera, G., Cirella, A., Spagnuolo, E., Herrero, A. & Pacor, F., 2010. Variability of kinematic parameters and its implication on the choice of the design scenario, *Bull. seism. Soc. Am.*, **100**, 941–953.

* Cuppillard, P., Delavaus, E., Burgos, G., Festa, G., Vilotte, J.P., Capdeville, Y. & Montagner, J.P., 2012. RegSEM: a versatile code based on the spectral element method to compute seismic wave propagation at the regional scale, *Geophys. J. Int.*, **188**(3), 1203–1220.

* Das, S. & Henry, C., 2003. Spatial relation between main earthquake slip and its aftershock distribution, *Rev. Geophys.*, **41**(3), doi:10.1029/2002RG000119.

* Di Stefano, R., Chiarabba, C., Chiaraluce, L., Cocco, M., De Gori, P., Piccinini, D. & Valoroso, L., 2011. Fault zone properties affecting the rupture evolution of the 2009 (Mw 6.1) L'Aquila earthquake (central Italy): insights from seismic tomography, *Geophys. Res. Lett.*, **38**, L10310, doi:10.1029/2011GL047365.

* Donahue, J. & Abrahamson, N., 2014. Simulation-based hanging-wall effects, *Earthq. Spectra*, **30**(3), 1269–1284in-press.

* Festa, G. & Zollo, A., 2006. Fault slip and rupture velocity inversion by isochrone back-projection, *Geophys. J. Int.*, **166**, 745–756.

* Festa, G. & Zollo, A., 2012. From data to source parameters: kinematic modeling, in *The Mechanics of Faulting: From Laboratory to Real Earthquakes*, Vol. 37/66, Issue 2, eds Bizzarri, A. & Bhat, H.S., Research Signpost.

* Gallovič, F. & Brokešová, J., 2004. On strong ground motion synthesis with the $k^{-2}$ slip distribution, *J. Seism.*, **8**, 211–224.

* Gallovič, F. & Brokešová, J., 2007. Hybrid k-squared source model for strong ground motion simulations: introduction, *Phys. Earth planet. Inter.*, **160**(1), 34–50.

* Gallovič, F. & Brokešová, J., 2008. Probabilistic aftershock hazard assessment II: Application of strong ground motion modeling, *J. Seism.*, **12**(1), 65–78.

* Gallovič, F. & Zahradník, J., 2011. Complexity of the M6.3 2009 L'Aquila (central Italy) earthquake: 1. Multiple finite-extent source inversion, *J. geophys. Res.*, **117**, B04307, doi:10.1029/2011JB008709.

* Gallovič, F., Imperatori, W. & Mai, P.M., 2015. Effects of three-dimensional crustal structure and smoothing constraint on earthquake slip inversions: case study of the $M_w$ 6.3 2009 L'Aquila earthquake, *J. geophys. Res.*, **120**, 428–449.

* Graves, R.W., 1996. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences, *Bull. seism. Soc. Am.*, **86**, 1091–1106.

* Graves, R.W., Aagaard, B.T., Hudnut, K.W., Star, L.M., Stewart, J.P. & Jordan, T.H., 2008. Broadband simulations for $M_w$ 7.8 southern San Andreas earthquakes: ground motion sensitivity to rupture speed, *Geophys. Res. Lett.*, **35**, L22302, doi:10.1029/2008GL035750.

* Hall, J.F., Heaton, T.H., Halling, M.W. & Wald, D.J., 1995. Near-source ground motion and its effects on flexible buildings, *Earthq. Spectra*, **11**, 569–605.

* Hartzell, S.H., 1978. Earthquakes aftershocks as Green's functions, *Geo-phys. Res. Lett.*, **5**, 1–4.
Hartzell, S., Harmsen, S., Frankel, A. & Larsen, S., 1999. Calculation of broadband time histories of ground motion: comparison of methods and validation using strong-ground motion from the 1994 Northridge earthquake, *Bull. seism. Soc. Am.*, **89**, 1484–1504.

Herrero, A. & Bernard, P., 1994. A kinematic self-similar rupture process for earthquakes, *Bull. seism. Soc. Am.*, **84**, 1216–1228.

Herrero, A. & Bernard, P., 1996. Modeling directivity of heterogeneous earthquake ruptures, *Bull. seism. Soc. Am.*, **86**(4), 1149–1160.

Heuze, F., Archuleta, R., Bonilla, F., Day, S., Doroudian, M., Elgamal, A. & Yang, Z., 2004. Estimating site-specific strong earthquake motions, *Soil Dyn. Earthq. Eng.*, **24**, 199–223.

Hutchings, L., 1994. Kinematic earthquake models and synthesized ground motion using empirical Green's functions, *Bull. seism. Soc. Am.*, **84**, 1028–1050.

Hutchings, L. & Wu, F., 1990. Empirical Green's functions from small earthquakes: a waveform study of locally recorded aftershocks of the 1971 San Fernando earthquake, *J. geophys. Res.*, **95**, 1187–1214.

Hutchings, L., Ioannidou, E., Foxall, W., Voulgaris, N., Savy, J., Kalogeras, I., Scognamiglio, L. & Stavrakakis, G., 2007. A physically based strong ground-motion prediction methodology: application to PSHA and the 1999 $M_w = 6.0$ Athens earthquake, *Geophys. J. Int.*, **168**, 659–680.

Irikura, K. & Kamae, K., 1994. Estimation of strong ground motion in broad-frequency band based on a seismic source scaling model and an empirical Green's function technique, *Ann. Geophys.*, **37**(6), 1721–1743.

Ishii, M., Shearer, P.M., Houston, H. & Vidale, J.E., 2005. Extent, duration and speed of the 2004 Sumatra–Andaman earthquake imaged by the Hi-Net array, *Nature*, **435**, 933–936.

Kohrs-Sansorny, K., Courboulex, F., Bour, M. & Deschamps, A., 2005. A two-stage method for strong ground motion simulations using stochastic summation of small earthquakes, *Bull. seism. Soc. Am.*, **94**, 314–319.

Koketsu, K. & Miyake, H., 2008. A seismological overview of long period ground motion, *J. Seismol.*, **12**, 133–143.

Komatitsch, D. & Vilotte, J.P., 1998. The spectral element method; an efficient tool to simulate the seismic response of 2D and 3D geological structures, *Bull. seism. Soc. Am.*, **88**, 368–392.

Lucca, E., Festa, G. & Emolo, A., 2012. Kinematic inversion of strong motion data using a Gaussian parameterization for the slip: application to the 2008 Iwate-Miyagi, Japan, earthquake, *Bull. seism. Soc. Am.*, **102**(6), 2685–2703.

Maercklin, N., Festa, G., Colombelli, S. & Zollo, A., 2012. Twin ruptures grew to build up the giant 2011 Tohoku, Japan, earthquake, *Scientific Reports*, Vol. 2, 709 pp.

Mai, P.M. & Thingbaijam, K.K.S., 2014. SRCMOD: an online database of finite-fault rupture models, *Seism. Res. Lett.*, **85**, 1348–1357.

Mai, P.M., Spudich, P. & Boatwright, J., 2005. Hypocenter locations in finite-source rupture models, *Bull. seism. Soc. Am.*, **95**, 965–980.

Mai, P.M., Burjanek, J., Delouis, B., Festa, G., Francois-Holden, C., Monelli, D. & Zahradnik, J., 2007. Earthquake source inversion blindtest: initial results and further developments, in *AGU Fall Meeting Abstracts*, Vol. 1, 08 pp.

McKay, M.D., 1988. Sensitivity and uncertainty analysis using a statistical sample of input values, in *Uncertainty Analysis*, pp. 145–186, ed. Ronen, Y., CRC Press.

Miyake, H. & Koketsu, K., 2005. Long-period ground motions from a large offshore earthquake: the case of the 2004 off the Kii peninsula earthquake, *Japan, Earth Planet Space*, **57**, 203–207.

Miyake, H., Iwata, T. & Irikura, K., 2003. Source characterization for broadband ground-motion simulation: kinematic heterogeneous source model and strong motion generation area, *Bull. seism. Soc. Am.*, **93**, 2531–2545.

Nozu, A. & Irikura, K., 2008. Strong-motion generation areas of a great subduction-zone earthquake: waveform inversion with empirical Green's functions for the 2003 Tokachi-oki earthquake, *Bull. seism. Soc. Am.*, **98**, 180–197.

Oglesby, D.D. & Mai, P.M., 2012. Fault geometry, rupture dynamics and ground motion from potential earthquakes on the North Anatolian Fault under the Sea of Marmara, *Geophys. J. Int.*, **188**(3), 1071–1087.

Olsen, K.B., 2000. Site amplification in the Los Angeles basin from three-dimensional modeling of ground motion, *Bull. seism. Soc. Am.*, **90**, S77–S94.

Pavic, R., Koller, M.G., Bard, P.Y. & Lacave-Lachet, C., 2000. Ground motion prediction with the empirical Green's functions technique: an assessment of uncertainties and confidence level, *J. Seismol.*, **4**, 59–77.

Scognamiglio, L., Tinti, E., Michelini, A., Dreger, D., Cirella, A., Cocco, M., Mazza, S. & Piatanesi, A., 2010. Fast determination of moment tensors and rupture history: application to the April 6th 2009, L'Aquila earthquake, *Seismol. Res. Lett.*, **81**(6), 892–906.

Shao, G.F. & Ji, C., 2012. What the exercise of the SPICE source inversion validation BlindTest I did not tell you, *Geophys. J. Int.*, **189**, 569–590.

Somerville, P., 2003. Magnitude scaling of the near fault rupture directivity pulse, *Phys. Earth planet. Inter.*, **137**, 201–212.

Somerville, P, et al., 1999. Characterizing crustal earthquake slip models for the prediction of strong ground motion, *Seismol. Res. Lett.*, **70**, 59–80.

Stupazzini, M., Paolucci, R. & Igel, H., 2009. Near-fault earthquake ground-motion simulation in the Grenoble Valley by a high-performance spectral element code, *Bull. seism. Soc. Am.*, **99**, 286–301.

Xu, J., Bielak, J., Ghattas, O. & Wang, J., 2003. Three-dimensional nonlinear seismic ground motion modeling in basins, *Phys. Earth planet. Inter.*, **137**, 81–95.

Zollo, A., Orefice, A. & Convertito, V., 2014. Source parameter scaling and radiation efficiency of microearthquakes along the Irpinia fault zone in southern Apennines, Italy, *J. geophys. Res.: Solid Earth*, **119**, 3256–3275.

## SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this paper:

**Figure S1.** Fourier amplitude displacement spectra of the east-west component of EGF1 and EGF6 recorded at the AQK station. Left-hand panel shows the original spectra and right-hand panel shows the same spectra after the correction procedure, assuming a Brune's source model.

**Figure S2.** Number of times that a given EGF is used for ground-motion simulations at stations AQA, AQG, AQM, AQV, MTR, ORC, RM04 and RM13.

**Figure S3.** Set of 25 k-square slip distributions used in section 'Effect of the EGF selection on ground motion predictions'. Distributions are obtained by fixing the low-frequency slip according to our kinematic inversion, and by randomly choosing high frequency slip phases.

**Figure S4.** Acceleration records of the EGF4 at station AQU (left-hand panel) and AQK (right-hand panel). Amplitudes at AQK are around twice larger as at AQU.

**Figure S5.** Comparison between Fourier amplitude acceleration spectra of EGF4 recorded at stations AQU (black lines) and AQK (red lines).

**Figure S6.** Comparison between blind simulations (grey curves) and real observations, when available (red curves) for stations AQA, AQG, AQK, AQM, AQV, CAMP, FAGN, GSA, LG12, LG13, MTR, ORC, RM04 and RM13 in terms of Fourier acceleration spectra (left-hand panels) and response spectra (right-hand panels). Black solid lines represent the mean spectra while black dashed lines represent the mean spectrum $\pm 1$ standard deviation (http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggv325/-/DC1).

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.