Computation of amplification factor of earthquake ground motion for a local sedimentary structure

We present methodology of calculating acceleration and corresponding earthquake ground motion characteristics at a site of interest assuming acceleration at a reference site for two basic configurations. In one configuration we assume that the reference ground motion is not affected by the local structure beneath the site of interest. In the other configuration we assume that the reference ground motion is affected by the local structure. Consequently, the two configurations differ from each other by the presence of the reference site within the computational model. For each of the two configurations we assume two wavefield excitations: a vertical plane-wave incidence and a point double-couple source. We illustrate the methodology on the example of the Grenoble valley. The extensive investigation of effects of local surface sedimentary structures based on the developed methodology is presented in the accompanying article by Moczo et al. (Bull Earthq Eng, 2018) in this volume.


3 Introduction
SIGMA was a Research & Development program of EDF (Electricité de France), AREVA, CEA (Commissariat à l'énergie atomique et aux énergies alternative) and ENEL in 2011(Senfaute et al. 2015Pecker et al. 2016Pecker et al. , 2017. Being focused on characterization of seismic ground motion in France and nearby countries, its main goal was to obtain robust and stable estimates of the seismic hazard based on proper characterization of uncertainties. The Work Package 3 (WP3) of the project aimed in developing methods of predicting whether a specific site needs a special investigation with respect to its site conditions in view of including site effects in the seismic hazard estimation. One of specific items of WP3 was to investigate potential of a few typical, virtual site configurations to undergo specific site effects associated to the 2D or 3D geometry of the underground structure. This was done using 1D, 2D and 3D forward numerical simulations to identify key parameters controlling the site amplification, and to quantitatively characterize their amount, with a special focus on the modifications to "classical" 1D effects by the 2D or 3D geometry.
A set of 7 nominal models for typical, virtual underground configurations has been defined (including a few real sites such as the Grenoble Alpine valley), together with a set of modifications of the nominal models to investigate effects of variations of a few, carefully selected, structural parameters on earthquake ground motion. Forward numerical simulations were performed in the linear domain with the finite-difference (FD) method Chaljub et al. 2010Chaljub et al. , 2015. In addition to 3D simulations for 3D models assuming a vertical plane wave incidence and/or point double-couple (DC) source, 2D simulations were performed for selected 2D profiles in the 3D models and several 2D nominal models assuming the vertical plane wave incidence, together with 1D simulations for 1D models for selected receiver positions along 2D profiles and in 2D models.
If no records or insufficient number of records are available for the investigated site it is necessary to account for a potential variability of earthquake ground motion. It is reasonable to use a set of properly selected accelerograms recorded at different locations to represent the ground motion variability. We may assume that the records represent ground motions at a free surface of a halfspace. This assumption corresponds to the typical situation in probabilistic seismic hazard assessment (PSHA): the estimated characteristics of the earthquake ground motion relate to outcropping "standard" rock. The question is then how the presence of a local surface sedimentary structure (LSSS) modifies the ground motion. In a computational (numerical-modelling) approach we can quantify an effect of LSSS by considering two separate computational models-a model of LSSS (spatial domain for numerical modelling) that does not include the reference site, and the model of the reference site.
In this article we present methodology for computational estimate of acceleration and other ground motion characteristics at a site of interest for a set of accelerograms corresponding to the ground motion at the reference site. We derive formulas for the site acceleration if the reference site is not in the model. First we assume that the ground motion is due to a plane-wave excitation. This assumption is very common but not very realistic. We thus also assume that the ground motion is due to a point DC source. The point-source assumption has also limitation because it is not valid for larger magnitudes. The point-source assumption is inadequate if there is a reasonable argument to assume an extended source near the site of interest. Then it is more appropriate to perform a case study. We complete the methodology by defining the frequency-dependent and single-valued amplification factors for five selected earthquake 1 3 ground-motion characteristics. We illustrate the methodology by a numerical example for the typical Alpine Grenoble valley in France. Finally we present (in Appendices 1 and 2) derivations of formulas for site accelerations in situation when we cannot exclude the effect of LSSS on ground motion at a reference site and both the reference site and LSSS have to be included in one computational model.
The application of the methodology is presented in the accompanying article by Moczo et al. (2018).

Computation of amplification factors of earthquake ground motion
A very common assumption is that the reference motion is not affected by LSSS. It is not correct as there always exists some diffraction from local heterogeneities, and as "reference" motion is derived from instrumental recordings at real sites where such diffracted waves certainly exist. However, the associated artefacts can reasonably be considered negligible. Such a situation can thus be computationally implemented using a model that includes LSSS but does not include the site at which the reference accelerogram was recorded. The configurations are shown in Fig. 1 for a ground motion due to a vertically incident plane wave and in Fig. 2 for a ground motion due to a DC point source.

Reference site is not in the model, plane-wave excitation
For this configuration, shown in Fig. 1, we assume ground motion due to a vertical incidence of a plane wave. Fig. 1 Illustration of input and output signals. ⃗ a i (t);i ∈ {1, … , n} is the reference acceleration at reference site HAL at the free surface of a homogeneous halfspace. Index i denotes the i-th of a set of n selected reference accelerograms. The corresponding acceleration at site of interest SIT is ⃗ s i (t) . If p ξ is a pseudoimpulse input signal (in particle velocity) of an incident plane wave polarized in the ξ-direction, r ξη is the η-component of the pseudoimpulse response (in the particle velocity) at site SIT

Fig. 2
Illustration of input and output signals. ⃗ a i (t);i ∈ {1, … , n} is the reference acceleration at reference site HAL at the free surface of a homogeneous halfspace. The corresponding acceleration at site of interest SIT is ⃗ s i (t) . If an elementary source e causes particle velocity ⃗ s e,HAL at HAL, ⃗ s e,SIT is the particle velocity at SIT 1 3 Pseudoimpulse input signal As all site responses are performed in the linear domain, for computation of the pseudoimpulse responses it is possible to use a Gabor signal as time function of the incident plane wave in the particle velocity: Here ω p = 2πf p , γ s controls the width of the signal, θ is a phase shift. For obtaining the transfer properties at a site it is reasonable to assume Denote the Fourier spectrum of the input signal by  p(f ).
Matrix of the time-domain pseudoimpulse responses at a site A plane wave polarized in the ξ-direction, ξ ∊ {x, y, z}, causes pseudoimpulse responses (in the particle velocity) r ξx (t), r ξy (t) and r ξz (t) at site SIT. The second index indicates the component. The matrix of the time-domain pseudoimpulse responses is then and its Fourier transform is  . Paolucci (1999), the matrix defined as characterizes transfer properties of the model between the horizontal plane at which the excitation by the plane wave is applied and site SIT. Paolucci (1999) has already given examples that the off-diagonal terms of this matrix may be significant, and that the 2D or 3D effects also induce a significant scattering on the diagonal terms.

Matrix of the Fourier transfer functions As proposed in
Acceleration at the free surface of a halfspace Assume acceleration ⃗ a i (t) at site HAL, that is, at the free surface of a homogeneous halfspace. This means, that 1 2 ⃗ a i (t) is the acceleration of the vertically incident plane wave.
Note, however, that in the numerical simulations we cannot use exactly 1 2 ⃗ a i (t) for a convolution in the local structure. This is because the numerically evaluated transfer function includes effects of grid dispersion (although optionally small, corresponding to a chosen spatial discretization). Consequently, if we replaced the local structure by a homogeneous medium (getting thus the model of a homogeneous halfspace), we would not get exactly ⃗ a i (t) at the free surface for 1 2 ⃗ a i (t) in the incident wave. Therefore we simulate propagation of the plane wave in the grid and obtain 1 2 ⃗ a i (t) affected by the grid dispersion. We then apply such numerically obtained 1 2 ⃗ a i (t) in the convolution.
Site acceleration If ⃗ a i (t) is the acceleration at site HAL, then the corresponding acceleration at site SIT is ⃗ s i (t): (1) r xx r yx r zx r xy r yy r zy r xz r yz r zz

Reference site is not in the model, excitation by a double-couple point source
In the configuration shown in Fig. 2 we assume that a reference acceleration ⃗ a i (t);i ∈ {1, … , n}, at the free surface of a homogeneous halfspace is caused by some DC point source at a specified hypocentre position. We want to know the corresponding acceleration ⃗ s i (t) at a site of interest SIT. First we have to find the DC mechanism. Because an arbitrary DC point source can be expressed as a linear combination of 6 independent elementary sources (canonically oriented dipoles and couples) we have to numerically simulate surface motion separately for each of the 6 elementary sources at the hypocentre position. From the 3 components of the reference acceleration and 6 × 3 components of the accelerations caused by the 6 elementary sources we have to find 6 coefficients of the linear combination of the elementary sources that gives the sought DC mechanism producing ⃗ s i (t) at SIT. Let us note why we have to determine specific mechanism if we know the source of each of the accelerograms considered for representing scatter in the wavefield excitation. We want to assume, e.g., 3 hypocentre positions for the Grenoble valley (the case investigated in the companion paper). We want to assume for each chosen point-source hypocentre all reference accelerograms. Therefore we have to find specific mechanism for each of the accelerograms separately. Such a procedure also allows to compare the amplification factors for point sources and incident plane waves.
If we simulate surface motions for each of 6 elementary sources also in the model with LSSS, we can use coefficients of linear combination for ⃗ a i (t) for evaluating ⃗ s i (t) corresponding to the ith DC source.
In the following we present the detailed theory.
Elementary sources e 6 elementary sources comprise 3 dipoles and 3 double couples.

Their moment tensors are
Particle velocity at halfspace due to an elementary source An elementary source e, with the assumed normalized source-time function, acting in the homogeneous halfspace causes particle velocity ⃗ s e, HAL ;e ∈ {1, 2, … , 6} at HAL at the free surface of a homogeneous halfspace. The matrix of the elementary particle-velocity seismograms at HAL for all 6 elementary sources is then defined as Let HAL elem denote the matrix of the corresponding acceleration seismograms.
Particle velocity at a site due to an elementary source An elementary source e acting in the model with LSSS causes particle velocity ⃗ s e, SIT ;e ∈ {1, 2, … , 6} at SIT. The matrix of the elementary particle-velocity seismograms at SIT for all 6 elementary sources is then defined as Let SIT elem Denote the matrix of the corresponding acceleration seismograms.

Acceleration at the free surface of a halfspace Consider acceleration
At each frequency f we want to find coefficients c 1,i , c 2,i , c 3,i , c 4,i , c 5,i , c 6,i of the linear combination of the elementary solutions such that that is, Matrix  HAL elem is the 3 × 6 matrix with complex elements. Because an inverse matrix to the 3 × 6 matrix does not exist, we only can find a pseudo-inverse matrix and consequently determine coefficients c 1,i , c 2,i ,…, c 6,i . Matrix  HAL elem can be decomposed using the singular value decomposition (SVD) method: Here U is the unitary 3 × 3 matrix and V is the unitary 6 × 6 matrix. † is the Hermitian conjugate (or adjoint) matrix to V: Here * indicates a complex conjugate matrix. The unitary matrices satisfy relations where I stands for the 3 × 3 and 6 × 6 identity matrices, respectively. S is the 3 × 6 matrix with elements S ij = 0 if i ≠ j. The three diagonal elements may be denoted by S j for j = 1, 2, 3 and are termed the singular values of matrix  HAL elem : Substituting  HAL elem in Eq. (10) by the r.h.s. of Eq. (11), and then by sequential multiplying Eq. (10) by † , S −1 and V we obtain where Site acceleration If ⃗ a i (t) is the acceleration at HAL, then the corresponding acceleration at SIT is ⃗ s i (t) . Matrix MESH represents relation between the ground motions at HAL and SIT assuming that the wavefield and ground motion were generated by a point DC source. The matrix has the meaning of the spectral matrix ratio  SIT elem  HAL elem −1 . Therefore it is equivalent to the analogous spectral matrix ratio  SIT elem  HAL elem −1 for the particle velocities. Consequently, where subscript S indicates that the three matrices relate to decomposition In fact, given our numerical-modelling method, the velocity-stress FD scheme, we primarily obtain  HAL elem and  SIT elem , and therefore we make use of relation (21) for determining MESH.

Amplification factors
The amplification factors are looked for on several ground motion characteristics (single-valued or frequency-dependent vector values as specified in Sect. 2.2.2) that are not related linearly with their values for the input motion (unlike a Fourier spectral ratio). It is thus needed to consider several realistic input accelerograms, in order to get robust estimates on the corresponding average amplification factors (and their signal-to-signal variability).

Selection of reference accelerograms
It has been shown in previous studies (e.g., Biro and Renault 2012;Bora et al. 2015Bora et al. , 2016 that the amplification factors of response spectral ordinates are sensitive to the frequency contents of the input motion. Consequently, it is reasonable to select the input accelerograms on the basis of their frequency contents. The selection should be based on the following criteria: • motion recorded on rock or stiff soil sites, • motion recorded in the near source area, • magnitude range, • very good signal-to-noise ratio over a wide frequency band with sufficiently low high-pass frequency, • wide distribution of peak frequencies of the peak acceleration response spectrum within the considered frequency range, • sufficient number of records for a meaningful statistical analysis.
Specific choice and a table of selected accelerograms will be given in Sect. 3 presenting a numerical example.

Frequency-dependent amplification factors
Let s ξ,i be the ξth component of a site acceleration corresponding to the reference ith accelerogram ⃗ a i . Then we can define the following frequency-dependent and single-valued amplification factors between the site and reference accelerations.
Amplification factor It is defined as the ratio of the relative displacement response spectra S D of the site acceleration, s ξ,i (t), and acceleration taken as a reference, a ξ,i (t): The amplification factor for the horizontal component may be defined as S D a x,i (f ;5%) S D a y,i (f ;5%)

3
Average amplification factor It is defined for a set of n accelerograms for the ξ-component as where AF ξ,i (f) is the amplification factor for the ith accelerogram and ξ ∊ {x, y, z, h}.
Standard deviation of the amplification factor It is defined as

Single-valued amplification factors
In addition to the frequency-dependent response spectra, a few single-valued ground motion intensity measures are often used to characterize the level of ground motion: amplification factors can therefore be defined also for such parameters. Here we consider some of them, such as the average amplifications factors in some absolute or site-related frequency ranges. (1994) to characterize the short and mid period amplification factors for building codes. Considering the fact that such factors were never given an exact, unambiguous mathematical definition (see for instance the variability of their definitions in the summary performed by Power et al. (2004) within the NGA framework, we adopted here the following definitions, to ensure full reproducibility and consistency in terms of frequency bandwidth. For the ξ-component at a site they are derived from the average, frequency-dependent amplification factors defined in Eq. (25) using a geometrical average over a 2-octave bandwidth centred on 0.1 s (10 Hz) for short period and 1 s (1 Hz) for mid-period.

and
The central short period (0.1 s) was deliberately considered lower than in Borcherdt (1994) (0.3 s), considering the importance of high-frequency equipment in the nuclear industry, and the fact that in the response spectra domain, as shown by Bora et al. (2016), the "high-frequency" values are controlled not only by the Fourier content at the same frequency, but by the whole spectral contents below that frequency.
Average amplification factor for [0.75, 3.0] f 0 and [0.75, 3.0] f 00 F L and F 0 . Let f 0 be the fundamental resonant frequency and f 00 ≡ min sites f 0 . Then the average amplification factors in the vicinity of f 0 and f 00 are defined using relations Such definitions are consistent with those the short and mid-period amplification factors (2-octave bandwidth), and are centred on 1.5f 0 and 1.5f 00 , respectively, instead of f 0 and f 00 , to account for the possibility of slight frequency shifts due to 2D or 3D effects (Bard and Bouchon 1985), and for the absence of symmetry around resonance frequencies, as amplification is larger at f 0 and beyond, than below f 0 .

Average amplification factor for a single-valued earthquake ground motion characteristic Let ,i ⃗
x be a single-valued earthquake ground motion characteristic at site ⃗ x and ψ ,i characteristic of the ith reference accelerogram. Then the amplification factor for characteristic ψ can be defined as and the average amplification factor for a set of n accelerograms as Appendix 3 defines characteristics of earthquake ground motion calculated in this study.

Numerical examples
Results of extensive numerical simulations and analysis for a set of local surface sedimentary structures is presented in the accompanying article by Moczo et al. (2018). Here we restrict to illustration of the described methodology on the example of the Grenoble valley. The Grenoble valley is a typical deep sediment-filled Alpine valley. Two aspects make it important: (1) Grenoble urban area with significant population, modern industry and research facilities. (2) Such "alpine valley" configuration occurs in different other areas within the European Alps, and in other mountainous areas with embanked valleys filled with young, post-glacial lacustrine sediments.
The numerical simulations of seismic motion are performed using the Fortran95 computer codes FDSim3D and FDSim2D . The computational algorithm is based on the (2,4) velocity-stress staggered-grid FD explicit heterogeneous scheme on the Cartesian discontinuous spatial grid. Here, (2,4) means the 2nd-order accuracy in time and 4th-order accuracy in space. In the FD method both medium and wavefield are represented by values in the discrete space-time grid. An explicit scheme for updating a particle velocity at a spatial position is obtained by a discrete approximation of the equation of motion and linear stress-strain relation formulated in the particle-velocity vector and stress tensor. The method was concisely described in the SIGMA deliverable D3-97 (Kristek et al. 2013). The basic references are the book by Moczo et al. (2014), and articles Etemadsaeed et al. (2016) and Kristek et al. (2017).

Model of the Grenoble valley, France
The Grenoble valley is a junction of three large valleys with complex geometry of the sediment-basement interface. The junction mimics letter Y. The sediments are made of the Quaternary fluvial and post-glacial deposits. The valley is surrounded by relatively high mountain ranges. The topography is neglected in this example (as it was shown to have only limited effects, see Chaljub et al. 2009;Tsuno et al. 2009;Garofalo et al. 2016). Figure 3 shows geometry of the sediment-basement interface. Table 1 lists values of mechanical parameters of the computational model (Chaljub et al. 2009(Chaljub et al. , 2010Garofalo et al. 2016).

Selected accelerograms
Following Sect. 2.2.1, 11 3-component accelerograms were selected from the RESORCE (Akkar et al. 2014) data base of accelerograms recorded on rock or stiff soil sites, in the near-source area (distance smaller than 40 km). The accelerograms have a very good signal-to-noise ratio over a wide frequency range with high-pass frequency smaller than 0.25 Hz. The accelerograms exhibit wide distribution of peak frequencies (F peak of the peak acceleration response spectrum), from around 1 Hz to beyond 16 Hz. The list of the accelerograms is given in Table 2 and the corresponding normalized spectra (PSA/pga) are illustrated in Fig. 4.

Computed acceleration and ground motion characteristics
In Fig. 5 we illustrate acceleration at a reference site (accelerogram 15560-Thjorarbru, Table 2), and accelerations at the two selected receiver positions along the P1 and P4 profiles across the Grenoble valley, P1rec and P4rec. The accelerograms at P1rec and P4rec indicate considerable site amplification and prolongation of ground motion. In Fig. 6 we show average amplification factors AF and their standard deviations, defined by Eqs. (25) and (26), respectively, at receivers P1rec and P4rec. We may note the obvious differences between the two receiver positions in the entire frequency range but the similar shape at frequencies larger than 3 Hz. We may point out surprisingly large amplification of the vertical component. (For discussion we refer to the companion article by Moczo et al. 2018) In Figs. 7 and 8 we show average amplification factors AF of CAV along P1 and P4, respectively. We again see large amplifications of the vertical component. There is clear difference between 2D and 3D simulations for the plane wave incidence in all components: the amplifications corresponding to 3D simulations are systematically larger than those corresponding to 2D simulations. For the horizontal components along P1, the amplifications for the point sources are closer to those resulting from 3D simulations for the vertically incident plane wave. For the vertical component along P1, the amplifications for the point sources are closer to those resulting from 2D simulations for the vertically incident  between 2D and 3D simulations for the plane wave incidence in all components along profile P4: the amplifications corresponding to 3D simulations are systematically larger than those corresponding to 2D simulations. Along the profile amplifications for the point sources are scattered with respect to the amplifications resulting from the 3D and 2D simulations assuming the plane-wave incidence. Comparisons for profile P1 are more complicated.

Conclusions
We developed methodology of calculating acceleration time history and corresponding characteristics of earthquake ground motion at a site of interest assuming acceleration at a reference site for two basic configurations: a reference site is a part of the model with a local surface sedimentary structure, a reference site is not part of the model with a local surface sedimentary structure. For each of the two configurations we assumed two wavefield excitations: a vertical plane-wave incidence and a point DC source. We illustrated the methodology of computation on the example of the Grenoble valley. and characterizes transfer properties of the model between the horizontal plane at which the excitation by the plane wave is applied and the reference site REF.

Matrix of the Fourier transfer functions for SIT
The matrix is defined as and characterizes transfer properties of the model between the horizontal plane at which the excitation by the plane wave is applied and the site of interest SIT.

Appendix 2
Reference site is in the model, excitation by a point DC source In the configuration shown in Fig. 12 we assume ground motion due to a point DC source.
Recall that an arbitrary point double-couple source can be obtained by a linear combination of 6 independent elementary sources.
Elementary sources e We consider the elementary sources represented by Eq. (6).

3
where φ denotes an average amplification factor of a characteristic of earthquake ground motion at the site, and 3D, 2D and 1D indicate dimension of a medium-wavefield configuration.