E El Haber
email: elias.el-haber@univ-grenoble-alpes.fr

C Cornou

D Jongmans

D Youssef Abdelmassih

Fernando Lopez-Caballero

T Al-Bittar

F Lopez- Caballero

Influence of 2D heterogeneous elastic soil properties on surface ground motion spatial variability

Keywords: Seismic response, spatial variability, probabilistic modeling

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.

Introduction

Spatial variation of earthquake ground motion (SVEGM) refers to the differences in amplitude and phase between recordings of the same earthquake at different locations. This variation is mainly caused by seismic source rupturing heterogeneities, regional wave propagation scattering and modification of ground motion by surface geology, most often referred as site effects [START_REF] Harichandran | Spatial variation of earthquake ground motion, what is it, how do we model it, and what are its engineering implications[END_REF]. One of the key factor controlling site effects is the presence of soft deposits overlying a more rigid geological formation. In sedimentary structures like valleys, soft surficial deposits may exhibit lateral variation in thickness, which lead to significant changes in surface ground motion even at close distance (e.g. [START_REF] Moczo | Wave diffraction, amplification and differential motion near strong lateral discontinuities[END_REF][START_REF] Field | Spectral amplification in a sediment-filled valley exhibiting clear basin-edgeinduced waves[END_REF][START_REF] Graves | Ground-motion amplification in the Santa Monica area: Effects of shallow basin-edge structure[END_REF][START_REF] Pagliaroli | Seismic microzonation of the central archaeological area of Rome: results and uncertainties[END_REF]. At local scale (tens to hundreds of meters), the spatial heterogeneity of near-surface material results from the surface natural processes of erosion and sedimentation [START_REF] Einsele | Sedimentary basins: evolution, facies, and sediment budget[END_REF] and human activities (construction, mining, …) generating anthropogenic deposits that can reach more than 10 m (e.g. [START_REF] Jongmans | The 1983 Liege earthquake: damage distribution and site effects[END_REF]Pagliaroli et al., 2014a). Additionally, intrinsic variation in properties within the soil layer contributes to the spatial variability [START_REF] Jenny | Factors of soil formation: A system of quantitative pedology[END_REF][START_REF] Burrough | Soil Variability:: A Late 20th Century View[END_REF]. The soil heterogeneity at various spatial scales may then have a significant effect on SVGEM, as frequently pointed out by the damage studies after earthquakes, especially on long-span structures such as dams, bridges and lifeline facilities [START_REF] Ariman | A review of the response of buried pipe-lines under seismic excitations[END_REF][START_REF] Anagnostopoulos | Pounding of buildings in series during earthquakes[END_REF][START_REF] Trifunac | Northridge, California, earthquake of 1994: density of pipe breaks and surface strains[END_REF][START_REF] Trifunac | Review: rotations in structural response[END_REF][START_REF] Bradley | Near-source strong ground motions observed in the 22 February 2011 Christchurch earthquake[END_REF]; for a summary of historic observations, refer to [START_REF] Kozák | Tutorial on earthquake rotational effects: historical examples[END_REF]. In particular, SVGEM has been put forward to explain the damage to the high rise buildings during the San Fernando earthquake in 1971 [START_REF] Hart | Torsional response of high-rise buildings[END_REF] and to the bridges during the earthquake of Loma Prieta in 1989 [START_REF] Kiureghian | Response spectrum me-thod for multisupport seismic excitations[END_REF]), of Northridge in 1994 [START_REF] Hall | Northridge earthquake[END_REF] and Christchurch in 2011 [START_REF] Chouw | Pounding damage to buildings and bridges in the 22 February 2011 Christchurch earthquake[END_REF].

There is therefore a need to quantify the soil parameter spatial variability and uncertainty for predicting SVGEM. To address this issue, probabilistic modeling approaches using random fields [START_REF] Popescu | Stochastic variability of soil properties: data analysis, digital simulation, effects on system behavior[END_REF] have been widely used in geotechnical engineering. In particular, such studies have been made on soil liquefaction (Lopez-Caballero and Modaressi, 2010), on the stability of soil dams and slopes [START_REF] Fenton | Statistics of free surface flow through stochastic earth dam[END_REF]Abdellah et al., 2000;[START_REF] Griffiths | Probabilistic slope stability analysis by finite elements[END_REF], on the analysis of shallow foundations (Youssef Abdel Massih and Soubra, 2008; Soubra and Youssef Abdel Massih, 2010; Al-bittar and Soubra, 2013 and Al-bittar and Soubra 2017) and on the analysis of equipartition time and diffusion regime in random heterogeneous media [START_REF] Khazaie | Influence of the spatial correlation structure of an elastic random medium on its scattering properties[END_REF]2017). In contrast, these techniques have aroused only limited interest among the seismological community, probably because of the difficulty of measuring small spatial scale variation in soil mechanical properties (Pagliaroli, 2014a;[START_REF] Salloum | The shear wave velocity structure of the heterogeneous alluvial plain of Beirut (Lebanon): combined analysis of geophysical and geotechnical data[END_REF]. The few 2D probabilistic approaches performed so far in earthquake engineering [START_REF] Assimaki | Effects of spatial variability of soil properties on surface ground motion[END_REF][START_REF] Nour | Finite element model for the probabilistic seismic response of heterogeneous soil profile[END_REF][START_REF] Thompson | Impediments to Predicting Site Response: Seismic Property Estimation and Modeling Simplifications[END_REF]Pagliaroli et al., 2014a) have highlighted the influence of the spatial variability of the soil properties on the surface ground motion. However, these papers mix effect of variable input ground motion with variability of linear and/or nonlinear soil properties on different surface ground motion indicators (spectral response, Peak Ground Acceleration, amplification, Housner intensity), making thus difficult the understanding of the respective impact of source of variability (input ground motion, linear and non-linear soil spatial variability) on SVGEM.

This paper aims at numerically studying the effect of the 2D soil spatial variability on SVGEM in the linear domain for a model consisting of one variable soil layer overlying an elastic bedrock, derived from a wellknown site (alluvial plain of Beirut, Lebanon). The chosen random soil parameter is the shear-wave velocity (Vs), the spatial variability being discretized by random fields using the Expansion Optimal Linear Estimation method (EOLE) [START_REF] Li | Optimal discretization of random fields[END_REF]. Synthetic seismograms at surface receivers are then computed using the FLAC2D Finite difference code. First, the effect of the 2D spatial variability of Vs at a single surface receiver is evaluated on the following scalar indicators: resonance frequency, corresponding amplification, as well as Arias based intensity and duration. Second, impact of spatially variable ground structure on the spatial correlation of surface ground motion is studied. Since some authors (Rahtje et al., 2010; Rodriguez-Marek et al., 2014; Haji-Soltani et al., 2017) chose to replace the 2D spatially variable seismic velocity profiles by 1D spatially variable ones for site-specific hazard assessment, we also compare surface ground motion indicators and related variability inferred from 2D and 1D modelling by using 1D Vs profiles extracted from the 2D Vs model.

Introducing ground structure spatial variability in wave propagation modeling

Strategy

Spatial variability of soil properties can be modeled by using the random field theory [START_REF] Vanmarcke | Random fields: analysis and synthesis[END_REF], defined by two functions (1) a Probability Density Function (pdf), which is usually positive-skewed and non-Gaussian in soils [START_REF] Popescu | Stochastic variability of soil properties: data analysis, digital simulation, effects on system behavior[END_REF], and (2) an autocorrelation function ρ(x,z). In this paper, the 2D shear-wave velocity structure, V s (x,z), is modeled as a random field with µ Vs and σ Vs being the mean and the standard deviation of the PDF, respectively. The coefficient of variation, COV= µ Vs /σ Vs , quantifies the range of dispersion of Vs around the mean. The autocorrelation function is defined by its spatial fluctuation scale described by the correlation distances θ x and θ z along the horizontal and vertical directions, respectively, which are the distances over which the soil property values are not correlated any more.

A number n of discretized 2D Vs models (called realizations) are generated using the probability density function and then used to simulate seismic wave propagation. Like any probabilistic modeling, n should be large enough to ensure the statistical convergence of the surface ground motion indicators (resonance frequency, spectral amplification, Arias intensity, duration …) in terms of average value and standard deviation. Among the numerous methods of random field discretization [START_REF] Sudret | Stochastic finite element methods and reliability: a state-ofthe-art report[END_REF], the Expansion Optimal Linear Estimation method (EOLE) belonging to the family of the series expansion methods (Li and Der [START_REF] Li | Optimal discretization of random fields[END_REF]) is used in this study. The main reason is that the convergence of ground motion indicators is reached after a lower number of seismic wave propagation simulations, compared to traditional Monte Carlo approaches [START_REF] Sudret | Stochastic finite element methods and reliability: a state-ofthe-art report[END_REF].

In order to understand near-surface spatial Vs variation effects on surface ground motion, we focus on a parametric sensitivity study for a simple case of a spatially variable soil layer overlying a bedrock with constant elastic properties. The Vs structure of the alluvial plain in Beirut (Lebanon) is chosen as a test site since detailed geotechnical and geophysical campaigns carried out at this site have allowed to image spatial variation of the underground geological structure (Saloum et al., 2014; Salloum, 2015).

Soil statistical parameters

In the approach adopted in this paper, the three statistical parameters describing the soil spatial variability are the coefficient of variation COV and the two autocorrelation distances θ x and θ z . Numerous studies attempted to determine these parameters from geotechnical laboratory or in situ tests (for a recent synthesis, see [START_REF] Salloum | Evaluation de la variabilité spatiale des paramètres géotechniques du sol à partir de mesures géophysiques: application à la plaine alluviale de Nahr-Beyrouth (Liban)[END_REF]. In contrast, very few studies were made on the quantification of the Vs variability. This probably results from the difficulty to obtain reliable Vs values through in situ tests.

In geotechnics, penetration tests (SPT -Standard penetration Tests or CPT -Cone penetration Tests) are widely used for soil investigation [START_REF] Cornforth | Landslides in Practice: Investigation, Analysis and Remedial/Preventative Options in Soils[END_REF]. They provide resistance parameters (N: number of blows for SPT; q c : cone resistance for CPT) that, like Vs, mainly increase with soil density or compactness if the effect of the effective stress due to overburden pressure is removed [START_REF] Karray | Influence of particle size on the correlation between shear wave velocity and cone tip resistance[END_REF]. For these two parameters, the coefficient of variation was found to cover a wide range between 5 and 43% (e.g Phoon and Kulhawy, 1999; [START_REF] Toro | Probabilistic models of site velocity profiles for generic and site-specific groundmotion amplification studies[END_REF] found COV ranging between about 30% and 45% for various NEHRP site classes. The autocorrelation parameters in a soil layer θ x and θ z can be estimated by using several methods (e.g. [START_REF] Popescu | Stochastic variability of soil properties: data analysis, digital simulation, effects on system behavior[END_REF][START_REF] Assimaki | Effects of spatial variability of soil properties on surface ground motion[END_REF]Thomson et al., 2009), depending on the available data and the chosen autocorrelation function. As most sedimentary layers are characterized by a strong anisotropy resulting from depositional processes, the horizontal autocorrelation distances are expected to be greater than vertical ones. Insitu geotechnical and Vs tests are usually performed along vertical profiles and most of the available information is on the vertical autocorrelation distance θ z . Thus, θ x is usually poorly constrained. Penetration test analyses provide a relatively wide range for the vertical autocorrelation distances, from a few cm to a few m, while the horizontal distances were estimated from around 1 m to more than 100 m (e.g., [START_REF] Alonso | Stochastic formulation of soil properties[END_REF][START_REF] Phoon | On quantifying inherent soil variability[END_REF][START_REF] Fenton | Estimation for stochastic soil models[END_REF] in the clay and gravel layers, respectively [START_REF] Salloum | Evaluation de la variabilité spatiale des paramètres géotechniques du sol à partir de mesures géophysiques: application à la plaine alluviale de Nahr-Beyrouth (Liban)[END_REF]. The vertical autocorrelation distance θ z was found to vary between 0.5 m and 2 m in the soil layers, using the same data. ERT (Electrical Resistivity Tomography)

profiles were used to evaluate the horizontal autocorrelation distance θ x , which was bracketed between 3.8 and 10.6 m [START_REF] Salloum | Evaluation de la variabilité spatiale des paramètres géotechniques du sol à partir de mesures géophysiques: application à la plaine alluviale de Nahr-Beyrouth (Liban)[END_REF], assuming that the variability is the same for the two geophysical parameters (shear wave velocity and electrical resistivity). For the purpose of simplifying the 2D modeling, the two soil layers were merged into a unique 15.5 m thick bed with a mean value µ Vs = 220 m/s overlying a homogeneous sound bedrock (Vs = 1000 m/s) (Figure 1a). The statistical parameters in the sedimentary layer are chosen in a range covering the values found during the survey and in the literature (5% < COV < 40%; 5 < θ x < 20 m; 1 < θ z <2 m). The 9 probabilistic models used for the sensitivity analysis are given in Table 1. No attenuation or damping is considered in the analysis in order to focus only on the effects of elastic properties variability. The fundamental resonance frequency of the soil layer for the mean Vs is f 0D = 3.54 Hz, with a theoretical amplification AF D = 7.1 at the resonance frequency. This case will be referred to as the deterministic model in the following.

Random field discretization method

The Expansion Optimal Linear Estimation method (EOLE), originally proposed by Li and Der [START_REF] Li | Optimal discretization of random fields[END_REF] is used herein to discretize the random field of Vs(x,z). It is described by a lognormal probability density function f (Vs) (Eq. 1):

( ) √ ( ( ( ) ) ) (1)
where  ln and σ ln are the mean and the standard deviation of Vs natural logarithm, respectively, and by an anisotropic square exponential autocorrelation function [( ) ( )] (Eq. 2), which gives the values of the correlation function between two arbitrary points (x, z) and (x', z'):

[( ) ( )] ( | | | | ) (2)
where θ x and θ z are the autocorrelation distances along x and z, respectively.

In this discretization method, one should first define a stochastic mesh composed of Nq grid points (or nodes)

and determine the autocorrelation matrix ∑ that gives the correlation between each grid point of the stochastic mesh and the other grid points using Eq. 2. The stochastic mesh dimensions ( ) are taken equal to half of the autocorrelation distance of the model in each direction ( , ). The random field Vs(x,z) at any point is estimated by the following equation, using the autocorrelation matrix :

̃( ) ∑ √ ( ) ∑ ( ) (3)
where ̃( ) is the approximate of the random field, and are the mean and standard deviation values of Vs, is the vector of independent standard gaussian random variables, and are the eigenvalues and eigenvectors of the autocorrelation matrice ∑ , ∑ ( ) is the correlation vector between the values of the random field at the different nodes of the stochastic grid and its value at the arbitrary point (x, z) as obtained using Eq. 2. The series of Eq. 3 is truncated to a number of terms s (expansion order), which is smaller than the number of grid points, after sorting the eigenvalues (j=1… S) in a descending order. This number should ensure that the variance of the error is smaller than a prescribed tolerance (10% in our study). The variance of the error on the random field in the EOLE method is given by (Li and Derkiureghian, 1993):

[ ( ) ̃( )] ∑ (( ) ∑ ( ) ) (4)
where ( ) and ̃( ) are the exact and the approximate values of the random fields at a given point (x, z), respectively.

Each vector of standard Gaussian random variables (j=1… S) provides (when substituted into Eqs. 3.3 and 3.4) spatial variations in Vs that fulfill the correlation structure of this field and are called Vs realizations. This is performed by computing, for this vector, the values of Vs at the centroids and the nodes of the different elements of the stochastic mesh. For each Since the finite difference mesh of FLAC 2D used to simulate wave propagation is much smaller than the stochastic mesh, the kriging method [START_REF] Sacks | Design and analysis of computer experiments[END_REF]) is used for determining V s at the centre of each finite difference mesh. In this study, the ordinary kriging and an anisotropic square exponential function for the correlation function are used. For each model (Table 1), the autocorrelation distances used for the kriging method are similar to the ones defined in the EOLE method. The lower and upper boundaries of the autocorrelation distances were taken equal to 0.5 m and 100 m, respectively. Thus, in both directions, 5 Vs are determined within one autocorrelation distance using the EOLE method and the others are extrapolated using the ordinary kriging method. By the combination of these two discretization methods, Figures 2d and 2h show example of one Vs realization for the 2 models #5 and #9 (Table 1), which have the same correlation distances (θ x =10 m and θ z =2 m) and differ by the COV value (20% and 40%, respectively). In order to check the accuracy of the generated 2D Vs realizations, the average shear wave velocity ( ) and the coefficient of variation ( ) are calculated for the two examples. The autocorrelation distances and are then evaluated respectively for the horizontal and vertical 1D samples extracted from the 2D realizations (refer to the black dashed rectangles in Figures 2d and2h). The re-estimated statistical parameters are summarized in Table 2. The results show that the generated 2D realizations have about the same statistical parameters initially introduced during the random field discretization phase. 2d and2h. For each realization, and are estimated, respectively, for the horizontal and vertical rectangle of each figure. 1) is shown in Figure 1a. The source time function is a pseudo-Dirac having a flat Fourier amplitude spectrum between 1 and 25 Hz (Figures 1b and1c). According to the frequency and V s ranges, the minimum and maximum wavelengths in the sediment layer are 2 m and 750 m, respectively. The mesh size Δl is chosen less than one tenth of the minimum wavelength to avoid numerical dispersion phenomena (Eq. 5) (Kuhlemeyer, R. L. and Lysmer, 1973):

Table 2: Statistical parameters: the average shear wave velocity ( ), the coefficient of variation ( ) and the autocorrelation distances ( , ) computed for the two V s realizations shown in Figures

Statistical parameters

(5

)
where ∆l max is the maximum size of the finite difference mesh, λ min is the minimum wave length, Vs min is the minimum shear wave velocity, and f max is the maximum wave frequency.

With the parameters of this study, the maximum mesh size, ∆l max , is 0.2 m ( ).

Zero horizontal displacements are applied along lateral boundaries of the model, while the horizontal and vertical movements are fixed at its base. Free field boundaries are applied to the side edges of the model. A quiet boundary (absorbing effect defined in FLAC 2D (Itasca, 2011)), is applied to the model base in order to model a flexible base that absorbs the energy emitted by the waves reflected on the surface and arriving to the model base. After initializing the stresses, a seismic shear stress excitation consisting of a vertically incident plane SV-type wave is applied to the base of the model:

( ) (6)
Where σ s is the applied shear stress, ρ is the soil density, V s is the medium shear wave velocity and υ s is the horizontal component of the velocity signal.

The 2D soil structure used in the simulation is a 22 m x 165 m model. The width of the model is relatively small compared to the maximum wavelength of 750 m. Therefore, in order to check that the boundary conditions are efficient and that such a model size does not affect the wave propagation, we compared seismograms computed for two different model sizes (750 m and 165 m wide) keeping fixed all the other parameters. The simulations performed over a duration of 2 seconds show similar results (Appendix 1). Finally, the computation time for one simulation and 10-second seismograms is around 5 hours on a single processor and a PC having an i7-core and 2.4 Hz CPU frequency, leading to a total computation time of about 4500 hours for all models.

Characterization of the ground motion spatial variability

In this section, we present an analysis of the synthetic seismograms computed for two probabilistic models (one realization), discuss the probalistic simulation stability and provide the definition of the ground motion indicators extracted from the seismograms. and θz=2m) and #9 (COV=40% θx=10m and θz=2m) (see Table 1). The white dashed lines in (a) and (e) correspond to the fundamental frequency of the deterministic model, f oD .

3.1. Ground motion synthetics and probabilistic simulation stability

Ground motion synthetics

Figure 2 shows the simulation results of one Vs realization for the 2 models #5 and #9 (Table 1). For each case, the Vs model, the vertical and horizontal surface velocity signals, and the Fourier amplitude spectra of the horizontal seismograms are shown. The increase in COV (40%) is clearly visible on the Vs model #9 (figure 2h), which exhibits larger velocity contrasts than model #5 (Figure 2d; COV= 20%). The Vs spatial variation turned out to have a strong influence on the surface ground motion characteristics: seismograms exhibit larger amplitude and longer duration at receivers located over near-surface low velocity zones (see Figures 2b,c, f and g), especially for the horizontal seismograms and for COV=40%. The Fourier amplitude spectra (Figures 2a and2e) show that the amplification at the fundamental resonance frequency occurs at or a little lower the deterministic 1D resonance frequency (f 0 =3.54 Hz), with lateral variations resulting from the surface wave propagation. Amplification at higher frequencies are also spatially variable (Figure 2e). In order to determine the wave types, the particle motions are plotted in Figures 3c and3d vanishing below the soil fundamental frequency (see also the Fourier amplitude spectrum in Figure 2a).

Although the theoretical curves from the deterministic model (first three modes, black curves) do not perfectly match the dispersion maxima caused by Vs differences between the 2D and the deterministic cases, the comparison suggests that the first two Rayleigh wave modes have been excited in the 2D probabilistic model.

Convergence of probabilistic ground motion simulations

Probabilistic modeling approaches require simulating ground motion for a large enough number of realizations of discretized Vs models to ensure the statistical convergence of the average and standard deviation estimators of any surface ground motion parameter. We define the convergence of a ground motion parameter α as a fraction of the relative variation of α between realization i and realization i+1:

( ) | | | | (7)
The convergence is considered to be reached when CONV(α) < 5% [START_REF] Haldar | Effect of soil spatial variability on the response of laterally loaded pile in undrained clay[END_REF]. For both the average and standard deviation estimators, the convergence was tested on both surface ground motion parameters in time (Arias Intensity and duration) and frequency (Fourier amplitude spectra, site fundamental frequency and corresponding amplification) (see section 3.2 for definition). As an example, the convergence for the Fourier amplitude spectra (horizontal component) at different frequencies and the duration are shown in Appendix 2 for model #9 that exhibits the largest Vs range. According to the defined criterion, the convergence is reached after 86 realizations for all frequencies for Fourier amplitude spectra and after 82 realizations for the duration. All the tests performed on other indicators showed that 100 realizations ensure the convergence. Finally, due to the presence of absorbent boundaries, receivers located less than 30 m from the borders were omitted in the analysis.

Surface ground motion indicators

In this paper, we consider surface ground motion indicators in time and frequency domain, as described hereafter. In the frequency domain, the spectral amplification AF(f) is defined as the ratio between the Fourier amplitude spectra of the signals recorded at the soil surface and at the outcropping bedrock. For vertically incident waves, this latter is 2 times the spectrum shown in Figure 1c. Two parameters are extracted from the spectral amplification: the fundamental frequency f 0 and the corresponding amplification, AF f0 . In the time domain, we consider the Arias based intensity parameter (A b I) that provides a measure of the shaking intensity [START_REF] Arias | Measure of earthquake intensity[END_REF], here defined as the integral of the square surface horizontal velocity as function of time. The second time parameter is the effective duration (D AbI ) defined as the difference between the time where 5% and 95% of the total A b I is reached. The durations inferred for models #5 are indicated in Figure 3a.

Influence of ground structure variability on surface ground motion indicators at a single station

Resonance frequency and related amplification

The amplification AF was calculated at the central receiver (located at x=83 m) for all models (Table 1) from the 10-second seismograms after smoothing the Fourier amplitude spectra by using the Konno and Ohmachi algorithm with b=50 (Konno and Ohmachi, 1998). Spectral amplifications for the 100 simulations are shown in Figures 4a and4b for models #5 (COV=20%) and #9 (COV=40%), respectively. The variability of the amplification curves increases with COV, both in terms of maximum amplification and corresponding frequency. We computed the probability density functions of f 0 and AF fo (Figures 4c and4d) using the nonparametric kernel density estimation technique [START_REF] Rosenblatt | Remarks on some nonparametric estimates of a density function[END_REF][START_REF] Parzen | On estimation of a probability density function and mode[END_REF]. Compared to the deterministic fundamental frequency f 0D (dashed red line), the probabilistic f 0 values tend to be lower than f 0D by about 2% and 20% for COV=20% and 40%, respectively (Figure 4c). The distribution of amplifications at the fundamental resonance frequency peaks (AF fo ) (Figure 4d) is slightly higher by about 6.7% than the deterministic one (AF foD ) for COV = 20%, with AF fo values ranging between 3.15 and 12.0, while it is significantly shifted to a lower amplification values for COV=40% with AF fo values ranging between 1.6 and 17.2.

Figure 4 : Amplification functions at the central receiver. (a) and (b) Amplification curves (100 realizations in grey) for parameter set #5 (COV=20% θx=10m and θz=2m) and for the parameter set #9 (COV=40% θx=10m and θz=2m), respectively. The black curve is the theoretical transfer function for the deterministic model (Vs=220 m/s). (c) KSdensity of the fundamental resonance frequency f0 computed by considering all realizations for parameter set #5 (in blue) and #9 (in black). The red dashed line indicates the fundamental resonance frequency f0D of the deterministic model. (d) KSdensity of the amplification AFf0 at the fundamental resonance frequency for parameter set #5 (in blue) and #9 (in black). The red dashed line indicates the amplification AFD for the deterministic model.

The average amplification and the standard deviation computed on the natural logarithm of amplifications at the central receiver (x= 83 m) are shown for all models in Figures 5a and5b, respectively. The average AF curves provide fundamental resonance frequencies close to the deterministic one, except for the higher COV (40%). For that COV, the fundamental frequency f 0 is shifted down to values around 3 Hz, as already highlighted in Figure 4c. At frequencies equal or over f 0 , standard deviations increase with frequency, the higher standard deviation being found for larger COV values as a consequence of stronger scattering of the AF curves (Figure 5b), which in turns leads to a decrease in average amplification values at frequencies larger than f 0 (Figure 5a). In order to evaluate the effect of the receiver position on amplification, we calculated the standard deviation σ(ln(AF)) for all receiver locations. The curves are presented in Figure 5c at the same scale as the ones in Figure 5b. The standard deviation is slightly affected by the receiver location at frequencies higher than 15 Hz for the greatest COV (40%). The main impact of considering all receivers is however observed for frequencies lower than f 0 , a range for which the standard deviation is significantly greater, whatever the COV value. Such increase results from the filtering effect of the sedimentary layer, which yields little energy below the fundamental frequency [START_REF] Scherbaum | Determination of shallow shear wave velocity profiles in the Cologne, Germany area using ambient vibrations[END_REF] and greater variability in amplification. Moreover, Appendix 3 compares the standard deviation of ln(AF) computed by considering, firstly, all the receivers of all the realizations (in dashed lines) and secondly, the central receiver of each realization (in continuous lines) for the model #6 having COV=20%, θz=2m and θx=10m and model #9 having COV=40%, θz=2m and θx=10m. The standard deviations calculated by the two methods are very similar. This observation emphasis the stationarity of ln(AF) due to the stationary random field generated by the EOLE discretization method.

The simulation results are synthesized for the 9 probabilistic models in Figure 6. First, the average f 0 values along with error bars for all parameter sets are shown as a function of COV at the center receiver in Figure 6a.

As already evidenced in Figure 4 for two models, the average f 0 at the central receiver is close to the deterministic value f 0D = 3.54 Hz for low COV but sharply decreases to around 2.8-3 Hz for COV=40%. The two other probabilistic parameters (θ x and θ z ) do not seem to significantly influence neither the average nor the dispersion of f 0 . This significant influence of COV on average f 0 values is observed at all receivers (Figure 6e), θ x and θ z having little effect. The variations in f 0 over distance are small et remain below 10%.

The same kind of analysis are presented for the average amplification AF fo in Figures 6b andf. At the central receiver (Figure 6b), the average amplification is close to the deterministic value AF D and, again, θx and θz seem to have little influence. However, the average amplification variability increases largely for COV=40% with values being higher or lower than AF D , depending on the combination of θx and θz. Figure 6f shows that, in contrast to f 0 , AF fo is spatially variable, especially for COV=40%. Such spatial variation is due to large lateral Vs variation and/or to the presence of a double amplification peak close to the fundamental resonance frequency on the AF curve at some receiver location. The second peak generated by low Vs superficial zones may exhibit an amplitude higher than the first one identified as the amplification at the fundamental resonance frequency (see Figure 2e for receivers located between X=30 m and X=120 m).

Arias based intensity and duration

Similarly to the spectral parameters, the results for temporal parameters (Arias intensity A b I and duration A b ID) are synthesized in Figure 6c, d, g andh. The average Arias intensity and duration (Figures 6c and6d, respectively) are systematically larger than those corresponding to the deterministic Vs structure, and increase with COV, as well as the standard deviations. This effect is particularly clear for A b ID (Figure 6d), which may reach values two to three times higher than the deterministic one for COV=20% and 40%, respectively. Also, these two temporal parameters seem to be slightly influenced by the receiver location, as shown in Figures 6g and6h, except A b I for COV=40%. In this case, large lateral variations of 50% can be observed over short distances.

In summary, this study outlines the predominant influence of COV on the surface ground motion indicators in terms of average and standard deviation compared to the correlation distances. However, for the largest COV (40%), ground motion average indicators related to the amplitude of seismic signals (amplification AF fo , Arias intensity A b I) exhibit larger variation over distance for θx=5 m and θz=2 m suggesting then an influence of the size of the heterogeneity, while the average fundamental frequency f 0 and duration A b ID turned out to be spatially independent, whatever the COV value.

Effects of spatial ground variability on spatial correlation of indicators

In previous sections, the variability of indicators was analyzed at a single station. However, studying the impact of spatial ground variations on the spatial correlation of ground motion is also of interest for understanding or predicting damage to long-span civil engineering structures (e.g. [START_REF] Abrahamson | Spatial variation of strong ground motion for use in soil structure interaction analyses[END_REF]Schneider et

Scalar indicators

We first computed the mean spatial correlation of the 4 scalar indicators for the 100 Vs realizations. Figures 7a,7b, 7c and 7d display the spatial correlation of the resonance frequency (Δf o ), the corresponding amplification (ΔAFf o ), the Arias based intensity (ΔA b I) and the Arias based intensity duration (ΔA b ID), respectively, as a function of ΔX/θ x for all the models (Table 1).

The variability of ΔA b I and ΔA b ID increases linearly until reaching a constant value for ΔX larger than 2θ x , a larger variability being systematically observed for larger COV values. This general trend is explained by the random field discretization in Vs (see section 2). For a given realization, the Vs vertical profiles below receivers located at X+/-θ x are correlated, leading to close seismic responses at the surface. In contrast, the seismic responses at two receivers with a spacing larger than 2θ x are not expected to be similar. Thus, the seismic indicators ( A b I and A b ID) that are strongly controlled by the most energetic seismic phase (the incident SV wave (Figure 3a), mostly itself influenced by the locally 1D Vs structure) will lead to indicator differences not correlated anymore for ΔX larger than 2θ x and thus constant whatever the receiver spacing.

The spatial correlation of the fundamental resonance frequency (Δf 0 ) and the corresponding amplification 7a, 7b andAppendix 4) indicates an increase in variability with distance, reaching constant or maximum values for ΔX/θ x ranging between 2.5 and 4 for θ x =20 m, between 6 and 8 for θ x =10 m and between 13 and 16 for θ x =5 m. These ratios correspond to an absolute value ΔX ranging between 50 and 80 m, i.e. close to the wavelength at the fundamental resonance frequency for the deterministic Vs structure (λ=62 m). As for other indicators, the f 0 and AFf o variability increases with COV and their convergence seems to be independent of θ x like the time domain parameters (see also Appendix 4 where results are plotted in terms of non-normalized ΔX). 8a to 8d). For small frequencies (<f 0D ), the wavelengths are much larger than θx (5 and 10 m) and θz (1 and 2 m) and the analyzed inter-receivers distances (ΔX). Then, σ(Δln(AF)) increases with the increase of ΔX and no convergences is observed with ΔX in the range between 0 and 100m. The correlation distances θx and θ z appear to have little influence on the variability, except for one frequency (7 Hz) for COV=40%. The spatial amplitude variability with frequency shows a different picture depending on the receiver spacing ΔX (Figures 8e and8f). For ΔX =5 m, the variability is weak at low frequency and progressively increases with frequency (see for instance the curve for COV=40%). This shows the effect of the wavelength. Indeed, for COV=40%, V s varies from 50 to 750 m/s, corresponding to minimum wavelengths of 2, 5, 7 and 14 m for V s = 50 m/s and at 24, 10, 7 and 3.54 Hz, respectively. For wavelengths much larger than the size of the heterogeneity (i.e. at low frequency) the amplitude variability is low, while maximum when wavelength is about the size of the heterogeneity. It could explains that maximum variability is observed for θx=5m at high frequency (Figure 8e). For large spacing (ΔX =50m > θ x .), the Vs structure below the receivers is not correlated anymore and σ(Δ|ln(AF)|) is maximum at the resonance frequencies (f 0D , 3 f 0D and 5 f 0D ), with no clear dependency on θ x (Figure 8f).

(ΔAFf o ) (Figures

Frequency-dependent indicator

Can 2D modeling be replaced by 1D modeling?

In order to evaluate the importance of accounting for 2D wave propagation on surface ground motion for site specific hazard assessment compared to alternative modelling approaches that simplify the 2D spatial variation of ground structure into a suite of 1D variable soil profiles (e. showing the parameter ratios and standard deviation differences for all realizations at the central receiver, the spatial evolution of the four parameters along the profile for one given realization , and the average and standard deviation of 2D/1D ratio for all realizations, respectively 6.1. Fundamental frequency, corresponding amplification and spectral amplification Figure 6.1a and b compare in terms of ratio the average fundamental resonance frequency and amplification together with related standard deviation derived from 1D and 2D modeling approaches. Average resonance frequencies are very similar for both modeling types (Figure 9a), while standard deviation is larger for 1D modelling (Figure 9b). Such higher standard deviation comes from 1D resonance frequencies being only controlled by SH body waves propagating through locally 1D soil profiles, while 2D wave propagation modeling also incorporates locally diffracted surface waves that sample larger underground volume and lead to homogenize seismic response. Figure 10a illustrates, for the ground model realization shown in Figure 2d, differences between spatial distribution of resonance frequency derived from 1D and 2D wave propagation modeling. The resonance frequency computed from the 2D ground motion synthetics is very weakly sensitive to the spatial variation of ground velocity structure compared to that computed from 1D synthetics, except when strong lateral discontinuities of Vs occur (e.g. at X=65 m in Figure 10e). Although average amplification computed from 1D and 2D ground motion synthetics are very similar, especially for COV of 5% and 20%, standard deviation of amplification computed from 2D synthetics are significantly larger (by 15% to 40%) than the one inferred from 1D synthetics, whatever the COV value. Such a large increase in standard deviation is easily explained when looking at the spatial variation of 2D amplification for a single Vs realization (Figure 10b). Amplification variation is either due to large lateral Vs variation prone to efficiently generate (or not) surface waves (Figures 2b and10b) or the presence at some receivers of a double amplification peak close to the resonance frequency, leading in most cases to an amplitude at the fundamental resonance frequency much lower than the one at the secondary peak (see Figure 2e for receivers located between X=30 m and X=120 m).

In order to analyze the effect of 2D wave propagation modeling approach compared to the 1D one over the whole frequency range from 1 Hz to 25 Hz, we computed the ratio between Fourier amplitude spectra obtained from 2D and 1D ground motion synthetics. Accounting for 2D wave propagation modeling in spatially variable Vs structure obviously leads to larger amplification at frequencies higher than the fundamental resonance frequency of the probabilistic model, except for the frequency range in-between the fundamental and first higher mode (Figure 11a). This increase in amplification, especially significant for COV=20% and COV=40% (Figure 11a), is also associated with a higher variability (Figure 11b), which increases with COV.

These observations clearly outline that, although the 1D probabilistic approach may correctly reproduce average fundamental resonance frequencies and corresponding amplification, it under-predicts ground motion amplification at high frequencies and related variabilities whose origin relates to the presence of locally diffracted surface waves. 1 (COV=20% θx=10m and θz=2m).

Arias based intensity and duration

Average Arias based intensity and duration inferred from 2D ground motion synthetics are systematically higher than the ones inferred from 1D synthetics, by a factor of 1.2 to 1.6, whatever the considered COV and autocorrelation distances (Figures 9c and9d). The increase in intensity and duration for 2D modeling is simply explained by the presence of the locally generated surface waves that contaminates the strongest S-wave phase used in the computation (e.g. see Figure 2b for receivers close to X = 50 m, X = 100 m and X = 130 m and the corresponding Arias intensity increase intensity in Figure 10c). As a consequence, standard deviation of Arias based intensity is larger for 2D synthetics compared to 1D ones. Interestingly however, the relative difference between standard deviation of duration is decreasing with COV, reaching values below one for COV=40%. This feature is explained by the increase of the duration standard deviation in 1D modelling due to the presence of spatial variability in A b I computed from 2D synthetics is greater than the one derived from 1D synthetics, while the opposite is observed for A b ID. The larger variability in A b I (Figure 7c) results from the presence of body and surface waves diffracted at near-surface heterogeneities in the early seismic phase (see Figure 3a). The lower variability in AbID (Figure 7d) most probably results from the presence of diffracted surface waves in the analysis time window leading to large duration as shown in Figure 7d.

Discussion and conclusion

In this paper, the effect of the Vs spatial variability on the surface ground motion is assessed through a set of numerical computations. A simple elastic 2D velocity structure is used for surface ground motion simulations in order to solely focus on the effects of Vs variability. Vs is modeled as a random field, considering the coefficient of variation (COV) and the horizontal and vertical autocorrelation distances (θx and θz, respectively) as parameters.

We show that spatial variation of shear-wave velocity structure efficiently generate locally diffracted surface waves. Hence, the fundamental resonance frequency is only very weakly sensitive to the local Vs heterogeneities while the Arias Intensity and derived duration are clearly increased due to the locally diffracted body and surface waves that contaminate the most energetic phase. Analysis of synthetic seismograms at a single station clearly outlines that the coefficient of variation on the shear-wave velocity (COV) mostly controls the variability of the four ground motion indicators (resonance frequency and corresponding amplification, Arias intensity and duration), as compared to the horizontal and vertical correlation distances. As regards spectral amplifications, considering shear-wave spatial variability in the modeling leads to reduce mean amplification, and conversely to increase variability, at frequencies higher than the site fundamental resonance frequency, as also observed in [START_REF] Assimaki | Effects of spatial variability of soil properties on surface ground motion[END_REF]. A striking result is the 20% decrease in mean fundamental resonance frequency (f o ) for COV = 40% (Figures 5a and6a) compared to the resonance frequency of the deterministic model. This reduction of f 0 , similarly observed on 1D and 2D ground motion probabilistic modeling, is related to the Vs random field discretization. Indeed, we calculated the SV waves travel time for the 1D soil columns extracted from 2D Vs realizations for two parameter sets with COV = 20% and COV=40%, respectively (probabilistic models #5 and #9, Table 1). The SV wave travel time for the deterministic sedimentary layer is t d = 0.0705 s. The travel time for the parameter set having a COV= 20 % (Figure 12a) is found to be almost normally distributed around t d , while log-normally distributed for the parameter set with COV = 40 % (Figure 12b). In this latter case, a larger number of realizations has travel times exceeding the t d , leading to a reduction of the mean fundamental resonance frequency for the 1D ground motion synthetics.

Although the calculation of travel times for heterogeneous 2D soil models from 1D soil columns is approximate, these results suggest that long travel times in 2D realizations with COV=40% may be the cause of the resonance frequency drop. Like other authors for 2D probabilistic modeling [START_REF] Assimaki | Effects of spatial variability of soil properties on surface ground motion[END_REF][START_REF] Nour | Finite element model for the probabilistic seismic response of heterogeneous soil profile[END_REF]2003;[START_REF] Thompson | Impediments to Predicting Site Response: Seismic Property Estimation and Modeling Simplifications[END_REF] Comparison between spectral amplification obtained from 2D and 1D probabilistic modeling synthetics clearly outlines that, although both modeling approaches predict similar mean fundamental resonance frequency and corresponding amplification, 2D wave propagation modeling leads to larger amplification and related variability at frequencies higher than the fundamental resonance frequencies, especially for large COV, compared to 1D wave propagation probabilistic modeling approach (Figures 9 and11). This increase in amplification and standard deviation is definitely caused by the presence of locally diffracted surface waves in 2D wave propagation modeling and, hence, questions the reliability of using a population of 1D Vs profiles to account for spatial variability of the seismic properties as proposed for site-specific hazard assessment (Rathje et 1 whatever the site condition, i.e. larger than the constant value obtained in our numerical simulations. This may be due to the fact that our modeling is not accounting for 3D wave propagation, attenuation and broadband seismograms.

To sum up, our study clearly outlines the main control of COV on the spatial variability of surface ground motion indicators and that stochastic ground motion modeling by using 1D variable soil profiles to account for spatial variation of ground structure properties will underestimate surface ground motion spatial variability as compared to full 2D/3D wave propagation modeling. However, the simple case of a sedimentary layer over a bedrock does not represent the real site conditions. To confirm our results, further numerical simulations should be carried out for a set of typical ground structure, considering also attenuation and non-linear soil behavior.

Appendices

Figure 1 :

 1 Figure 1: (a) Two-layered ground model composed of a 15,5m thick soil layer overlying bedrock (Vs=1000 m/s). The Vs distribution in the soil layer was computed using the EOLE method with the following statistical parameters: µVs = 220 m/s; COV= 40%; θx= 10 m; θz =2 m (model # 9 in Table1). refers to the distance of the receiver from the left side of the model. refers to the distance between receiver and . Synthetic

Figure 2 :

 2 Figure 2: Seismic responses for the two Vs realizations of models #5 (left column) and # 9 (right column). (a) and (e) Fourier amplitude spectra (FAS) for the horizontal velocity seismograms shown in (b) and (f); (b) and (f) 3-second horizontal seismograms (H) computed for the Vs realizations shown in (d) and (h), respectively. (c) and (g) 3-second vertical seismograms (V) ; (d) and (h) Vs realizations for models #5 (COV=20% θx=10m

Figure 3 :

 3 Figure 3: (a) 3-second seismograms (horizontal and vertical velocities) simulated for the Vs model shown in Figure 2d at 3 different locations spotted by black arrows in Figure 2d (curves in black). Seismograms computed for the deterministic model are indicated in red. The vertical green lines represent the duration corresponding to 5 and 95% of the total energy calculated for the seismograms computed for the probabilistic model (see text for details). (b) Dispersion curves for the vertical component of the synthetic ground motion obtained for distance between 65 and 115 m between t=0 s and 8 s (refer to rectangle A1 in Figure 2c); the theoretical dispersion curves for the first three modes of the deterministic model are shown in black. (c) Particle motion of the 1-5 Hz band-pass filtered velocities recorded at X= 83 m between t=0.1 s and t=0.8 s (top) and 7-11 Hz band-pass filtered velocities recorded at X = 83 m between t=1s and 1.8s (bottom).

Figure 3a shows the

  Figure 3a shows the first three seconds of the ground surface seismograms (horizontal and vertical components) calculated for the model #5 at three distances indicated in Figure 2d. The horizontal signals are compared with the ones computed for the deterministic model considering no variability in the soil properties (µ Vs = 220 m/s; σ Vs =0).Figure 3a illustrates the motion spatial variability generated by the 2D probabilistic models in terms of amplitude of the first arrivals (e.g. see the signals at 45 m and 116m) and generation of strong later seismic waves (e.g. at 45 m) both on the horizontal and vertical components.

  Figure 3b indicates two branches of dispersion curves between 4Hz and 20 Hz, while the image is blurred below 4 Hz. This results from the high-pass filtering of the soil layer (Scherbaum et al. 2003), which leads the energy

Figure 5 :

 5 Figure 5 : (a) Average amplification at the central receiver (X=83 m) for different values of COV, θx and θz obtained by computing the geometric mean of the transfer functions for the 100 realizations. The amplification and the fundamental frequency for the deterministic model are shown by the black curve and black dashed line, respectively. (b) Standard deviation of the natural logarithm of the Amplification (σln(AF)) at the central receiver. (c) Standard deviation of σln(AF) obtained by considering all receivers location. The values of COV, θx and θz are shown with different colors, symbols and lines, respectively.

Figure 6 :

 6 Figure 6: Synthesis of simulation results for the 9 parameter sets. Left panel: Average values and error bars (+one standard deviation) for the four spectral and temporal parameters at the central receiver: (a) average fundamental frequency (fo), (b) average amplification at the fundamental frequency (AFfo), (c) average of the Arias based intensity (AIb) and (d) average of the Arias based duration (AbID). The black dashed lines indicate values of, f0, AFfo, AIb and the related duration of the deterministic model. Right panel: (a) average f0, (b) average AFfo, (c) average AbI and (d) duration at different locations on the surface for the 100 simulations. The values of COV, θx and θz are shown with different colors, symbols and lines, respectively.

  al., 1992; Ancheta et al., 2011; Goda and Atkinson, 2008; Liu and Hong, 2015; Koufoudi et al., 2018). Here, we focus on investigating the spatial correlation of the four ground motion scalar indicators defined in section 4 (resonance frequency and corresponding amplification, Arias based intensity and duration) and one frequencydependent indicator (Fourier amplitude spectrum). We estimated the spatial correlation by computing the absolute differences of indicators between two receivers separated by a distance ΔX, ΔX spanning all possible inter-receiver distances and being normalized by the horizontal correlation distance θ x .

Figure 7 :

 7 Figure 7: Spatial variability of the four scalar indicators as a function of the receiver spacing normalized by the horizontal autocorrelation distance ΔX/θx. (a), (b), (c) and (d), Variation of the differences in fundamental resonance frequency f0, related amplification AF, Arias based Intensity A b I and Arias based Intensity duration A b ID using 2D ground motion synthetics. (e), (f), (g) and (h), the same for 1D ground motion synthetics. The values of COV, θx and θz are shown with different colors, symbols and lines, respectively.

Following

  [START_REF] Abrahamson | Spatial variation of strong ground motion for use in soil structure interaction analyses[END_REF] and[START_REF] Schneider | The spatial variation of earthquake ground motion and effects of local site onditions[END_REF], we computed the standard deviation of the absolute difference of the natural logarithms of the Fourier amplitude spectra (σ(Δ|ln(AF)|)). Figures8a, 8b, 8c and 8d display σ(Δ|ln(AF)|) versus the receiver spacing ΔX for four frequencies (3.54 Hz, 7 Hz, 10 Hz and 24 Hz, respectively), while Figures 8e and 8f show σ(Δ|ln(AF)|) as a function of frequency for two receiver spacing values (5 and 50 m, respectively). The lowest frequency (3.54 Hz) corresponds to the fundamental resonance frequency f 0D previously defined. For all frequencies, σ(Δ|ln(AF)|) increases with ΔX until reaching a constant value that significantly increases with COV (Figures

Figure 8 :

 8 Figure 8: Frequency-dependent indicator variability computed from 2D ground motion synthetics. The standard deviation of the absolute difference of natural logarithm of Fourier amplitude spectra (σ(Δln(AF))) as a function of receiver spacing at frequencies equal to (a) 3.54 Hz, (b) 7 Hz, (c) 10 Hz and (d) 24 Hz. σ(Δln(AF)) as a function of frequency for two receiver spacing values (d) ΔX=5 m and (f) ΔX= 50 m. The values of COV, θx and θz are shown with different colors, symbols and lines, respectively.

  g. Rahtje et al., 2010; Rodriguez-Marek et al., 2014; Haji-Soltani et al., 2017), we extracted 100 1D soil columns at each receiver from the Vs 2D discretization. 1D wave propagation simulation was performed and spectral amplification, Arias based intensity and duration are computed at each receiver and compared to the values obtained from 2D wave propagation modeling. The 1D and 2D modeling results are compared from various perspectives in Figures 9, 10 and 11,

Figure 9 :

 9 Figure 9: Ratio between the parameter extracted from 2D wave propagation modeling and 1D ones at the central receiver (left column), and the relative difference between standard deviation of this parameter inferred from 1D and 2D wave propagation modelling (right column): (a) and (b), average fundamental resonance frequency (c) and (d) average amplification at the fundamental resonance frequency, (e) and (f) average of the Arias based intensity, (g) and (h) average of Arias based duration.: (e) fundamental resonance frequency, (f) average amplification at the fundamental resonance frequency, (d) average of the Arias based intensity, e) average of Arias based duration. The values of COV, θx and θz are shown with different colors, symbols and lines, respectively.

Figure 10 :

 10 Figure 10: Comparison of 2D (black line) and 1D modeling (blue line) results for the Vs realization shown in (e). (a) Fundamental resonance frequency (f 0 ), (b) amplification at f 0 (AF f0 ), (c) Arias based Intensity (A b I) and (d) Arias based Duration. The test parameter is #5 inTable 1 (COV=20% θx=10m and θz=2m).

Figure 12 .

 12 Figure 12. Travel-time distribution (blue thick line) in the soil column below the central receiver for the 100 realizations for (a) model #5 (COV=20%) and (b) model #9 (COV=40%) listed in Table 1. The travel time computed from the deterministic model (td=0.0705 s) is shown a black dashed line.

  or 1D modeling (Li and Assimaki, 2010; Rathje et al., 2010; Rodriguez-Marek et al., 2014; Haji-Soltani et al., 2017), we discretized Vs because it is the physical parameter that can be measured in the field. Considering however that seismic wave propagation is primarily controlled by travel time, random field discretization on travel time would be the most natural and relevant choice. Although such discretization is easily achievable for 1D body or surface wave propagation modeling (e.g. Pilz and Parolai, 2014; Teague and Cox, 2016; Pilz and Fäh, 2017; Teague et al., 2018), it remains more complex for body and surface wave propagation in 2D heterogeneous medium.

  al., 2010; Rodriguez-Marek et al., 2014; Haji-Soltani et al., 2017). Our results also indicate that the indicator spatial variability increases linearly until reaching a constant value, with greater variability being systematically observed for larger COV. Effect of horizontal correlation distance, θx, is only observed for frequency dependent (and thus wavelength dependent) indicators (Fourier amplitude spectra) when receiver spacing, wavelength and lateral size of heterogeneity (θx) are of the same order. Interestingly, ground motion indicators (Arias Intensity and derived duration) computed on the strongest seismic phase dominated by S waves are not anymore spatially correlated for receiver spacing larger than two times the horizontal correlation distance. This opens up interesting perspectives in measuring horizontal autocorrelation distances from available seismological recordings at dense array data. Maximum of the spatial correlation of the natural logarithms of the Fourier amplitude spectra obtained on real data (Abrahamson et al., 1990; Schneider et al., 1992; Ancheta et al., 2011; Goda and Atkinson, 2008; Liu and Hong, 2015; Koufoudi et al., 2018) is close to

Appendix 1 :

 1 Synthetic horizontal surface obtained for a homogeneous layer over a half-space for a vertically incident SV plane wave excitation at the base of the homogeneous layer. Source time function is a delta-like with flat Fourier amplitude spectrum between 1 and 25 Hz. Shear-wave velocity in the homogeneous layer is 750 m/s and 1000 m/s in the halfspace. The computational model is 15.5 m along depth axis and (a) 165 m and (b) 750 m width along the horizontal direction. In (b), the computed seismograms are only shown for the 265 central receivers out of 750 receivers.

  Srivastava and Babu, 2009; Zhao et al., 2018). For getting the coefficient of variation (COV) of Vs, Thomson et al. (2007) compiled results of Vs measurements using different techniques (down-hole tests, seismic cone tests) in various geologic units (Holzer et al, 2005; Wills and Clahan, 2006; Moss, 2008). They found that COV values are generally contained within the interval from14 to 46% for Vs measurements in the same geologic unit. COV values in a similar range are also reported by Haji-Soltani and Penscheck (2017) for a very soft sediment (5 to 15%) and by Rodriguez-Marek et al. (2014) for very hard rock site (10 to 25%). Making use of the Electric Power Research Institute (EPRI, 1993) Vs profile database,

Table 1 :

 1 Properties in the two layers (soil and bedrock) for the probabilistic modeling. The bedrock properties are fixed, as well as the density ρ and the P-wave velocity Vp in the soil layer. The soil S-wave velocity is considered as a random field and characterized by four statistical parameters (µ Vs , COV, θx, θz). Nine models are defined with various COV, θx and θz (see text for details). The so-called deterministic model is defined by the values shown in italics in the Table.

  ; Jaksa et al., 1997; Assimaki et al., 2003; see Salloum, 2015 for a synthesis).

	Using suspension log slowness data, Thomson et al. (2009) determined a maximum θ z value of about 5 m; θ x
	was estimated to 50 m to 100 m, assuming an anisotropy factor of 10 to 20. These results illustrate the wide
	scattering of soil statistical parameter values, resulting from the wide variety in geological depositional contexts
	(marine, lake, glacial, alluvial sediments), and from the large uncertainty on the horizontal autocorrelation
	distance quantification.
	As our parametric study is built on the alluvial plain of the river Nahr Beirut (Beirut, Lebanon), the statistical
				Properties
	parameters determined from the extensive near-surface geotechnical (borehole measurements, SPT, laboratory
	Layers measurements) and geophysical (seismic, resistivity) campaigns (Salloum et al. , 2014; Salloum, 2015) are used Model # ρ (Kg/m3) Vp (m/s) θx θy µ Vs (m/s) COV (%)
						(m)	(m)
	hereafter. These experiments showed the presence of interbedded Quaternary layers of pebble, gravel, sand and
						5	2	1
	5 clay overlying marly limestone of Tertiary age, with strong vertical and horizontal variability. The typical soil
						10	2	2
	column is, from the surface to depth: (1) a 7.5 m thick gravel layer with Vs=350 m/s, (2) a softer clay layer of 8 m 5 2 3 15.5 m in thickness with Vs=150 m/s and (3) weathered limestone (Vs>400 m/s). Fitting a theoretical lognormal
	a thickness of probability density function on down-hole tests data led to a Vs coefficient of variation (COV) of 13% and 44% 20 2 6 1600 1500 2 5 220 20 10 1 4
	Soil with				40	5	2 1	7 8
						10
							2	9
	Bedrock	2500	3000	1000	---	---	---	---

Table 1

 1 The finite difference code FLAC 2D is used to propagate seismic waves in the 2D Vs structures[START_REF] Bouckovalas | Numerical evaluation of slope topography effects on seismic ground motion[END_REF][START_REF] Pagliaroli | Studio numerico e sperimentale dei fenomeni di amplificazione sismica locale di rilievi isolati[END_REF][START_REF] Stamatopoulos | Mitigation of the seismic motion near the edge of cliff-type topographies[END_REF][START_REF] Stamatopoulos | Mitigation of the seismic motion near the edge of clifftype topographies using anchors and piles[END_REF]. After discretizing V s for all the models listed in

		Model #5	Model #9
	(m/s)	222	217
	(%)	20.93	38.81
	(m)	11.09	10.85
	(m)	1.91	1.98
	2.4. Wave propagation simulation

, the minimum value of Vs is found to be 50 m/s and the maximum value is 750 m/s. An example of one realization of the Vs structure obtained for model #9 (Table

overall very stiff or very soft 1D Vs profiles leading to significant change in the duration of S-wave main pulse at the surface. This is illustrated in Figure 10d with duration derived from 1D synthetics being short above overall stiff areas while large above overall soft areas (e.g. at receivers located at X = 10 m or X = 45 m). 1).

Spatial correlation of indicators

Spatial correlation of scalar indicators for the 100 Vs realizations using 1D motion synthetics have been computed. Figures 7e,7f, 7g and 7h display the spatial correlation of the resonance frequency (Δf o ), the corresponding amplification (ΔAFf o ), the Arias based intensity (ΔA b I) and the Arias based intensity duration (ΔA b ID), respectively, as a function of ΔX/θ x for all models (Table 1). The variability of all indicators increases linearly until reaching a constant value for ΔX larger than 2θ x , over which correlation between signals does not exist anymore. The larger variability is systematically observed for larger COV values. This is consistent with the control of 1D ground motion synthetics by the locally 1D Vs structures. The 1D f 0 variability for all COV values is larger than that computed for 2D synthetics (compare Figures 7a and7e) in the range ΔX/θ x < 10. This is due to the presence of surface waves in the 2D seismograms, which are little sensitive to short scale variations of locally 1D velocity structure but rather to the homogenized ground structure elastic properties. In contrast, the variability in AFf o is smaller than that observed for 2D synthetics (compare Figures 7b and7f) for ΔX/θ x > 2, whatever the value of COV. This could be related to energetic diffracted surface waves and the double amplification peak close to the resonance frequency, as already mentioned in section 4.1. For ΔX>2θ x , the

[1] Abdallah, I., Malkawi, H., Hassan, W., Abdulla, F. (2000). Uncertainty and reliability analysis applied to slope stability. Structural Safety 22,