Method to simulate wavefields from ambient-noise sources

The shear (SH)-wave transfer function and the horizontal-to vertical (HV) spectral ratio are fast methods to estimate the S-wave velocity profile and thickness of surface layers overlying a bedrock on the basis of resonance frequencies. Understanding the behaviour of the respective frequency spectra is essential to obtain reliable results. In this work, we propose a full-wave numerical method, based on a pseudospectral spatial differentiation, to simulate SH and P-S waves generated by random sources distributed spatially and temporally (ambient noise). The modeling allows us to implement seismic attenuation, surface waves and casual source radiation patterns, based on random values of the angles of the moment tensor at each source location. We focus on the location of the resonance peaks, without considering the amplitudes, since it is known that the HV method cannot predict the amplitudes. First, we analyze Lamb's problem for which an analytical P-S solution exists. The modeling algorithm is verified for a Ricker time history, but the analysis can be performed by using spikes as sources. The experiments based on ambient noise are compared to those of a coherent line source as a reference spectrum (e.g., and earthquake event far away from the receivers). SH-wave resonance frequencies can be identified in the spectra only when the random sources are located below the bedrock. In the case of P-S waves, the SH-wave transfer function is a good approximation to the HV spectrum, mainly when the noise is generated in the bedrock. Finally, we have assumed a square basin and found that coherent (e.g., earthquake-type) sources may yield identiable peaks but ambient noise gives unreliable results.


Introduction
There is nowadays a growing consensus that the most significant source of ambient seismic noise in the Earth is produced by wind-generated ocean gravity waves and their interactions, the ensuing storms and the coupling with the solid earth (Ardhuin et al., 2011).In fact, Ali et al. (2012) observed double-frequency microseisms peaks in the frequency range of 0.15-0.4Hz, generated by the nonlinear interactions of ocean waves with the shoreline along the coasts of the Arabian Sea and the Arabian Gulf.
The use of seismic noise to assess the seismic motion at a site was pioneered in Japan by Kanai in the early 20th century.It is no surprise as the first recordings of strong ground motion were obtained in that country (see e.g.Ishimoto, 1932).Kanai measured microtremors in several sites and noticed the sundry characteristics of noise depending on surface geology.The measurement of microtremor-horizontal-to-verticalspectral ratio or MHVSR (for short HV) has been proposed by Nogoshi and Igarashi (1971) and then applied by Nakamura (1989) to obtain the spectral ratio by dividing the spectrum of the horizontal component by that of the vertical component at the surface, either displacement, particle velocity or acceleration, since the results are equivalent.This spectral ratio is the result of averaging the ratios of particular windows.The most salient feature of HV is the prominence of a peak (usually one) that appears at a frequency that is a reliable indicator of the dominant frequency of the site.This empirical fact led Nakamura (1989) to consider the HV to be the SH-wave transfer function.The proposal was soon followed by controversy, mainly generated by the loose theoretical arguments employed.
The use of HV was supported by the measurements made in Mexico City in very soft ground (Lermo and Chávez García, 1993) They also found that for Rayleigh waves propagating in a layer over a half space, the HV yields the fundamental resonance frequency and the related amplitude is acceptable.This in fact, links HV with ellipticity and, no doubt is related to the special conditions in Mexico City.Other researchers suggested that the HV ratio is closely related to the site SH-wave transfer function (Bonnefoy-Claudet et al., 2008;Oubaiche et al., 2016).
The use of HV became soon a popular tool to assess the dominant frequency and its amplitude has been used an indicator of strong ground motion amplification.This is not widely accepted but it is an easy way to have a proxy.Rigurously, the transfer function is obtained from two receivers at di↵erent depth levels in a vertical seismic well for earthquake sources detected along a borehole.With luck, one can measure the shear waves horizontally polarized (Ohrnberger et al., 2004).Borcherdt (1970) proposed an empirical technique known as the standard spectral ratio (SSR), which requires a suitable reference rock site in the vicinity of the sediment site of interest.In practice it is di cult to find an appropriate reference site, since the response of these rock sites can widely vary.The first numerical test of the HV measurement has been performed by Lachet and Bard (1994), who concluded that the HV peak is independent of the source function and coincides with that of incident shear waves, while the shape depends on the polarization of the fundament Rayleigh wave.In contrast, the HV cannot be used to predict the amplitude of the resonance peaks.
The HV spectral ratio method has been used to estimate the height of sand dunes in the desert.Hanssen and Bussat (2008) investigated a site characterized by recent deposits of sand and sabkha (flat coastal plain with a salt crust) above the bedrock consisting of hard carbonates.Aldahri et al. (2018) used the HV method to determine the resonance frequency and the maximum amplification factor at the Ubhur district, a northern extension of Jeddah in Saudi Arabia.A similar work in the city of Damman has been performed by Al-Malki et al. (2014).In a recent work, this technique has been used to estimate the thickness of glaciers, as well as the basal conditions (Picotti et al., 2017).Lermo and Chávez García (1993) suggested that HV can also be obtained with earthquakes.In fact, in many applications they are used together with ambient noise and active sources (see Alajmi et al., 2016).The HV of deep earthquakes is being used in Japan to identify the velocity structure (Kawase et al., 2011;Nagashima et al., 2014).
In an attempt to model ambient noise and compare techniques to assess siteresponse, within a 2D setting, Coutel and Mora (1998) generate synthetic seismograms with a Chebyshev pseudospectral method for diverse configurations, and test four estimation techniques.Basically, they consider incident SV plane waves (earthquake) and micro tremors (randomly oriented surface sources, i.e., noise) combined with the HV and the HH sediment-to-bedrock ratio.The combinations are: SBSR [SV waves, vx (surface)/vx (bedrock)], SBNR [noise, vx (surface)/vx (bedrock)], HVSR [SV waves, vx (surface)/vz (surface)], and HVNR [noise, vx (surface)/vz (surface)], where v denotes the particle velocity.They conclude that the HVSR, HVNR, and SBNR are unreliable, i.e., yield di↵erent results from the true response as measured by the more reliable SBSR.
Later, Konno and Ohmachi (1998) show that the HV ratio is related to the resonant frequency of the fundamental-mode Rayleigh wave and to the SH-wave transfer function.The ellipticity is invoked again.More numerical tests were performed in the framework of the SESAME project, whose results can be found in Bonnefoy-Claudet et al. (2008).They have considered 1D plane-layered models and sources randomly distributed near the surface, i.e., impulsive and continuous.The HV ratio predicts the resonance frequency of the 1D transfer function corresponding to a vertically incident SH wave.However, it overestimates the site amplification.For high impedance contrasts between the surface layer and the bedrock, the contribution to the peak amplitude comes from the fundamental Rayleigh and Love waves, while for moderate and low contrast the fields are mainly Love and shear body waves.In general, the HV peak corresponds to Rayleigh-wave, Love-wave and/or SH-wave resonances.Van der Baan (2009) explains the resonances obtained from the HV ratio as due to SH and Love waves but in general these depend on several factors, such as the type of source, medium properties, interface geometry, etc.
On the basis of theoretical modeling studies, Albarello and Lunedei (2010) found that surface waves contribute to frequencies larger than the fundamental resonance frequency, whereas body waves contribute to the resonance frequency.Recently, Oubaiche et al. (2016) claimed that the HV peak frequency is better explained by some of the SH-wave transfer function peaks than by the Rayleigh-wave ellipticity, a result that seems to be confirmed by our work, at least for the particular model considered here.
Each specific case has to be analyzed with numerical modeling to obtain a precise interpretation.
A recent development is the discovery of the imaging power of seismic ambient noise.
The works of Shapiro and Campillo (2004) and Shapiro et al. (2005) clearly show the retrieval of surface waves from on a regional scale.Campillo and Paul (2003) obtained the Green's tensor from the correlation of coda waves and Sánchez- Sesma and Campillo (2006) demonstrated the exact relationship between elastodynamic Green's functions and cross-correlations of a uniform set of equipartitioned plane waves.Sánchez-Sesma et al. (2006) extended these results to inclusions in 2D configurations.Perton et al. (2009) proposed the idea of directional energy density related to the imaginary part of the Green function at the source.This led Sánchez-Sesma et al. (2011a,b) to formulate the di↵use field assumption and propose to relate the HV ratio with a square root of the ratio of the sum of imaginary parts of the horizontal Green function over the imaginary parts of the vertical Green function.Some problems with lateral heterogeneity are treated by Matsushima et al. (2014).Rong et al. (2017) compared introduced the empirical transfer function (ETF), defined as the spectral ratio of the records at the surface to the records at a borehole, in order to describe the site amplification of vertical borehole arrays.Shear-wave motion measurements of the ETF is somehow closely related to the S-wave transfer function.These authors show that the amplitude discrepancy is primarily due to two factors, the vertical site response and the HVSR at the bedrock.
The purpose of this work is to explore a full-wave modeling approach to study, by way of numerical experiments, the performance of the SH and HV spectral responses on the basis of the location of the resonance frequencies.We discard any analysis based on amplitudes.We compute the wavefield, based on a modeling method developed by Car-cione (1992Car-cione ( , 2014) ) and verified with the method of generalized reflection/transmission coe cients by Coutel and Mora (1998), who compute HV ratios for several angles of incidence of the SV wave.The computations are based on a Fourier-Chebyshev pseudospectral method, which employs global di↵erential operators in which the field is expanded in terms of Fourier and Chebyshev polynomials along the horizontal and vertical directions, respectively.The proposed algorithm can obtain solutions for general heterogeneous media because the space is discretized on a mesh whose grid points can have varying values of the elastic properties, i.e., the medium can be inhomogeneous.The SH and P-SV(S) formulations are solved in the presence of the free surface, so that the modeling simulates surface waves as well.Wave attenuation is considered, where the quality factors are related to the wave velocities by an empirical relation.

SH-wave transfer function and HV ratio
Let us consider the model shown in Figure 1a, where h is the thickness of the sediment layer, and define by (vx, vz) the in-plane horizontal and vertical particle-velocity components, and by vy the cross-plane component.In plane layered media, these components are decoupled and describe P-S and SH waves, respectively.Let us denote the sediment layer with i = 1 and the bedrock (half space) with i = 2, and define the corresponding complex (Zener) shear-wave velocities as v Si , where where ! is the angular frequency, c S is the high (unrelaxed)-frequency limit velocity, Q S is the minimum quality factor at the frequency f , which is the centre frequency of the relaxation peak, ⌧ = 1/(2⇡f ) and i = p 1. The quantities f , c S and Q S define the media (e.g., Carcione, 2014).When Q S = 1, a S =1 and v S = c S , i.e., the lossless case.The P-wave properties are based on the same equations, replacing S by P .The equations of motion are given in the next section.
In the frequency domain, the HV ratio is In practice, this calculation is rather complex.For real data, it is necessary to perform a statistical analysis of the recorded wavefield in the frequency domain, by computing the amplitude spectra of the three components in a number of selectable time windows.
On the other hand, the body SH-wave transfer function, F , for a viscoelastic sediment layer (sand) of thickness h over a viscoelastic bedrock describes the ratio of the horizontal cross-plane displacements between the top and bottom of the layer due to horizontal harmonic motions of the bedrock.In the literature, it is assumed that A justification of this approximation is given in Lermo and Chávez García (1993).The SH-wave transfer function is (Takahashi and Hirano, 1941;Kramer, 1996), where ⇢ i denotes the mass density.The site transfer function is merely |F |.A rigid bedrock is obtained for ⇢ 2 v S2 ! 1.In this case and in the absence of loss, we have the following resonance frequencies when the cosine vanishes, Infinite amplitude values are obtained at the previously indicated resonance frequencies.

Modeling method
The synthetic seismograms are computed with P-S and SH modeling codes based on an isotropic and viscoelastic stress-strain relation.Sources can be body forces or momenttensor components with random properties.

SH waves
The propagation of viscoelastic SH waves in the (x, z)-plane describes the behaviour of the horizontal cross-plane particle velocity, vy.Euler equation and Hooke law yield the particle-velocity/stress formulation of the SH equation of motion, where denotes stress, e is memory variable, µ = ⇢c 2 S is the shear modulus, fy is a body force, mxy and myz are moment-tensor components, and a dot above a variable denotes time di↵erentiation (Carcione, 2014).
To model surface (Love) waves, free-surface boundary conditions are implemented with the non-periodic Chebyshev operator by using a boundary treatment based on characteristics variables (e.g., Carcione, 1992Carcione, , 2014)).At every time step, the field variable in the free surface are modified as: v ).

P-S wave equation
The time-domain equations for wave propagation in a 2D heterogeneous viscoelastic medium can be found in Carcione and Helle (2004) and Carcione (1992Carcione ( , 2014)).The two-dimensional velocity-stress equations for anelastic propagation in the (x, z)-plane, assigning one relaxation mechanism to dilatational anelastic deformations (l = 1) and one relaxation mechanism to shear anelastic deformations (l = 2), can be expressed by i) Euler-Newton's equations: where xx, zz and xz are the stress components, and fx and fy are external body forces.
ii) Constitutive equations: where e 1 , e 2 and e 3 are memory variables, mxx, mzz and mxz are moment tensor components defining the radiation patterns of the source mechanism: (e.g., Carcione et al., 2015), where M 0 is the moment magnitude, is the dip angle, and k = ⇢(c 2 P c 2 S ) and µ are the unrelaxed (high-frequency) bulk and shear moduli, respectively, where c P is the P-wave velocity.The moment tensor (10) represents a fault whose plane is perpendicular to the (x, z)-plane, where the strike and rake angles are both equal to 90 o .
In nD-space numerical modeling, the dilatational and shear quality factors are functions of the complex bulk and shear moduli, K and µ, respectively.These are , where the subindices denote real and imaginary parts.The quality factor of the P waves is A low-loss relation between these quality factors can be obtained.It is Then, we obtain In 2D space we have and ⌧ (1) 7) holds.Then, we define the unrelaxed velocities, Q P and Q S and obtain the relaxation times with the previous formulae.
The numerical algorithm is based on the Fourier-Chebyshev pseudospectral method for computing the spatial derivatives and a 4th-order Runge-Kutta technique for calculating the wavefield recursively in time.
The Chebyshev method is used along the vertical direction and because it is nonperiodic, it allows us the implementation of the free-surface boundary conditions, to model surface waves (Love and Rayleigh waves).At every time step, the field variable in the free surface are modified as: v At the bottom of the mesh, the implementation of non-reflecting boundary conditions requires: v x zz x ), with = (old) Since the wave equation is linear, we implement time spikes as sources, since seismograms with di↵erent source time histories can be obtained with only one simulation by convolving the source time history with the recorded trace.

Examples
The example considers a sediment layer overlying a sti↵ formation (see Figure 1a), The SH-wave transfer function, |F |, is shown in Figure 2, where h = 100 m and the Zener S-wave relaxation peak is located at a frequency of f = 5 Hz.The black, red and blue lines correspond to a lossy layer over a lossy bedrock, a lossless layer over a lossless bedrock and a lossless layer over a rigid bedrock, respectively.Actually, the blue peaks reach infinite values in the case of a rigid bedrock and a lossless layer, while the location of the peaks are those predicted by equation ( 5).The case of a deformable (non-rigid) bedrock generates finite-amplitude peaks, since energy is lost through transmission in the bedrock, and as can be seen, higher resonances are damped in the lossy case.More details about the e↵ects of attenuation can be found in Carcione et al. (2016).

Lamb's problem
In order to understand the physics of the problem and verify that the full-wave modeling algorithm is correctly approximating the spectrum of the wavefield, we perform simulations for the so-called Lamb's problem, i.e., the propagation of surface (Rayleigh) First, we verify the numerical solution with the analytical solution of Lamb's problem, that is, the response of an elastic (lossless) half-space bounded by a free surface to an impulsive vertical force, fy.In this case there is no bedrock (the properties of the half-space are those of the layer).The analytical solution is obtained by the method of Cagniard-De Hoop (Berg et al., 1994).The equation to solve are ( 8) and ( 9), where the source time history (a Ricker wavelet) is s h (t) = (a 0.5) exp( a), a = [⇡fp(t ts)] 2 , ts = 1.4/fp, with fp = 10 Hz, the source central frequency.The source and the receiver are located at 1.8 m depth and the traveled distance is 700 m. Figure 3 shows the comparison between solutions, where the dots correspond to the numerical solution, whereas Figure 4 shows the analytical (solid line) and numerical (dots) hodograms.
The agreement is very good.The comparison between the anelastic (black line) and elastic (red line) solutions is shown in Figure 5, indicating attenuation and velocity dispersion of the wavefield.
Another test can be performed if one considers the spectra of the traces shown in Figure 3.We consider 4098 (2 12 ) samples and apply a five-point triangular weighted smoothing function to the numerical HV ratio.The comparison with the analytical results is shown in Figure 6, whereas Figure 7 compares the elastic and anelastic spectra.It can be seen that the centroid of the spectrum moves to the low frequencies in the anelastic case.Thus, the numerical code has a good performance and can be used to solve more general problems, such as the layer over half space with the presence of random sources, as shown in Figure 1a.
Before we proceed to attack these problems, we further study the characteristics of the HV spectrum for Lamb's problem.Next, we consider sets of receivers around 700 m o↵set and sum the single HV ratios to obtain an average HV ratio.The average HV spectra for 10 and 30 receivers, compared to a single receiver spectrum located at 700 m from the source, are displayed in Figure 8.The receivers are taken symmetrically around 700 m o↵set.As can be seen, averages can smooth and/or remove possible peaks, because the shape of the curve depends on the source-receiver o↵set as it is illustrated in Figure 9, where the HV spectrum for 100, 300, 700 and 1000 m o↵set is represented.The oscillations increase with o↵set.This e↵ect of averaging has to be tested for each specific geological model.
When solving the velocity-stress formulation with pseudospectral algorithms, the sources can be implemented in one grid point in view of the accuracy of the di↵erential operators.Here, we compute the seismograms by implementing discrete deltas as sources [ (t) in the continuum].The mesh supports a maximum frequency fmax = c min /2dmax, where c min is the minimum velocity and dmax is the maximum grid spacing, so that frequency components greater than fmax will be aliased.However, if the source time-function s h (t) is band-limited with cut-o↵ frequency fmax, those anomalous frequency components will be removed after the time convolution between the seismograms and h(t) or alternatively multiplication with the source spectrum, H(!), in the time-frequency domain.Equivalently, the spectrum of the seismogram obtained with the delta function will be valid till fmax.This is shown in Figure 10, where we compare the HV spectrum obtained from the delta function with that of the bandlimited Ricker function.Then, it is enough to use delta functions as sources, since the spectrum does not depend on the time history, as already found in the literature (e.g., Lachet and Bard, 1994).

The layer-bedrock case
Let us now consider the presence of the layer-bedrock interface.In order to define better the location of the interface and reduce discrepancies with the theoretical transfer function, we consider a 255 ⇥ 217 mesh, a horizontal grid spacing dx = 10 m, a maximum vertical grid spacing of 5.185 m (at the mesh centre) and a vertical extent of 1080 m (including the absorbing boundaries).Fifty grid points of absorbing strips at the sides and bottom of the model yield an e↵ective physical model of 1250 m ⇥ 835 m.The interface is located at a depth of 100 m (the vertical grid points 25 and 26 correspond to the layer and the bedrock, with z = 100 and 105 m, respectively, and the interface is located at grid point 25).
To model the ambient noise, we consider 100 source locations (at single grid points) randomly distributed below the interface (see Figure 1a).Each source location triggers 40 sources (spikes) randomly distributed between 0 and 16.4 s, and each spike has a random amplitude between 0 and 1.Moreover, SH waves are generated with the body force fy, while in the P-S case, the moment tensor (10) is the source, which has random values of the dip angle .In addition, we compute the spectrum corresponding to a horizontal line of sources for reference.Although unrealistic as ambient noise, it may represent an earthquake event far away from the surface layers or basin.
The required time step of the Runge-Kutta algorithm is dt = 1 ms and the solution is propagated 2 14 = 16384 time steps, i.e., ⇡ 16.4 s. Figure 1a shows a scheme of the model, where the dots correspond to the random sources (1000 ⇥ 40 = 40000 spikes, randomly distributed in space and time).Figure 1b shows a snapshot of the SH wavefield at 0.3 s.

SH waves
In the following simulations, we compute the numerical SH-wave transfer function |F |, where we consider the field vy at the first grid point (point 1, surface) and at the last grid point defining the layer (point 25, depth = 100 m), and perform the ratio of the respective frequency spectra.The results confirmed those of Oubaiche et al. (2016) that the HV peak frequency is explained by the SH transfer function.

Resonance frequencies of a basin
The HV method does not work properly in the presence of basins, if the width of the basin is comparable to its thickness.This problem has been analyzed by Zhu and Thambiratnam (2016) and Zhu et al. (2017), for SH and P-S waves, respectively.According to Bard and Bouchon (1985), the SH-body-wave resonance frequencies corresponding to a 2D basin of half-width w and thickness h are where m and n are associated with lateral and vertical interferences.The half-width is defined as the length over which the local sediment thickness is greater than half the maximum thickness.In the case of an infinite horizontal extent of the basin and m=n=0, w ! 1 and f SH = f 0 .For a square basin (h = 2w) and m = 0, the fundamental resonance frequency is giving the peak locations 6.5 Hz, 10.4 Hz, 15.5 Hz, 21 Hz, etc., for f 0 = 2.9 Hz.Zhu and Thambiratnam (2016) find that for high velocity contrasts between the basin and the bedrock, the fundamental frequency is predicted by equation ( 16) (their Table 1), but it is not clear if this equation is equally e↵ective for low contrasts (here the velocity contrast is 2.16) .
We assume a basin with h = 2w = 100 m, where the basin has the properties of the layer above, and consider a horizontal line source at a depth corresponding to the bottom of the basin (vertical grid point 25).In this case, the algorithm requires a time step of 0.5 ms. Figure 16 shows the result of the simulation (lossless case).We consider the field vy at the first grid point (point 1, surface) and at grid point 25 and perform an average of the transfer functions within the basin (10 receivers at the surface and at the bottom).It can be seen that the fundamental peak is not predicted by equation ( 17) (compare the symbols with the first vertical red line).The di↵erences could be due to the fact that Bard and Bouchon (1985) consider a rigid bedrock, so that equation ( 17) cannot be applied to predict resonance frequencies of a basin unless the bedrock is rigid.
In the case of P-S waves, we consider ambient noise below the bedrock, determined by the moment tensor.The results, as those of Figure 16, are shown in Figure 17.
Although there seems to be some apparent agreement with the analytical |F | function regarding the higher modes, the simulation cannot reproduce the fundamental mode, so that nothing can be obtained from this spectrum, unlike that of Figure 14 corresponding to a layer of infinite extent.This confirms the conclusion of Coutel and Mora (1998) that the estimation of site amplification spectra yields unreliable or incorrect results when subsurface basin structure is present, at least for the example of a square basin presented here.Fig. 8 Lamb's problem.Average HV spectra for 10 and 30 receivers compared to a single receiver spectrum located at 700 m from the source (black curve).The receivers are taken symmetrically around 700 m o↵set.

Time (s)
Fig. 9 HV spectrum at 100, 300, 700 and 1000 m o↵sets (analytical solution, the source is a vertical force).The oscillations increase with o↵set.Fig. 16 SH-wave site amplification of a square basin, corresponding to a line source at the bottom of the basin (symbols).Also shown are the peaks of a layer with infinite horizontal extent (blue line) and those of indicated by Bard and Bouchon (1985) for a basin embedded in a rigid medium [equation ( 17)] (vertical red lines).
and body (P-S) waves in the presence of a free surface.The simulations use a 135 ⇥ 81 mesh, with a horizontal grid spacing dx = 20 m, a maximum vertical grid spacing of 19.352 m (at the mesh centre) and a vertical extent of 1412 m (including the absorbing boundaries).Eighteen grid points of absorbing strips at the sides and bottom of the model yield an e↵ective physical model of 2000 m ⇥ 1100 m.
Figure11shows function |F | obtained from the simulations compared to the theoretical function, equation (4) (blue line), where the black dots correspond to random sources and the red dots to a horizontal line source below the interface.The response is the sum of the SH-wave transfer-function spectra of all the receivers and then a five-point triangular weighted smoothing function is applied to this average ratio.The location of the peaks is correctly reproduced and discrepancies may be due to the indetermination about the location of the interface (point sources emitting at all angles "see" this location di↵erently from a vertical propagating line source).The same simulations for sources located in the layer are shown in Figure12(line source at 0.5 m depth).The numerical peaks do not follow those of the SH transfer function, due to the presence of surface (Love) waves, which are not considered in the theoretical transfer function (4).4.2.2 P-S wavesWe conduct the same numerical experiments of those of Figures 11 and 12 in the P-S case to obtain the HV spectrum at the surface.The source of ambient noise is the moment tensor (10) and a line source of horizontal forces, fx (SV waves) at 50 m below the interface.The random sources are deeper than 50 m.Figure13shows the seismograms recorded at the surface, corresponding to the ambient noise.The energy of the wavefield has been attenuated after approximately 10 s, due to anelastic absorption and multiple reverberations within the layer.The results are shown in Figures14 and 15for sources below and above the layer-bedrock interface, respectively.In the latter case, the line source is located at 0.5 m depth and random sources occupy all the vertical extent of the layer, from grid point 2 to grid point 25.Figure14shows that the HV peak resonances can be approximated by the SH-wave resonances, in terms of location of peak frequencies, for this specific model configuration.The correspondence is weaker in Figure15, since the line-source response shows some agreement for the second and third peaks, whereas the random-noise spectrum has a more defined trend, compared to the SH-wave transfer function, mainly regarding the fundamental peak.

Fig. 1
Fig. 1 (a) Model and source-recording configurations.(b) Snapshot of the SH wavefield, where random sources simulating ambient noise (spatially and temporally distributed) are generated.The white line is the layer-bedrock interface.

Fig. 2
Fig.2Analytical SH-wave site transfer function.The black, red and blue lines correspond to a lossy layer over a lossy bedrock, a lossless layer over a lossless bedrock and a lossless layer over a rigid bedrock, respectively.

Fig. 6
Fig.6Normalized spectra of the particle velocity components (a) and HV spectrum (b) for Lamb's problem (elastic, lossless case), where the solid lines correspond to the analytical solution and the dots to the numerical solution.

Fig. 7
Fig. 7 Normalized spectra of the particle velocity components (a) and HV spectrum (b) for Lamb's problem, where the black and red curves correspond to the anelastic and elastic cases, respectively.

Fig. 10
Fig. 10 Lamb's problem.Comparison between the HV spectra obtained from the seismograms generated with the Ricker time history and a delta function.

Fig. 11
Fig. 11 SH-wave transfer function obtained from seismograms generated with a line source and random sources located below the interface (dots), compared to the analytical function (blue line).

Fig. 12
Fig.12SH-wave transfer function obtained from seismograms generated with a line source and random sources located above the interface (dots), compared to the analytical function (blue line).The line source depth is 0.5 m (grid point 2).

Fig. 13
Fig. 13 Seismograms recorded at the surface (center of the model).

Fig. 15
Fig. 15 Same as Figure 14, with the sources located above the layer-bedrock interface.The line source depth is 0.5 m (grid point 2).

Fig. 17 P
Fig.17P-S wave HV ratio of a square basin, corresponding to ambient noise (symbols) below in the bedrock.Also shown are the peaks of a layer with infinite horizontal extent (blue line).