The contribution of multiple scattering to Love wave generation in the secondary microseism

Love to Rayleigh energy ratios in the secondary microseism; 1st row: frequency-wavenumber analysis at the Gräfenberg array (GRF); 2nd row: frequency-slowness analysis using a tiltmeter array in Japan ; 3rd row: ring laser analysis at the Wettzell (WET) ring laser; 4th row: array-derived rotation analysis at Piñon Flat, USA; 5th row: beamforming analysis at several locations in Europe.

Reference	L/R ratio	Location
Friedrich et al. (1998)	0.25	Gräfenberg, Germany
Nishida et al. (2008)	0.6	Hi-net, Japan
Tanimoto et al. (2015)	1.25	Wettzell, Germany
Tanimoto et al. (2016)	0.4–0.5	Piñon Flat, USA
Juretzek & Hadziioannou (2016)	0.6–1.2	Europe

Reference	L/R ratio	Location
Friedrich et al. (1998)	0.25	Gräfenberg, Germany
Nishida et al. (2008)	0.6	Hi-net, Japan
Tanimoto et al. (2015)	1.25	Wettzell, Germany
Tanimoto et al. (2016)	0.4–0.5	Piñon Flat, USA
Juretzek & Hadziioannou (2016)	0.6–1.2	Europe

Table 1.

Love to Rayleigh energy ratios in the secondary microseism; 1st row: frequency-wavenumber analysis at the Gräfenberg array (GRF); 2nd row: frequency-slowness analysis using a tiltmeter array in Japan ; 3rd row: ring laser analysis at the Wettzell (WET) ring laser; 4th row: array-derived rotation analysis at Piñon Flat, USA; 5th row: beamforming analysis at several locations in Europe.

Reference	L/R ratio	Location
Friedrich et al. (1998)	0.25	Gräfenberg, Germany
Nishida et al. (2008)	0.6	Hi-net, Japan
Tanimoto et al. (2015)	1.25	Wettzell, Germany
Tanimoto et al. (2016)	0.4–0.5	Piñon Flat, USA
Juretzek & Hadziioannou (2016)	0.6–1.2	Europe

Reference	L/R ratio	Location
Friedrich et al. (1998)	0.25	Gräfenberg, Germany
Nishida et al. (2008)	0.6	Hi-net, Japan
Tanimoto et al. (2015)	1.25	Wettzell, Germany
Tanimoto et al. (2016)	0.4–0.5	Piñon Flat, USA
Juretzek & Hadziioannou (2016)	0.6–1.2	Europe

Apart from Friedrich et al. (1998), all studies shown in Table 1 convert the surface amplitude to kinetic energy by integrating over fundamental mode surface wave eigenfunctions. Friedrich et al. (1998) only uses the squared velocity to obtain a measure of energy.

Longuet–Higgins’ theory, which is shown to be responsible for the Rayleigh wave excitation, can not explain the amount of Love wave energy in the secondary microseism. Here, we propose a mechanism which does not happen at the source, but rather on the propagation path between source and observation location.

The seismic wavefield undergoes strong scattering and conversions through the heterogeneous structure of the Earth (Levander & Holliger 1992). Scattering is more important for the high frequency wavefield, but also plays a role in the microseismic frequency band (Lee et al. 2003). The energy partition of the different wave types is not a quantity only determined by the source but is also strongly influenced by the heterogeneity of the medium (Aki et al. 1976). Therefore, the observed amount of Rayleigh and Love wave energy measured in the secondary microseism is not necessary a source effect, but can also change during the propagation of the wavefield through scattering and conversion (Rind & Down 1979).

In this study, we aim to give a better insight into the contribution of this propagation effect. We try to quantify the change of the Love and Rayleigh wave energy which results from internal scattering and conversions of the propagating microseismic wavefield.

Only a few studies exist for surface wave scattering, mostly applied to Rayleigh to body wave conversion and single scattering. Some studies, such as Kennett (1984), consider conversions towards Love waves as well. Kennett (1984) calculated the reflection and transmission matrices of a surface wavefield that passes a finite heterogeneous area (e.g. sedimentary basin). He showed that the wavefield in the perturbed area can be calculated by the coupling of the surface wave modes of the background field (unperturbed field). This approach was extended to 3-D structures in Kennett (1998). They found that progressively more energy was coupled into the transverse component, implying that a significant portion of the wavefield ends up as Love waves. Snieder (1986) derived scattering coefficients for surface waves under the restriction of the Born- and far field approximation. By using a mountain root as an example, he calculated the 2-D surface wave scattering patterns for periods between 20 and 200 s. Besides the coupling of the surface wave types, the derivations are also usable for the coupling of different surface modes. Although large scale scatterers such as mountain roots or sea mounts could affect the secondary microseismic wavefield, the mainly single scattering approach in this study would not suffice to explain the relatively high Love wave content we observe. Maeda et al. (2008) used the results of Snieder (1986), and calculated the scattering coefficient for Rayleigh to body waves and vice versa. They used these to synthesize the coda envelope on the surface of a half-space. Love waves are not considered in this study. Maupin (2001) extended previous scattering theory to an anisotropic scatterer. Using a fundamental mode as an incident wave, she showed that less than 5% of the scattered energy is stored in higher modes. Some studies have addressed topographical scattering of body and surface waves (e.g. Bannister et al. 1990; O’Brien & Bean 2009; Kumagai et al. 2011; Takemura et al. 2015), but since we focus on internal scattering in this study, we will not discuss these studies in detail here.

We address two cases in this study. In the first part, we try to get better understanding of the Rayleigh and Love wave generation in random media. To this end, we perform a parameter study where we quantify the Love-to-Rayleigh (L/R) ratio depending on different scattering parameters. In a second part, we investigate the role of scattering in the secondary microseism frequency band, using realistic scattering properties of the Earth. Since there is no analytical solution for this type of problem, we use numerical simulations to calculate the surface wavefield in a random medium, and compare our results to observed L/R ratios. For the numerical simulations we use the finite difference solver Sofi3D (Bohlen 2002), which was successfully tested for multiple scattering in random media (Bohlen 2002; Huang et al. 2006). We neglect intrinsic attenuation since we seek for the maximum amount of scattering (Bohlen 2002).

This paper structured as follows: in Section 2 we discuss the theoretical background of scattering of seismic waves and the rotation measurements that we use to separate the Rayleigh from the Love waves. The model setup is detailed in Section 3. There, we explain how we determine our numerical parameters and how we set up the receiver geometry to calculate the vertical and horizontal rotation. The data processing in Section 4 explains how we calculate the L/R ratio. For that, we use the total L/R ratio in dependence on the source–receiver distance. In Section 5, we show the results for the parameter study and the realistic model. Finally we discuss our results and give a short outlook in Section 7.

2 THEORY

2.1 Scattering

As a result of processes ranging from mineralogical alteration to tectonic deformations, the lithosphere is a quite complex medium with heterogeneties from millimetre to kilometre scales (Levander & Holliger 1992). Propagating seismic waves are affected by the heterogeneous medium and undergo strong scattering and mode conversions. The scattering properties of the crust can be mathematically described by a deterministic part (e.g. mean velocity) and a random fluctuation term (Wu et al. 1994; Holliger 1996)

\begin{eqnarray*} q(x)=q_0(1+\xi (x)). \end{eqnarray*}

(1)

Here q₀ is the mean value (e.g. density or velocity) and ξ(x) is the fractional fluctuation characterized by a Gaussian distribution. To calculate the fractional fluctuation, we can use the 3-D inverse Fourier transform of its squared power spectral density function (PSDF). Depending on the problem to solve, different types of PSDF’s can be used. Studies on borehole data (e.g. Holliger 1997; Dolan et al. 1998) and numerical simulations (Frankel & Clayton 1986) confirm that the von Karman function is a good choice to model the Earth’s crust. Its 3-D PSDF is given by (Sato et al. 2012):

\begin{eqnarray*} P({\bf m})=\frac{8\pi ^{3/2}\Gamma (\kappa +3/2)\sigma ^2a^3}{\Gamma (\kappa )(1+{\bf a}^2{\bf m}^2)^{\kappa +3/2}} . \end{eqnarray*}

(2)

a is referred to as the correlation length and gives the corner wave number of the PSDF, while the Hurst exponent (κ) gives the decay of the PSDF after the corner wave number. σ is the fluctuation strength, which determines the deviation from the mean elastic parameter (e.g. v_p). We choose m as the wave number of the medium to differ from the wave number k of the wavefield. Γ is the gamma function.

For am > 1, the PSDF exhibits self similar behaviour, which means that the fluctuations have the same characteristics independent of the wavenumber m. a can be a vector, where different entries a_i(i=x,y,z) lead to an anisotropic medium with different scattering properties depending on the propagation direction of the seismic field.

The random medium can be classified according to different scattering regimes (Pyrak-Nolte 2002). To distinguish between these regimes we can use the non-dimensional parameter ka. If ka ∼ 1, the wavelength of the seismic wavefield is similar to the size of the scatterers. The scattered wavefield is radiated with large angles relative to the incident direction. Therefore, we call it the large angle scattering regime. The interaction between the wavefield and the scatterer is maximal and we observe a large coda (Pyrak-Nolte 2002; Sato et al. 2012). If ka > 1, the scatterers are larger than a wavelength. The wavefield is scattered in forward direction and diffraction becomes more important. The third case is the Rayleigh scattering regime, where ka < 1. The fluctuations are weak and scattering can be described by the Born approximation.

In the first part of this study, we concentrate mainly on the large angle scattering case. We aim for a maximum amount of scattering, in order to obtain information about the L/R ratio for relatively short propagation paths. In reality, it is more likely that the secondary microseismic wavefield is mainly situated in the Rayleigh scattering regime. We discuss this case in the second part of this study.

The scattering coefficient is the scattering power per unit volume and is reciprocal to the scattering mean free path (SMFP) (Sato et al. 2012). The SMFP gives distance after which the wavefield is statistically scattered the first time. If the SMFP is larger than the propagation path, the wavefield is in the single scattering domain, whereas for a propagation of several SMFP the wavefield enters the diffusive regime (Sato et al. 2012). In Section 4, we show a method to determine the SMFP, which in turn helps us quantify the scattering strength of the medium.

2.2 Rotation

In this study, we exploit the rotational components of the wavefield to separate the Rayleigh from the Love waves. Rayleigh waves are elliptically polarized in the vertical plane, and have a rotational component around the transverse axis. Love waves, which result from the interaction of SH waves, are transversely polarized and have a rotational component around the vertical axis. Given a flat, layered half-space, we can use the measurements of rotational motion around vertical and horizontal axes as a tool to separate the two wave types.

To calculate the rotations, we use the array-derived rotation (ADR) approach presented in Spudich et al. (1995). The displacement is recorded by a three-component micro-array as shown in Fig. 2. The distance ∂x between two receivers is given by the distance of two gridpoints in our model. The different micro-arrays are separated by a larger distance, which is used later in equations (6) and (7). First, we calculate the displacement gradients |$\frac{\partial u_i}{\partial x_j}$| where u_i is the displacement in i direction, and where i, j are 3-D Cartesian components. The rotation for a micro-array can then be calculated using:

\begin{eqnarray*} \Omega (t) = \frac{1}{2}\left( \frac{\partial u_j(t)}{\partial x_j}-\frac{\partial u_i(t)}{\partial x_j} \right), \end{eqnarray*}

(3)

In the ADR approach, the aperture of the array limits the available wavelengths. A limitation is that h < λ/4, where h is the array aperture and λ is the minimum wavelength of the wavefield (Spudich & Fletcher 2008).

3 MODEL SETUP

We solve the 3-D elastic wave equation in a semi-infinite medium with the physical parameters (velocity, density) described by a von Karman correlation function (see Fig. 1). As this study focuses on conversions of Rayleigh to Love waves due to internal scattering, we use a model without surface topography or oblique internal layers. An efficient way to solve the wave equation in such cases is through finite difference solvers. For all calculations, we use the SOFI3D solver (Bohlen 2002), which has the advantage of having a parallel MPI environment. The solver has already been successfully tested for seismic wave propagation in random media (Bohlen 2002).

Figure 1.

Two realizations of a random medium with different correlation lengths; top: correlation length a = 50 m; bottom: correlation length a = 500 m; the different colours show the random distribution of the density [kg m^–3] given by the von Karman PSDF (see eq. 2). Both plots show a single layer (SL) layering type, as detailed in Fig. 3.

The models consist of an absorbing frame inside the medium and a free surface at the top. To mimic the absorbing frame we use the perfectly matched boundary condition (PML) which is implemented in the solver. The solver documentation suggests using between 10 and 20 gridpoints to avoid edge reflections. In tests with and without random media, the choice of 20 gridpoints showed negligible boundary reflections, even in strongly scattering media.

Additionally, we choose a 8th order Holberg operator for the spatial discretization. According to Holberg (1987) this gives a minimum amount of 3.19 gridpoints per wavelength. However, the strong scattering of the wavefield requires more gridpoints than in normal propagation conditions. Therefore, we choose a relatively dense grid to prevent grid dispersion: in order to have 20 gridpoints per central Rayleigh wavelength, the grid size is chosen at 10 m. In a second part, we consider a model which approaches realistic conditions. The exact model setup of the realistic model is shown and discussed in Section 5.4.

All models considered in this study have a similar source type and receiver arrangement. To mimic the 'Longuet-Higgins'-type source mechanism, we apply a vertical point force at the surface. While this is a simple interpretation of the source mechanism, it ensures that the source generates no transversal polarized waves. The source time function is given by a Ricker wavelet with a central frequency of 15 Hz. The source and receiver geometry can be seen in Fig. 2(a). The source is located at the centre of the surface and the receivers are placed at the free surface, radially in lines around the source. We place 32 micro-arrays of receivers, which consist of a tetrahedral shape with three receivers at the surface and one at depth. These micro-arrays, with inter-receiver distance of 20 m, are shown in Fig. 2(b). The radial 'star'-like shape of the receivers is chosen to allow the calculation of the L/R ratio as a function of source–receiver distance. Moreover, since the medium is isotropically scattering, we can average over all directions to calculate the mean wavefield, instead of simulating multiple realizations of the random medium.

Figure 2.

(a): Source and receiver arrangement at the free surface, as seen from above; red circle: source location; blue triangles: receivers; (b): zoom of the top plot (green square), illustrating the micro-array setup (red squares).

The two sections of the study use two different model setups. In a first part, we perform a parameter study where we consider the effect of different scattering conditions on the Love-to-Rayleigh ratio. We use a smaller model with a size of 20 x 20 x 10 km, and a micro-array distance of 300 m.

We determine the L/R ratio for varying values of the fluctuation strength σ and the correlation length a (see Table 2). While we vary one parameter, we fix the others. For example, we vary σ while we use a constant correlation length of 200 m and a multilayer plus gradient (MLGR) layering type. This way, we can study the L/R ratios for each scattering parameter individually. Table 2 lists all combinations of scattering properties and layering types that we use for this part of the study. The scattering properties are the same for the velocities and the density, which ensures the best conditions for the wave type conversion (Bannister et al. 1990). To create the von Karman media we use the same method as Frankel & Clayton (1986).

Table 2.

Statistical properties of the random media for the different model implementations. All models use a Hurst exponent κ=0.3. a is the correlation length; σ is the fluctuation strength; the acronyms for the layering types are explained in Fig. 3; λ_max/a is the ratio between the central wavelength and the correlation length a; SMFP is the scattering mean free path around the central frequency of 15 Hz; L/R is the maximum Love to Rayleigh ratio reached for the corresponding model.

σ [%]	a [m]	Layering	λ_max/a	SMFP [m]	L/R
15	50	MLGR	∼ 4	3855	0.21
5	200	MLGR	∼ 1	9536	0.06
10	200	MLGR	∼ 1	3266	0.24
15	200	MLGR	∼ 1	1741	0.44
20	200	MLGR	∼ 1	1196	0.47
30	200	MLGR	∼ 1	613	0.49
15	500	MLGR	∼ 0.4	1896	0.37
15	200	SL	∼ 1	1948	0.35
15	200	ML	∼ 1	1946	0.36
15	200	GR	∼ 1	1940	0.34
5	2000	ML	∼ 8	1.222 · 10⁶	0.11
10	2000	ML	∼ 8	6.2 · 10⁵	0.17

σ [%]	a [m]	Layering	λ_max/a	SMFP [m]	L/R
15	50	MLGR	∼ 4	3855	0.21
5	200	MLGR	∼ 1	9536	0.06
10	200	MLGR	∼ 1	3266	0.24
15	200	MLGR	∼ 1	1741	0.44
20	200	MLGR	∼ 1	1196	0.47
30	200	MLGR	∼ 1	613	0.49
15	500	MLGR	∼ 0.4	1896	0.37
15	200	SL	∼ 1	1948	0.35
15	200	ML	∼ 1	1946	0.36
15	200	GR	∼ 1	1940	0.34
5	2000	ML	∼ 8	1.222 · 10⁶	0.11
10	2000	ML	∼ 8	6.2 · 10⁵	0.17

Table 2.

Statistical properties of the random media for the different model implementations. All models use a Hurst exponent κ=0.3. a is the correlation length; σ is the fluctuation strength; the acronyms for the layering types are explained in Fig. 3; λ_max/a is the ratio between the central wavelength and the correlation length a; SMFP is the scattering mean free path around the central frequency of 15 Hz; L/R is the maximum Love to Rayleigh ratio reached for the corresponding model.

σ [%]	a [m]	Layering	λ_max/a	SMFP [m]	L/R
15	50	MLGR	∼ 4	3855	0.21
5	200	MLGR	∼ 1	9536	0.06
10	200	MLGR	∼ 1	3266	0.24
15	200	MLGR	∼ 1	1741	0.44
20	200	MLGR	∼ 1	1196	0.47
30	200	MLGR	∼ 1	613	0.49
15	500	MLGR	∼ 0.4	1896	0.37
15	200	SL	∼ 1	1948	0.35
15	200	ML	∼ 1	1946	0.36
15	200	GR	∼ 1	1940	0.34
5	2000	ML	∼ 8	1.222 · 10⁶	0.11
10	2000	ML	∼ 8	6.2 · 10⁵	0.17

σ [%]	a [m]	Layering	λ_max/a	SMFP [m]	L/R
15	50	MLGR	∼ 4	3855	0.21
5	200	MLGR	∼ 1	9536	0.06
10	200	MLGR	∼ 1	3266	0.24
15	200	MLGR	∼ 1	1741	0.44
20	200	MLGR	∼ 1	1196	0.47
30	200	MLGR	∼ 1	613	0.49
15	500	MLGR	∼ 0.4	1896	0.37
15	200	SL	∼ 1	1948	0.35
15	200	ML	∼ 1	1946	0.36
15	200	GR	∼ 1	1940	0.34
5	2000	ML	∼ 8	1.222 · 10⁶	0.11
10	2000	ML	∼ 8	6.2 · 10⁵	0.17

Following Imperatori & Mai (2012), we choose a Hurst exponent of κ=0.3 for all models. We tested the effect of Hurst exponents in the range of 0.1–0.5 (Takemura et al. 2015), but the influence on the results was negligible compared to that of the other discussed parameters (σ, a, layering).

The central source frequency of 15 Hz leads to central wavelength of 200 m. For a correlation length of 500 m it means that we are in the forward scattering regime. For the correlation length of 50 and 200 m we are in the large angle scattering regime.

Underlying each random medium is a layered velocity model. Several different layering types, such as single, multiple, gradient layers, are tested, as illustrated in Fig. 3(a). Here, the depth dependence of the mean velocity <v_p> and the mean density <ρ> are shown for the different layering types. <v_s> is determined with a Poisson ratio |$v_p/\sqrt{3}$|⁠. The single layer (SL) type contains only one layer, while the multiple layering (ML) type contains four layers. We also test a gradient layer (GR), where the physical properties of the model increase gradually as well as a a multiple layer plus gradient (MLGR) type, where we combine both.

Figure 3.

Mean velocity profiles 〈v_p〉 and mean density profiles 〈ρ〉 for the different models that we use in the parameter study (two left-hand panels), and for the real Earth-like crustal model (two right-hand panels).

In the second part of the study, we consider scattering conditions which are closer to typical crustal scattering conditions in the secondary microseismic frequency range. Several observations show correlation lengths below 1 km (Holliger & Levander 1992; Sens-Schönfelder et al. 2009), but there are also studies with much larger correlation lengths (e.g. Hillers et al. 2013).

We choose a correlation length of 2 km. In combination with a scattering strength σ of 5 and 10% and a Hurst exponent of 0.3, the modeled SMFP are comparable to those observed in the secondary microseism frequency band. For a σ of 5% this leads to a SMFP of 1222 km. Lee et al. (2003) calculate a SMFP of 1605 km for a frequency of 0.25 Hz. For a σ of 10% we obtain a SMFP of 620 km, similar to the SMFP of 690 km observed by Sens-Schönfelder et al. (2009) for a frequency of around 0.2 Hz.

The Ricker source time function used in this part of the study has a central period of 5 s, to mimic the central period of the secondary microseism. This leads to a central wavelength of 16 km and a λ_C/a ratio of 8. This ratio is larger than in the smaller models in the first part, and leads to less scattering events, which in turn explains the significantly larger SMFP (see Table 2). We are still in the large angle scattering regime. Since the central wavelength is quite large, we now use a model size of 1000 x 1000 x 500 km to allow for >50 wavelengths propagation distance. We use a coarser grid, with an element size of 500 m.

The physical properties of this model can be seen in Fig. 3(b). Again, <v_s> is determined with a Poisson ratio of |$v_p/\sqrt{3}$|⁠. The receiver arrangement is similar to Fig. 2(a), but with the micro-arrays placed at a distance of 15 km. The receiver–receiver distance inside each micro-array is 1000 m.

4 DATA PROCESSING

We solve the elastic wave equation for several random model implementations, as listed in Table 2. For that we use a Ricker source time function with a central frequency of 15 Hz and a maximal frequency of 30 Hz. In order to separate the wavefield into Love and Rayleigh waves, we calculate the rotation around the vertical axis (corresponding to Love waves) and around the transverse horizontal axis (corresponding to Rayleigh waves). This is done for each micro-array, using the approach described in eq. (3). Fig. 4 shows the seismograms for the vertical (red colour) and transverse horizontal (blue colour) rotations for two different realizations of a random model. Fig. 4(a) shows the seismograms for a correlation length of the same order as the central wavelength (a = 200 m) and a relatively weak fluctuation strength (⁠|$\sigma =5\%$|⁠). Fig. 4(b) shows seismograms for the same correlation length (a = 200 m), but for a much stronger fluctuation strength (⁠|$\sigma =20\%$|⁠).

Figure 4.

Seismograms for the rotations around the transverse axis Ω_T (Rayleigh waves, blue) and the vertical axis Ω_V (Love waves, red) at different distances from the source for two models with a varying σ of 10% (top panel) and 20% (bottom panel) (see Table 2). The correlation length a is 200m and the layering type is MLGR for both models.

For a weakly heterogeneous medium with a fluctuation strength of |$\sigma =5\%$|⁠, the Rayleigh wave first arrival is well pronounced for increasing distance and the coda is weak. This confirms the small amount of scattering that we expect for such a medium. For a σ of 20%, the direct Rayleigh wave has disappeared after several wavelengths propagation distance and the coda is well developed for both Rayleigh and Love waves. Another interesting observation is that, for distances larger than approximately 3 km, the coda of the Love wave starts earlier than that of the Rayleigh waves. The Love wavefield propagates faster and overtakes the Rayleigh wavefield after a certain distance.

According to Sato et al. (2012), the energy density decrease of a surface wavefield in a random isotropic medium is proportional to:

\begin{eqnarray*} E(t,r)\propto 1/r*\exp (r/l), \end{eqnarray*}

(4)

where l is the SMFP and r is the source–receiver distance. In eq. (4), we distinguish between the geometrical spreading of surface waves, given by the factor 1/r, and scattering attenuation, given by the exponential factor. For simplicity we neglect intrinsic attenuation. The energy density of the field can be calculated by:

\begin{eqnarray*} E(t,r)=\frac{\rho _0[\dot{u}(t,r)^2+\mathcal {H}^2(\dot{u}(t,r))]}{2}, \end{eqnarray*}

(5)

where ρ₀ is the mass density of the medium and |$\mathcal {H}$| is the Hilbert transform. Fig. 5 shows the energy density envelope for two different realizations of the medium with fluctuation strengths σ of 10 and 20%. Similar to Fig. 4 we can see how the energy is redistributed from the direct wavefield to the scattered wavefield in the strongly scattering medium.

$Total mean squared (MS) envelope of the vertical velocity component $\dot{u}_z$ at source–receiver distance of 1.2 km for a fluctuation strength of 10% (top panel) and 20% (bottom panel). The correlation length a is 200 m and the layering type is MLGR for both models.$

Figure 5.

Total mean squared (MS) envelope of the vertical velocity component |$\dot{u}_z$| at source–receiver distance of 1.2 km for a fluctuation strength of 10% (top panel) and 20% (bottom panel). The correlation length a is 200 m and the layering type is MLGR for both models.

Eq. (4) can be used to calculate the SMFP. For that we take the natural logarithm on both sides of the equation and apply a linear regression of the log energy density versus the distance r. The SMFP is then given by the inverse of the slope. Fig. 6 shows the linear regression for four different models. Small σ corresponds to weak scattering, and a high SMFP (Fig. 6, cyan line), while high σ corresponds to strong interaction between wavefield and scatterer and results in a small SMFP (Fig. 6, blue line). In our different media, we have SMFP ranging from 0.5 to 10 km, which ensures a wide range of scattering properties for the parameter study.

Figure 6.

Linear regression of the normalized energy density peak of the coherent Rayleigh wave versus the distance. As an example we see the linear regressions of four different σ’s (for MLGR). The dots show the peak coherent Rayleigh wave for different distances and angles from the source for σ = 20%. The central frequencies for all these calculations is 15 Hz. The correlation length a is 200 m and the layering type is MLGR for all models.

Before we calculate the L/R ratio, we first have to verify that the seismograms contain surface wave energy, and that both Rayleigh and Love waves are present in the wavefield. As already mentioned, we use the transverse and vertical rotation rates to separate the Rayleigh and Love waves. To identify the individual surface wave modes, we use a delay and sum method given by the following equation:

\begin{eqnarray*} \Phi ^{V/T}(\omega ,v)=\frac{1}{N}\sum _{l=1}^N \hat{\Omega} ^{V/T}_l(\omega )e^{i \omega r_l/v}. \end{eqnarray*}

(6)

|$\hat{\Omega }_l$| are the Fourier transformed rotation rate around the vertical (V) and transverse (T) axis of the lth microarray (eq. 3) and N is the number of micro-arrays in the radial direction. The exponent applies a phase shift to the trace Ω_l. For certain phase shifts, the traces are in phase and add up constructively, which maximizes Φ. Repeating this procedure for different frequencies and velocities, we can find the phase-velocity frequency representation of the surface wave modes. The parameters inside the exponent are the angular frequency ω, the distance r_l between the first receiver and current receiver l, and the velocity v. The velocity range can be chosen freely, but should contain the velocities of the modes that we intend to study. In our case, we choose a minimum velocity of 2000 m s^–1, which is a little bit lower than the fundamental Rayleigh mode, and a maximum of 5000 m s^–1 which is enough to see the first four higher modes of both Love and Rayleigh waves.

The centre and the right plot of Fig. 7 shows the result of such a calculation for the transverse and the vertical rotation rates, for a model without layering. The mean physical properties of this medium are given by v_p = 5000 m s^–1, v_s = 3500 m s^–1 and ρ = 2900 kg m^–³. In the centre plot, we see a good match between the analytical fundamental Rayleigh mode, indicated by the dotted line, and the strongly coherent energy. Since no layering is present, the velocity is constant for all frequencies. The coherent energy in the right-hand plot corresponds to the SH body waves. As expected, we observe no Love wave energy in this case. Note that the Love-to-Rayleigh ratio (left) is negligibly small.

$Model without layers and with the scattering parameters $\sigma =15\%$ and a = 200 m; left-hand panel: L/R ratio as a function of the source–receiver distance. Note that the vertical scale is much smaller than in subsequent plots. The black line shows the mean over the 14 angular directions (see Fig. 2) and the yellow area shows one standard deviation; centre panel: dispersion curve for the rotations around the transverse axis ΦT (perpendicular to the source–receiver direction). The black dotted line shows the theoretical Rayleigh wave fundamental mode; right-hand panel: dispersion curve for the rotation around the vertical axis ΦV. The colour scale corresponds to rotation rate.$

Figure 7.

Model without layers and with the scattering parameters |$\sigma =15\%$| and a = 200 m; left-hand panel: L/R ratio as a function of the source–receiver distance. Note that the vertical scale is much smaller than in subsequent plots. The black line shows the mean over the 14 angular directions (see Fig. 2) and the yellow area shows one standard deviation; centre panel: dispersion curve for the rotations around the transverse axis Φ_T (perpendicular to the source–receiver direction). The black dotted line shows the theoretical Rayleigh wave fundamental mode; right-hand panel: dispersion curve for the rotation around the vertical axis Φ_V. The colour scale corresponds to rotation rate.

The surface wave eigenfunctions are strongly perturbed due to the heterogeneous medium, preventing us from integrating over the depth sensitivity to obtain the wave energy. As a solution, we calculate the signal energy instead. We calculate the L/R ratios by averaging over all micro-arrays which share the same radial distance l with the following equation:

\begin{eqnarray*} \mu _l^{L/R}=\frac{1}{M}\sum _{k=1}^{M}\frac{\int ^{t_e}_{0} \Omega ^V_{kl}(t)^2dt}{\int ^{t_e}_{0} [\Omega ^T_{kl}(t)^2+\Omega ^R_{kl}(t)^2]dt}. \qquad \end{eqnarray*}

(7)

where M is the number of micro-arrays in angular direction (k), and l is the distance between the source and current micro-array (lth micro-array in radial direction). The superscripts in Ω stand for the vertical (V), the transverse (T) and the radial (R) components. We integrate over the whole time series (t_e is the end of the time-series), and do not separate the direct and the scattered wavefield, since we are mainly interested in the change in energy distribution from Rayleigh waves towards the transversal polarized wave types.

The L/R ratios will be slightly different for different realizations of the random medium, following a normal distribution. The exact L/R ratio (expectation value) can be approximated statistically by taking the mean over many simulations, calculating the wavefield for the same background model with a different realization of the medium each time. Since this would be computationally very expensive we follow another approach. Since our random velocity model has an isotropic autocorrelation function, the wavefield scattering only depends on the source–receiver distance and not on the propagation direction. We can make use of this simplification and calculate the wavefield only once in a single random medium. The wavefield is then recorded along 14 receiver lines, each emanating at a different angle from the central source location (see Fig. 2). Under the assumption that each receiver line is roughly equivalent to a different realization of a random velocity model, the L/R ratio is measured as a function of distance along each of these lines, subsequently the mean and standard deviation are calculated.

The left plot in Fig. 7 shows the L/R ratio calculated with eq. (7). The L/R is plotted in dependence of the distance r. The black line shows the mean L/R ratio over different receiver lines, the yellow filled area shows the standard deviation. The L/R ratio stays consistently low over the whole distance, showing that no Love waves are able to develop in the unlayered medium.

5 RESULTS

5.1 L/R ratio as a function of the fluctuation strength σ

In order to quantify the L/R ratio as a function of the fluctuation strength σ, we calculate the wavefields for a random medium with a constant correlation length of 200 m, such that a ∼ λ. The fluctuation strength σ is varied between 5 and 30%, which corresponds to a maximum SMFP of ∼9500–600 m. (see Table 2). For each test we use the multiple layers with gradient (MLGR) layering type (see Fig. 2a). We quantify the L/R ratio using eq. (6) to extract the surface wave modes from the signal, and eq. (7) to calculate the L/R ratios as a function of the distance.

The top row of Fig. 8 shows the dispersion curves for a σ of 10%, the bottom row for |$\sigma$|= 20%. The left-hand column shows the dispersion curves of the transversal curl axis, while the right-hand column shows the dispersion curve of the vertical curl axis. The analytical Rayleigh and Love wave fundamental modes are plotted as black dotted lines. In all dispersion plots, the maximum centred around the theoretical dispersion curves indicates the presence of Rayleigh and Love waves, respectively. We note that the dispersion curves become less coherent with increasing fluctuation strength.

Figure 8.

Dispersion curves for a σ of 10% (top row) and a σ of 20% (bottom row). The correlation length a is 200 m for both models. Both models use a MLGR layering type; left-hand panels: dispersion curve for the rotations around the transverse axis Φ_T (perpendicular to the source–receiver direction). The black dotted line shows the theoretical Rayleigh wave fundamental mode; right-hand panels: dispersion curve for the rotation around the vertical axis Φ_V. The dotted line indicates the fundamental Love wave mode. The colour scale corresponds to rotation rate. Each Φ_V dispersion curve is normalized by the maximum of the corresponding Φ_T dispersion curve.

Each dispersion plot is normalized by the maximum of the transverse dispersion curve, which turns the maximum of the colour-scale as a proxy for the L/R ratio. Similar to the L/R ratio calculation of Fig. 9, we observe an increasing L/R with increasing fluctuation strength. However, the ratio between the dispersion plot maxima is overall smaller than the L/R calculated with the overall signal energy. This could have two origins. A possible cause is that the rotations around the vertical axis include S-phase energy of the reflected and scattered waves, which would increase the total L/R ratio in the left-hand column. Another possibility is the presence of higher modes, which are not taken into account in the approximate L/R ratio determined from the dispersion curves.

Figure 9.

L/R ratios as a function of distance, normalized by the SMFP (see Table 2), for varying fluctuation strength σ (increasing from left to right). The correlation length a is 200 m and the layering type is MLGR for all models. The black line shows the mean over the 14 angular directions (see Fig. 2) and the yellow area indicates one standard deviation; The red dotted line indicates the stabilized L/R ratio; The scattering properties of the random media are described above the plots.

Higher modes are not visible, independently from the choice of σ. Some coherent energy is visible above 3000 m s^–1. This is an artefact of the method we used, and is a result of the periodicity of the phase.

Fig. 9 shows the L/R ratios as a function of the normalized distance r_s, which is the distance divided by the SMFP (see Table 2). This way, we can describe the L/R ratio depending on the number of scattering events. The fluctuation strength |$\sigma$| increases from left to the right.

Two major things are visible: first, the L/R ratio increases with increasing source–receiver distance; second, the maximum value of L/R ratio is strongly dependent on σ. The increase of Love to Rayleigh ratio with propagation distance points towards gradual conversions of Rayleigh wave energy into Love waves.

All four models show a stabilization of the L/R ratio after several SMFP. For a σ of 15 , 20 and 30% the stabilization is reached after three to four SMFP, and the stabilization ratio has a mean around 0.4 (see red dotted line in Fig. 9). For a σ of 10%, the total propagation distance, limited by the size of the model, only covers 2.5 SMFP, which seems to be sufficient to reach the onset of stabilization around a L/R ratio of approximately 0.2. For a |$\sigma =20\%$| or above (SMFP ∼1000–500 m), the L/R ratio reaches values of more than 0.6. While this ratio is comparable to the observations of Nishida et al. (2008), such scattering parameters are on the high end of what is observed for the Earth’s crust at frequencies in the secondary microseism range (see Table 3).

5.2 L/R ratio as a function of the correlation length a

Next, we consider the dependence of the L/R ratio on the correlation length a. We investigate three different correlation lengths of 50, 200 and 500 m and a constant σ of 15% (see Table 2). We use a MLGR layering type for all models. In Fig. 10, we show only the models with a correlation length of 50 and 500 m. The L/R ratio with a correlation length of 200 m is shown in the second plot of Fig. 9.

Figure 10.

Results of two simulations with varying correlation length a. The top row of plots corresponds to a = 50 m and the bottom row to a = 500 m. The fluctuation strength σ is 15% for both models. Both models have a MLGR layering type. Left-hand panels: L/R ratios as a function of distance normalized by the SMFP (see Table 2). The black line shows the mean over the 14 angular directions (see Fig. 2) and the yellow area shows one standard deviation; centre panels: dispersion curve for the rotations around the transverse axis Φ_T (perpendicular to the source–receiver direction). The black dotted line shows the theoretical Rayleigh wave fundamental mode; right-hand panels: dispersion curve for the rotation around the vertical axis Φ_V. The dotted line indicates the fundamental Love wave mode. The colour scale corresponds to rotation rate. Each Φ_V dispersion curve is normalized by the maximum of the corresponding Φ_T dispersion curve.

As explained in Section 2.1, we consider two different scattering regimes in these calculations (Pyrak-Nolte 2002). For correlation lengths of 50 m (a < λ) and 200 m (a ∼ λ) we are in the large angle scattering regime, while for a correlation length of 500 m (a > λ), we are in the forward scattering regime. As seen in Table 2, the media with a = 200 m or a ∼ λ have the smallest SMFP and therefore maximal interaction of the wavefield and the scatterers. The calculation of the L/R ratios is again performed as described by eqs (6) and (7). In Fig. 10, the L/R ratios as a function of the normalized distance r_S are shown in the left-hand column, while the dispersion curves are presented in the centre and right.

In Fig. 10, we observe that the L/R ratio increases again with longer propagation distances. The maximum L/R ratio of 0.4 (see Fig. 9), is only observed for the model with λ ∼ a (see Table 2). This indicates that the L/R ratio is related to the number of scattering events, where the number of scattering events is maximal for λ ∼ a.

For a a of 500 m the L/R ratio increases faster and reaches a ratio of 0.3. The L/R ratio for a a of 50 m is only 0.2. The larger L/R ratio correlates also with the amount of scattering events.

The dispersion curves exhibit some interesting features. For the models with |$\sigma =15\%; a=50\,{\rm m}$| and |$\sigma =15\%; a=200\,{\rm m}$|⁠, the fundamental Rayleigh and Love wave modes are visible, similar to Fig. 8. The fundamental Love mode is slightly weaker for the case where a < λ, than for a ∼ λ. For the case where the wavefield propagates mainly in the forward scattering regime (a > λ), neither fundamental Rayleigh nor Love waves are clearly pronounced.

5.3 L/R ratio as a function of layering structure

Each implementation of the random medium in Table 2 is calculated for different underlying velocity layer structures.

The simulations in the previous two sections all focused on the MLGR layering type. Here, we investigate the effect of three different layering types, as shown in Fig. 3(a). Fig. 11 shows the results for a gradient (GR), multiple layers (ML) and a single layer (SL). For comparison, the MLGR layering type is shown in Fig. 9.

Figure 11.

L/R ratio as a function of the distance, normalized by the SMFP (see Table 2) for three models with varying layering types (see titles above plots). All models have σ = 15% and correlation length a = 200 m. The black line shows the mean over the 14 angular directions (see Fig. 2) and the yellow area indicates one standard deviation.

The correlation length a and the fluctuation strength σ are kept constant at values of 200 m (a ∼ λ) and |$\sigma =15\%$|⁠. The calculation and the visualization of the L/R ratios are the same as in the previous sections. The dispersion curves for the SL the ML and the GR layering type are shown in Appendix A.

For each layering type we observe an increased L/R ratio with increasing source–receiver distance and a sharp fundamental mode in the dispersion curve, similar to the case where |$\sigma =15\%; a=200\,{\rm m}$| in Fig. 9. The main difference to the MLGR model is that a clear stabilization of the L/R ratio does not occur, even for similar SMFP. The L/R ratio has a similar initial increase rate for the SL- and the ML types, with a maximum L/R ratio between 0.3 and 0.4 after a propagation distance of 5 SMFP. The L/R ratio for the GR model increases slower and reaches a maximum value of around 0.3.

The SMFP is similar for all models, which indicates that the Rayleigh wavefield scattering is not strongly affected by the layering structure. The Love wave generation, on the other hand, does depend on the layering type. For layering types with sharp velocity contrasts, the wavefield is converted to Love waves more effectively, resulting in an increase and stabilization of the L/R ratio. Therefore, we can state that the MLGR model gives the best conditions for the generation of Love waves. This statement is also supported by the dispersion curves: the ratio of the maximum energy in the Love and Rayleigh dispersion curves is largest for the MLGR case.

This part of the study indicates that the layering structure has some effect on both surface wave scattering and on the emergence of stable Love waves. A thorough investigation using a real, 3-D velocity model, including strong lateral variations in layering, should be pursued in future studies.

5.4 Realistic crustal model

In a final model, we simulate a case that more realistically represents the scattering properties of the earth crust at frequencies in the range of the secondary microseism, that is 0.1–0.2 Hz. The strength of crustal scattering is very dependent on the study area and type the of structure, with, for example oceanic crust less strongly scattering than continental crust .

In Table 3, we summarize the scattering mean free path for different studies of crustal scattering properties in different study areas. The SMFP ranges between 125 and 1000 km (e.g. Bianco et al. 2002; Lacombe et al. 2003) which mostly corresponds to studies in the frequency band between 1 and 10 Hz. Some observations of the mean free path exist for lower frequencies. For example, Lee et al. (2003) compared observed with calculated mean squared envelopes, and found a mean free path of approximately 1605 km at a period of 4 s. An even larger SMFP was observed by Sato & Nohechi (2001) in the lower crust, with a mean free path of 10⁶ km determined at 100 s period.

Table 3.

Table summarizing different studies that characterized scattering mean free path (SMFP) in the Earth crust.

Reference	Location	SMFP (km)	Frequency (Hz)
Fehler et al. (1992)	Kanto-Tokai, Japan	∼250	1–2
Leary & Abercrombie (1994)	Southern California	∼500	0.1–1
Sato & Nohechi (2001)	around the world	∼100000	0.01
Bianco et al. (2002)	Apennines, Italy	∼1000	10
Lacombe et al. (2003)	central France	∼125	4–9
Lee et al. (2003)	central asia	∼1605	0.25
Sens-Schönfelder & Wegler (2006)	Germany	∼690	0.2–24
Hillers et al. (2013)	San Jacinto fault, USA	∼ 70–100	0.2

Reference	Location	SMFP (km)	Frequency (Hz)
Fehler et al. (1992)	Kanto-Tokai, Japan	∼250	1–2
Leary & Abercrombie (1994)	Southern California	∼500	0.1–1
Sato & Nohechi (2001)	around the world	∼100000	0.01
Bianco et al. (2002)	Apennines, Italy	∼1000	10
Lacombe et al. (2003)	central France	∼125	4–9
Lee et al. (2003)	central asia	∼1605	0.25
Sens-Schönfelder & Wegler (2006)	Germany	∼690	0.2–24
Hillers et al. (2013)	San Jacinto fault, USA	∼ 70–100	0.2

Table 3.

Table summarizing different studies that characterized scattering mean free path (SMFP) in the Earth crust.

Reference	Location	SMFP (km)	Frequency (Hz)
Fehler et al. (1992)	Kanto-Tokai, Japan	∼250	1–2
Leary & Abercrombie (1994)	Southern California	∼500	0.1–1
Sato & Nohechi (2001)	around the world	∼100000	0.01
Bianco et al. (2002)	Apennines, Italy	∼1000	10
Lacombe et al. (2003)	central France	∼125	4–9
Lee et al. (2003)	central asia	∼1605	0.25
Sens-Schönfelder & Wegler (2006)	Germany	∼690	0.2–24
Hillers et al. (2013)	San Jacinto fault, USA	∼ 70–100	0.2

Reference	Location	SMFP (km)	Frequency (Hz)
Fehler et al. (1992)	Kanto-Tokai, Japan	∼250	1–2
Leary & Abercrombie (1994)	Southern California	∼500	0.1–1
Sato & Nohechi (2001)	around the world	∼100000	0.01
Bianco et al. (2002)	Apennines, Italy	∼1000	10
Lacombe et al. (2003)	central France	∼125	4–9
Lee et al. (2003)	central asia	∼1605	0.25
Sens-Schönfelder & Wegler (2006)	Germany	∼690	0.2–24
Hillers et al. (2013)	San Jacinto fault, USA	∼ 70–100	0.2

In most studies that quantify the fluctuation strength, it is smaller than 3% (e.g. Sens-Schönfelder et al. 2009) but there are also some studies with a σ around 10% (e.g. Holliger & Levander 1992). These high values are rather exceptional, and found in regions with stronger scattering properties such as volcanoes and volcanic arcs.

Finding a suitable value for the correlation length is not straightforward, since not many studies quantify it. Different studies, mainly in the continental crust, find correlation lengths between 0.5 and 2 km. (e.g. Holliger & Levander 1992; Roth 1997; Sens-Schönfelder et al. 2009; Sato et al. 2012). Other studies show a directional dependence of the correlation length, which becomes more important for the oceanic crust.

Note that when numerically simulating a scattering medium, the effective scattering strength experienced by the seismic wavefield does not correspond exactly to the scattering parameters of the input medium. The effective scattering strength is strongly dependent on the grid spacing and on the numerical domain (see Frankel & Clayton 1986). In order to approach realistic crustal scattering values as best as possible, we have therefore tested several input models, and determine the resulting scattering strength by measuring the SMFP of the wavefield. We choose a σ that leads to a SMFP that can be compared with real values. The first σ is 5% which leads to a SMFP of 1222 km. This value is comparable to the SMFP of Lee et al. (2003), who use a similar frequency range to ours. We additionally show results for a σ of 10%, corresponding to a SMFP of ∼600 km, which should be similar to a tectonically active zone (e.g. transition zone between oceanic and continental crust). We choose a constant correlation length a of 2 km. This value is chosen to prevent that it is smaller than the grid spacing (Frenje & Juhlin 2000). We set a central period of 5 s, similar to the period of the secondary microseism. This leads to a central wavelength of 15 km, which puts us in the Rayleigh scattering regime (a < λ). The model size is 1000 km × 1000 km × 250 km, allowing for a propagation distance of 30 wavelengths. We choose a realistic layering structure, shown in Fig. 3(b). The receiver geometry is the same as shown in Fig. 2, but with a micro-array distance of 0.5 km. We also use the same methods to quantify the L/R ratios and the same method to calculate the SMFP as described in Section 4.

The top plot of Fig. 12 shows the results for a σ of 5%, the bottom plot for a σ of 10%. Compared to the results of the parameter study, we see that L/R ratios are much smaller. This is mainly caused by the small correlation length compared to the central wavelength: the seismic waves are only weakly influenced by the comparably small scatterers. Nevertheless, we see a steady increase in L/R for both models. For a σ of 5%, the mean L/R ratio reaches a value of around 0.1 after propagating the length of the model (approximately 30 λ). Due to the large SMFP, this corresponds to only 0.2 SMFP. The maximum L/R ratio for |$\sigma = 5\%$| would explain a quarter of the observed L/R ratio (in e.g. Nishida et al. 2008; Juretzek & Hadziioannou 2016). For the model with stronger scattering, with |$\sigma =10\%$| (Fig. 12, bottom), the seismic waves cover 0.8 SMFP. The L/R ratio seems to start stabilizing after ∼0.6 SMFP. The L/R ratio reaches values of 0.15 and a maximum of 0.2. This could even explain up to a third of the observed Love-to-Rayleigh ratios (Nishida et al. 2008). The dispersion curves show that most energy propagates in the fundamental modes, similar to the case in the parameter study.

Figure 12.

Results of simulations for the real model setup (see Table 2 and Fig. 3) with varying σ. The top row of plots correspond to σ = 5% and the bottom row to σ =10%. All plots are given for a correlation length a= 2 km. Left-hand panels: L/R ratios as a function of distance normalized by the SMFP (see Table 2). The black line shows the mean over the 14 angular directions (see Fig. 2) and the yellow area shows one standard deviation; centre panels: dispersion curve for the rotations around the transverse axis Φ_T (perpendicular to the source–receiver direction). The black dotted line shows the theoretical Rayleigh wave fundamental mode; right-hand panels: dispersion curve for the rotation around the vertical axis Φ_V. The dotted line indicates the fundamental Love wave mode. The colour scale corresponds to rotation rate. Each Φ_V dispersion curve is normalized by the maximum of the corresponding Φ_T dispersion curve.

6 DISCUSSION

In this paper, study the influence of internal scattering along the source–receiver path on the L/R ratio in the secondary microseism, in an effort to understand the origin of Love waves in the observed ocean-generated noise. This is a different approach than the hypothesis that the Love waves are generated directly at the source (e.g. Fukao et al. 2010; Saito 2010). Based on this study, it is likely that a significant contribution of scattering conversions through complex crustal structures is responsible for the differences in directional incidence of Rayleigh and Love waves that have been observed in several studies (e.g. Friedrich et al. 1998; Behr et al. 2013; Juretzek & Hadziioannou 2016).

Due to the nature of the ambient seismic wavefield, we are mainly interested in scattering conversions between surface wave types. In particular, we consider layered media and scattering conversions from Rayleigh towards Love waves. Much of the theory of seismic wave scattering is based on body wave scattering (Sato et al. 2012). Most of the applications of scattering in seismology use approximations to neglect the surface wave scattering part (e.g. Sens-Schönfelder & Wegler 2006), which would make the calculations far more complicated.

Here, we used numerical simulations of scattering in strongly heterogeneous 3-D media to investigate the emergence of surface waves through conversions. Due to the large computational cost, the number of studies relying on numerical simulation of high-frequency, deterministic seismological wave propagation in 3D random media has only recently started increasing (e.g. Hartzell et al. 2010; Imperatori & Mai 2012; Obermann et al. 2016).

Imperatori & Mai (2012) used a finite difference solver to simulate scattering due to wave propagation through a 3-D von Karman type random medium. They aim to explore the effect of scattering on ground motion in a region near a finite source. They find a good fit between numerically simulated and observed waveforms for correlation lengths between several hundred metres and a few kilometres, which is comparable to the correlation length used in our analysis of realistic crustal models, as is their standard deviation range of |$\sigma =5{-}10\%$|⁠. Obermann et al. (2016) simulated the 3-D wavefield random media in order to study the sensitivity of coda waves to localized velocity perturbations. Although they compute the SMFP with a method similar to ours, they used body waves to calculate the maximum of the coherent phases. Comparing SMFP between our two studies for similar scattering properties of the medium (e.g. |$\sigma =10\%$|⁠, correlation length a = 300 m, comparable to our model with a |$\sigma =10\%$| and a correlation length of a = 200 m), we find SMFP differing by 10–40%. This could be attributed to a difference in scattering behaviour between surface waves and body waves, which have different theoretical scattering coefficients (Snieder 1986). It could also be an effect of the different solvers used in both studies. Obermann et al. (2016) use SPECFEM, which is based on the spectral-element method, while we use a finite difference solver.

The first part of this paper is a general parameter study investigating scattering conversions towards surface waves. We chose higher fluctuation strengths than expected in reality, as well as a correlation length close to the central wavelength of the wavefield (a ∼ λ) to ensure strong scattering. We limit the fluctuation strength to a maximum value of 30%, since beyond that, the wavefield is too strongly scattered to allow the retrieval of clear Love wave dispersion curves. The parameter study shows that the fluctuation strength has a stronger influence on the development of the Love to Rayleigh (L/R) ratio than the correlation length or the type of layering. This can be seen in Figs 8, 10 and A1, where the Love wave dispersion curve energy varies stronger in dependence of σ than for the other parameters. We also find a dependence between the SMFP and the L/R ratios. The shorter the SMFP, the more scattering events for a given distance, and thus the larger the final L/R ratio. Generally, we observe that the L/R gradually increases as seismic waves propagate over several SMFP.

For a σ larger than 15% and a a ∼ λ we observe a stabilization of the L/R ratio at around 0.4 (see e.g. centre and bottom plots of Figs 9 and 10, bottom).

The observed stabilization of Love to Rayleigh ratio is reminiscent of the stabilization of the ratio in energy in the horizontal-to-vertical components which characterizes equipartitioning (e.g. Hennino et al. 2001; Sánchez-Sesma & et al. 2011). Equipartitioning denotes the equal energy distribution of the seismic wavefield to all possible modes of P and Swaves (Weaver 1982, 1985). The energy partitioning is no longer dependent on the original source radiation pattern. For an infinite space with a Poisson solid, this ratio is 10.4 (Weaver 1982). In Hennino et al. (2001), experimental and theoretical values of the S-wave energy over P-wave energy (S/P ratio) were found at S/P ≈ 7.2. Margerin et al. (2009) showed that this state can also be reached in a layered half-space. The equipartition ratio in this case decreases to a value of 7, while part of the body wave energy is stored in surface waves. Several studies (Shapiro et al. 2000; Hennino et al. 2001; Margerin et al. 2009) observe that stabilization of the body wave equipartition ratio occurs after a few (4–5) scattering mean free paths propagation distance. Although we are focusing on surface waves, our L/R ratio stabilizes after a similar propagation distance. This suggests that our observations are related to the equipartition regime. However, to explain the low Love-to-Rayleigh ratio would require a theoretical model for coupled body and surface wave multiple scattering, which is not available to date.

For very strong scattering (⁠|$\sigma =20{-}30\%$|⁠), we obtain a maximum L/R ratio of 0.6 (see e.g. Fig. 9). This ratio corresponds to the ratio observed by Nishida et al. (2008) for the secondary microseism around Japan. In a study by Gualtieri et al. (2013), the authors found that the horizontal power spectral energy predicted by their numerically calculated Rayleigh wave energy did not match the observed level. Introducing Love waves, with an energy ratio of approximately 0.6, would be one explanation for the gap.

6.1 Case study for secondary microseism and realistic Earth crustal scattering

As summarized in Table 3, the scattering properties of the crust are typically not very strong on length scales similar to ocean microseism wavelengths (0.1–0.2 Hz: approximately 18–35 km). In our parameter study, we have opted to consider extremely strong scattering properties, in order to make sure we include clear multiple scattering effects. With this strong scattering present, we find Love to Rayleigh ratios of the order of 0.5 (for a maximum fluctuation strength σ of 30%). As can be seen in Table 1, this is on the lower end of the L/R ratios typically observed in the secondary microseismic band. We conclude that even in cases with very strong scattering, we hardly obtain L/R ratios of the same order as those observed in reality. When we now consider a simulation with more realistic, crustal scattering parameters, we obtain a L/R ratio of at most 0.2, after a propagation distance of 0.8 SMFP.

We calculate two models with fluctuation strengths σ of 5 and 10%. For a |$\sigma= 5\%$|⁠, we measure a SMFP of 1222 km, which is close to the MFP of Lee et al. (2003). They found a SMFP of 1605 km for central period of 4 s. We will compare these L/R values to real observations. Juretzek & Hadziioannou (2016) measured the L/R ratio at the Gräfenberg array for noise sources in North Atlantic, and found a L/R ratio around 0.8. They point out the L/R ratio for stations far inland, or behind complex crustal structures (e.g. Alps) are likely more strongly affected by wave type conversions along the propagation path. The distance from the French Atlantic coast to the Gräfenberg array is around 1200 km. Therefore, assuming a weakly scattering continental crust and negligible scattering in the oceanic crust, the wavefield travels almost one SMFP through the continental crust from the coast to the receiver. Our realistic crustal models have a size of 500 km, which is a little bit less than half this distance. Our maximum L/R ratio, obtained for |$\sigma = 10\%$|⁠, is 0.15. Taking twice this value to account for the larger propagation distance, we would obtain a ratio of 0.3, still only less than the half of the ratio observed by Juretzek & Hadziioannou (2016).

For noise sources in the Pacific Ocean, recorded, for example in Japan, the wavefield has to pass the ring of fire to reach the continental crust. In such a region with many volcanic structures, the heterogeneity is stronger than in tectonically less active regions (Takahashi et al. 2007; Sato et al. 2012). To mimic wave propagation across the ring of fire, we suppose a |$\sigma =10\%$|⁠, a SMFP of 620 km, and 250 km for the width of the ring of fire. The wavefield must therefore propagate approximately half a SMFP of more strongly scattering crust to reach the continent. In Fig. 12, we find that the L/R ratio reaches a value of approximately 0.2 after half a SMFP propagation distance. Assuming that the wavefield again propagates about 1000 km through continental crust before reaching the seismometer, we would expect a L/R ratio of the order of 0.5. This value is quite close to the L/R ratio of Nishida et al. (2008) which has been observed in Japan (see Table 1).

It must be noted that in this study, we consider only fundamental mode Rayleigh and Love waves. Some studies have also observed higher modes in the secondary seismic noise (e.g. Kimman et al. 2012). The presence of higher modes would also influence the L/R ratio, by affecting the energy redistribution upon scattering events. Nevertheless, the study of Gualtieri et al. (2013) has shown that the fundamental modes stores a major portion of the energy.

We conclude that even considering relatively strongly scattering crustal features, wave-type conversions alone cannot account for the observed L/R ratios in the secondary microseism. Other factors to increase the L/R ratio must be taken into account. Possible mechanisms include conversions due to surface topography, complex layering structures or sedimentary basins. Wave type conversions as the wavefield reaches the continental slope could also have an effect on the L/R ratio observed on land.

7 CONCLUSION

The aim of this study is to investigate one of the possible mechanisms for generating Love waves in the secondary microseismic noise field. We use a vertical pressure source on a flat seafloor, to mimic the source mechanism proposed by Longuet-Higgins (1950), which mainly generates Rayleigh waves and other vertically polarized waves locally around the source. Through scattering conversions along the propagation path, these vertically polarized waves are converted towards Love waves. We aim to quantify the heterogeneity of the medium necessary to generate a similar ratio of Love to Rayleigh waves as that observed in the secondary microseism.

To that end, we model a vertical point source to initially generate a wavefield composed of vertically polarized waves, which then propagates in a random medium. We measure the wavefield with lines of micro-arrays which extend radially from the source. With these measurements, we evaluate the Love to Rayleigh ratio as a function of distance from the source, for different random media.

This study consists of two parts. The first part considers the influence of different scattering parameters (fluctuation strength σ, correlation length a, layering structure) on the L/R ratio. Here, we model wave propagation in random media with increasingly strong scattering, including cases with scattering stronger than what would be expected for realistic crustal structures. As can be seen from the dispersion curves measured in our various model implementations, the scattered wavefield still contains a significant portion of fundamental Love and Rayleigh modes.

Our parameter study shows that fluctuation strength σ has the largest influence on the L/R ratio, while variations in layering structure and the correlation length a are less effective. We show that for a strong fluctuation strength, the mean L/R ratio stabilizes around a value of 0.4 after approximately 4 SMFP. The maximum L/R obtained is approximately 0.6, for σ = 20%, and for a ∼ λ, where λ is the central frequency of the Rayleigh waves.

In the second part of the study, we consider realistic scattering parameters for the Earth’s crust, and propagation distances comparable to ocean—mid-continent paths. For a frequency range centred around 0.2 Hz, we obtain L/R ratios up to 0.2, accounting for less than half of the Love wave energy typically observed in the secondary microseism in most studies to date. We conclude that wave conversions and scattering off internal heterogeneities may contribute in a non-negligible way to the observed Love-to-Rayleigh ratios in the secondary microseism. However, scattering effects alone are not sufficient, and additional, different generation mechanisms for Love waves are necessary to explain the observed L/R values.

ACKNOWLEDGEMENTS

Computing resources were provided by the Leibniz Supercomputing Centre (LRZ, projects no. h019z and pr63qo on SuperMUC). The authors wish to thank our system administrator Jens Oeser for support. We are grateful to Stefan Wenk, Carina Juretzek, Josselin Garnier, Ludovic Margerin and Roel Snieder for constructive discussions and comments. We thank the editor, Anne Obermann and an anonymous reviewer for their thorough review and comments, which helped improve the paper.

This work was funded by the Emmy Noether program (HA7019/1-1) of the German Research Foundation (DFG), and has benefited from discussions within COST Action ES1401-TIDES, supported by COST (European Cooperation in Science and Technology).

REFERENCES

Aki

K.

,

Christoffersson

A.

,

Husebye

E.S.

,

1976

.

Three-dimensional seismic structure of the lithosphere under Montana LASA

,

Bull. seism. Soc. Am.

,

66

(

2

),

501

–

524

.

Bannister

S.C.

,

Husebye

E.S.

,

Ruud

B.O.

,

1990

.

Teleseismic P coda analyzed by three-component and array techniques: deterministic location of topographic P-to-RG scattering near the Noress array

,

Bull. seism. Soc. Am.

,

80

(

6B

),

1969

.

Behr

Y.

,

Townend

J.

,

Bowen

M.

,

Carter

L.

,

Gorman

R.

,

Brooks

L.

,

Bannister

S.

,

2013

.

Source directionality of ambient seismic noise inferred from three-component beamforming

,

J. geophys. Res.: Solid Earth

,

118

(

1

),

240

–

248

.

10.1029/2012JB009382

10.1046/j.1365-246X.2002.01696.x

Bianco

F.

,

Del Pezzo

E.

,

Castellano

M.

,

Ibanez

J.

,

Di Luccio

F.

,

2002

.

Separation of intrinsic and scattering seismic attenuation in the southern Apennine zone, Italy

,

J. geophys. Int.

,

150

(

1

),

10

–

22

.

10.1016/S0098-3004(02)00006-7

Bohlen

T.

,

2002

.

Parallel 3-D viscoelastic finite difference seismic modelling

,

Comput. Geosci.

,

8

,

887

–

899

.

10.1046/j.1365-246X.1998.00410.x

Dolan

S.

,

Bean

C.

,

Riollet

B.

,

1998

.

The broad-band fractal nature of heterogeneity in the upper crust from petrophysical logs

,

J. geophys. Int.

,

132

(

3

),

489

–

507

.

10.1111/j.1365-246X.1992.tb03470.x

Fehler

M.

,

Hoshiba

M.

,

Sato

H.

,

Obara

K.

,

1992

.

Separation of scattering and intrinsic attenuation for the Kanto-Tokai region, Japan, using measurements of S-wave energy versus hypocentral distance

,

J. geophys. Int.

,

108

(

3

),

787

–

800

.

Frankel

A.

,

Clayton

R.W.

,

1986

.

Finite difference simulations of seismic scattering: implications for the propagation of short-period seismic waves in the crust and models of crustal heterogeneity

,

J. geophys. Res.: Solid Earth

,

91

(

B6

),

6465

–

6489

.

10.1029/JB091iB06p06465

10.1016/S0926-9851(00)00003-3

Frenje

L.

,

Juhlin

C.

,

2000

.

Scattering attenuation: 2-D and 3-D finite difference simulations versus theory

,

J. Appl. Aeophys.

,

44

(

1

),

33

–

46

.

Friedrich

A.

,

Krüger

F.

,

Klinge

K.

,

1998

.

Ocean-generated microseismic noise located with the Gräfenberg array

,

J. Seismol.

,

2

(

1

),

47

–

64

.

Froment

B.

,

Campillo

M.

,

Roux

P.

,

Gouédard

P.

,

Verdel

A.

,

Weaver

R.L.

,

2010

.

Estimation of the effect of nonisotropically distributed energy on the apparent arrival time in correlations

,

Geophysics

,

75

(

5

),

SA85

–

SA93

.

Fukao

Y.

,

Nishida

K.

,

Kobayashi

N.

,

2010

.

Seafloor topography, ocean infragravity waves, and background Love and Rayleigh waves

,

J. geophys. Res.: Solid Earth

,

115

(

B4

),

B04302

.

Gualtieri

L.

,

Stutzmann

E.

,

Capdeville

Y.

,

Ardhuin

F.

,

Schimmel

M.

,

Mangeney

A.

,

Morelli

A.

,

2013

.

Modelling secondary microseismic noise by normal mode summation

,

J. geophys. Int.

,

193

(

3

),

1732

–

1745

.

10.1093/gji/ggt090

Gualtieri

L.

,

Stutzmann

É.

,

Farra

V.

,

Capdeville

Y.

,

Schimmel

M.

,

Ardhuin

F.

,

Morelli

A.

,

2014

.

Modelling the ocean site effect on seismic noise body waves

,

J. geophys. Int.

,

197

(

2

),

1096

–

1106

.

10.1093/gji/ggu042

Gutenberg

B.

,

1936

.

On microseisms

,

Bull. seism. Soc. Am.

,

26(2)

,

111

–

117

.

Hartzell

S.

,

Harmsen

S.

,

Frankel

A.

,

2010

.

Effects of 3D random correlated velocity perturbations on predicted ground motions

,

Bull. seism. Soc. Am.

,

100

(

4

),

1415

–

1426

.

10.1785/0120090060

Hasselmann

K.

,

1963

.

A statistical analysis of the generation of microseisms

,

Rev. Geophys.

,

1

(

2

),

177

–

210

.

10.1029/RG001i002p00177

10.1103/PhysRevLett.86.3447

Hennino

R.

,

Trégourès

N.

,

Shapiro

N.M.

,

Margerin

L.

,

Campillo

M.

,

van Tiggelen

B.A.

,

Weaver

R.L.

,

2001

.

Observation of equipartition of seismic waves

,

Phys. Rev. Lett.

,

86

,

3447

–

3450

.

Hillers

G.

,

Graham

N.

,

Campillo

M.

,

Kedar

S.

,

Lands

M.

,

Shapiro

N.

,

2012

.

Global oceanic microseism sources as seen by seismic arrays and predicted by wave action models

,

Geochem., Geophys., Geosyst.

,

13

(

1

),

Q01021

.

Hillers

G.

,

Ben-Zion

Y.

,

Landès

M.

,

Campillo

M.

,

2013

.

Interaction of microseisms with crustal heterogeneity: a case study from the San Jacinto fault zone area

,

Geochem. Geophys. Geosyst.

,

14

(

7

),

2182

–

2197

.

10.1002/ggge.20140

10.1111/j.1365-2478.1987.tb00841.x

Holberg

O.

,

1987

.

Computational aspects of the choice of operator and sampling interval for numerical differentiation in large-scale simulation of wave phenomena

,

Geophys. Prospect.

,

35

(

6

),

629

–

655

.

10.1111/j.1365-246X.1996.tb06025.x

Holliger

K.

,

1996

.

Upper-crustal seismic velocity heterogeneity as derived from a variety of P-wave sonic logs

,

J. geophys. Int.

,

125

(

3

),

813

–

829

.

10.1111/j.1365-246X.1997.tb04071.x

Holliger

K.

,

1997

.

Seismic scattering in the upper crystalline crust based on evidence from sonic logs

,

J. geophys. Int.

,

128

(

1

),

65

–

72

.

Holliger

K.

,

Levander

A.R.

,

1992

.

A stochastic view of lower crustal fabric based on evidence from the Ivrea Zone

,

Geophys. Res. Lett.

,

19

(

11

),

1153

–

1156

.

10.1029/92GL00919

Huang

J.-W.

,

Bohlen

T.

,

Milkereit

B.

,

2006

.

Numerical solutions of seismic scattering in heterogeneous media

,

Can. SEG Record.

,

31

(

9

),

44

–

47

.

10.1111/j.1365-246X.2011.05296.x

Imperatori

W.

,

Mai

P.M.

,

2012

.

Sensitivity of broad-band ground-motion simulations to earthquake source and Earth structure variations: an application to the Messina straits (Italy)

,

J. geophys. Int.

,

188

(

3

),

1103

–

1116

.

Juretzek

C.

,

Hadziioannou

C.

,

2016

.

Where do ocean microseisms come from? A study of Love-to-Rayleigh wave ratios

,

J. geophys. Res.: Solid Earth

,

121

(

9

),

6741

–

6756

.

10.1002/2016JB013017

Kedar

S.

,

Longuet-Higgins

M.

,

Webb

F.

,

Graham

N.

,

Clayton

R.

,

Jones

C.

,

2008

.

The origin of deep ocean microseisms in the North Atlantic ocean

,

Proc. R. Soc. Lond., A

,

464

,

777

–

793

., .

10.1098/rspa.2007.0277

Google Preview

10.1046/j.1365-246X.1998.1331488.x

Kennett

B.

,

1998

.

Guided waves in three-dimensional structures

,

J. geophys. Int.

,

133

(

1

),

159

–

174

.

10.1111/j.1365-246X.1984.tb02853.x

Kennett

B.L.N.

,

1984

.

Guided wave propagation in laterally varying media I. Theoretical development

,

Geophys. J. R. astr. Soc.

,

79

(

1

),

235

–

255

.

Kimman

W.

,

Campman

X.

,

Trampert

J.

,

2012

.

Characteristics of seismic noise: fundamental and higher mode energy observed in the northeast of the Netherlands

,

Bull. seism. Soc. Am.

,

102

(

4

),

1388

–

1399

.

10.1785/0120110069

Kumagai

H.

,

Saito

T.

,

O’Brien

G.

,

Yamashina

T.

,

2011

.

Characterization of scattered seismic wavefields simulated in heterogeneous media with topography

,

J. geophys. Res.: Solid Earth

,

116

(

B3

),

doi:10.1029/2010JB007718

.

10.1046/j.1365-246X.2003.01976.x

Lacombe

C.

,

Campillo

M.

,

Paul

A.

,

Margerin

L.

,

2003

.

Separation of intrinsic absorption and scattering attenuation from Lg coda decay in central France using acoustic radiative transfer theory

,

J. geophys. Int.

,

154

(

2

),

417

–

425

.

Leary

P.

,

Abercrombie

R.

,

1994

.

Frequency dependent crustal scattering and absorption at 5160 Hz from coda decay observed at 2.5 km Depth

,

Geophys. Res. Lett.

,

21

(

11

),

971

–

974

.

10.1029/94GL00977

Lee

W.S.

,

Sato

H.

,

Lee

K.

,

2003

.

Estimation of S-wave scattering coefficient in the mantle from envelope characteristics before and after the ScS arrival

,

Geophys. Res. Lett.

,

30

(

24

),

2248

.

Levander

A.R.

,

Holliger

K.

,

1992

.

Small-scale heterogeneity and large-scale velocity structure of the continental crust

,

J. geophys. Res.: Solid Earth

,

97

(

B6

),

8797

–

8804

.

10.1029/92JB00659

10.1111/j.1365-246X.2007.03414.x

Lin

F.-C.

,

Ritzwoller

M.H.

,

Townend

J.

,

Bannister

S.

,

Savage

M.K.

,

2007

.

Ambient noise Rayleigh wave tomography of New Zealand

,

J. geophys. Int.

,

170

(

2

),

649

–

666

.

Longuet-Higgins

M.

,

1950

.

A theory of the origin of microseisms

,

Phil. Trans. R. Soc. Lond. A

,

243

(

857

),

1

–

35

.

10.1111/j.1365-246X.2007.03603.x

Maeda

T.

,

Sato

H.

,

Nishimura

T.

,

2008

.

Synthesis of coda wave envelopes in randomly inhomogeneous elastic media in a halfspace: single scattering model including Rayleigh waves

,

J. geophys. Int.

,

172

(

1

),

130

–

154

.

10.1111/j.1365-246X.2008.04068.x

Margerin

L.

,

Campillo

M.

,

Van Tiggelen

B.A.

,

Hennino

R.

,

2009

.

Energy partition of seismic coda waves in layered media: theory and application to Pinyon Flats Observatory

,

J. geophys. Int.

,

177

,

571

–

585

.

10.1046/j.1365-246x.2001.01460.x

Maupin

V.

,

2001

.

A multiple-scattering scheme for modelling surface wave propagation in isotropic and anisotropic three-dimensional structures

,

J. geophys. Int.

,

146

(

2

),

332

–

348

.

Nishida

K.

,

Kawakatsu

H.

,

Fukao

Y.

,

Obara

K.

,

2008

.

Background Love and Rayleigh waves simultaneously generated at the Pacific Ocean floors

,

Geophys. Res. Lett.

,

35

(

16

),

L16307

.

10.1029/2008GL034753

Obermann

A.

,

Planès

T.

,

Hadziioannou

C.

,

Campillo

M.

,

2016

.

Lapse-time-dependent coda-wave depth sensitivity to local velocity perturbations in 3-D heterogeneous elastic media

,

J. geophys. Int.

,

207

(

1

),

59

–

66

.

10.1093/gji/ggw264

10.1016/j.jvolgeores.2009.03.004

O’Brien

G.S.

,

Bean

C.J.

,

2009

.

Volcano topography, structure and intrinsic attenuation: their relative influences on a simulated 3D visco-elastic wavefield

,

J. Volc. Geotherm. Res.

,

183

(

1-2

),

122

–

136

.

Pyrak-Nolte

L.

,

2002

.

Seismic imaging of fractured media

, in

Proceedings of the 5th International Workshop on the Application of Geophysics in Rock Engineering

,

Toronto, Canada

,

5-13

.

Google Preview

Rind

D.

,

Down

W.L.

,

1979

.

Microseisms at Palisades: 2. Rayleigh wave and Love wave characteristics and the geologic control of propagation

,

J. geophys. Res.: Solid Earth

,

84

(

B10

),

5632

–

5642

.

10.1029/JB084iB10p05632

10.1016/S0031-9201(97)00048-4

Roth

M.

,

1997

.

Statistical interpretation of traveltime fluctuations

,

Phys. Earth planet. Inter.

,

104

,

213

–

228

.

10.1111/j.1365-246X.2010.04695.x

Saito

T.

,

2010

.

Love-wave excitation due to the interaction between a propagating Ocean wave and the sea-bottom topography

,

J. geophys. Int.

,

182

(

3

),

1515

–

1523

.

10.1111/j.1365-246X.2011.05064.x

Sánchez-Sesma

F.J.

et al.,

2011

.

A theory for microtremor H/V spectral ratio: application for a layered medium

,

J. geophys. Int.

,

186

(

1

),

221

–

225

.

Sato

H.

,

Nohechi

M.

,

2001

.

Envelope formation of long-period Rayleigh waves in vertical component seismograms: single isotropic scattering model

,

J. geophys. Res.: Solid Earth

,

106

(

B4

),

6589

–

6594

.

10.1029/2000JB900383

10.1111/j.1365-246X.2006.03139.x

Sato

H.

,

Fehler

M.

,

Maeda

T.

,

2012

,

Seismic Wave Propagation and Scattering in the Heterogeneous Earth

,

Springer Verlag

.

Sens-Schönfelder

C.

,

Wegler

U.

,

2006

.

Radiative transfer theory for estimation of the seismic moment

,

J. geophys. Int.

,

167

(

3

),

1363

–

1372

.

Sens-Schönfelder

C.

,

Margerin

L.

,

Campillo

M.

,

2009

.

Laterally heterogeneous scattering explains Lg blockage in the Pyrenees

,

J. geophys. Res.: Solid Earth

,

114

(

B7

),

B07309

.

Shapiro

N.M.

,

Campillo

M.

,

Margerin

L.

,

Singh

S.K.

,

Kostoglodov

V.

,

Pacheco

J.

,

2000

.

The energy partitioning and the diffusive character of the seismic coda

,

Bull. seism. Soc. Am.

,

90

(

3

),

655

.

10.1785/0119990021

Shapiro

N.M.

,

Campillo

M.

,

Stehly

L.

,

Ritzwoller

M.H.

,

2005

.

High-resolution surface-wave tomography from ambient seismic noise

,

Science

,

307

(

5715

),

1615

–

1618

.

10.1111/j.1365-246X.1986.tb04372.x

Snieder

R.

,

1986

.

3-D linearized scattering of surface waves and a formalism for surface wave holography

,

Geophys. J. R. Astr. Soc

,

84

,

581

–

605

.

Spudich

P.

,

Fletcher

J.B.

,

2008

.

Observation and prediction of dynamic ground strains, tilts, and torsions caused by the Mw 6.0 2004 Parkfield, California, earthquake and aftershocks, aerived from UPSAR array observations

,

Bull. seism. Soc. Am.

,

98

(

4

),

1898

–

1914

.

10.1785/0120070157

Spudich

P.

,

Steck

L.K.

,

Hellweg

M.

,

Fletcher

J.

,

Baker

L.

,

1995

.

Transient stresses at Parkfield, California, produced by the M 7.4 Landers earthquake of June 28, 1992: observations from the UPSAR dense seismograph array

,

J. geophys. Res.

,

100

,

675

–

690

.

10.1029/94JB02477

Stutzmann

E.

,

Roult

G.

,

Astiz

L.

,

2000

.

GEOSCOPE station noise levels

,

Bull. seism. Soc. Am.

,

90

(

3

),

690

–

701

.

10.1785/0119990025

10.1111/j.1365-246X.2012.05638.x

Stutzmann

E.

,

Ardhuin

F.

,

Schimmel

M.

,

Mangeney

A.

,

Patau

G.

,

2012

.

Modelling long-term seismic noise in various environments

,

J. geophys. Int.

,

191

(

2

),

707

–

722

.

10.1111/j.1365-246X.2006.03197.x

Takahashi

T.

,

Sato

H.

,

Nishimura

T.

,

Obara

K.

,

2007

.

Strong inhomogeneity beneath quaternary volcanoes revealed from the peak delay analysis of S-wave seismograms of microearthquakes in northeastern Japan

,

J. geophys. Int.

,

168

(

1

),

90

–

99

.

Takemura

S.

,

Furumura

T.

,

Maeda

T.

,

2015

.

Scattering of high-frequency seismic waves caused by irregular surface topography and small-scale velocity inhomogeneity

,

J. geophys. Int.

,

201

(

1

),

459

–

474

.

10.1093/gji/ggv038

Tanimoto

T.

,

Hadziioannou

C.

,

Igel

H.

,

Wasserman

J.

,

Schreiber

U.

,

Gebauer

A.

,

2015

.

Estimate of Rayleigh-to-Love wave ratio in the secondary microseism by colocated ring laser and seismograph

,

Geophys. Res. Lett.

,

42

(

8

),

2650

–

2655

.

10.1002/2015GL063637

Tanimoto

T.

,

Lin

C.-J.

,

Hadziioannou

C.

,

Igel

H.

,

Vernon

F.

,

2016

.

Estimate of Rayleigh-to-Love wave ratio in the secondary microseism by a small array at Piñon Flat observatory, California

,

Geophys. Res. Lett.

,

43

(

21

),

11 173

–

11 181

.

Traer

J.

,

Gerstoft

P.

,

Bromirski

P.D.

,

Shearer

P.M.

,

2012

.

Microseisms and hum from ocean surface gravity waves

,

J. geophys. Res.: Solid Earth

,

117

(

B11

),

B11307

.

10.1111/j.1365-246X.2009.04239.x

Tsai

V.C.

,

2009

.

On establishing the accuracy of noise tomography travel-time measurements in a realistic medium

,

J. geophys. Int.

,

178

(

3

),

1555

–

1564

.

Weaver

R.

,

Froment

B.

,

Campillo

M.

,

2009

.

On the correlation of non-isotropically distributed ballistic scalar diffuse waves

,

J. acoust. Soc. Am.

,

126

(

4

),

1817

–

1826

.

10.1121/1.3203359

Weaver

R.L.

,

1982

.

On diffuse waves in solid media

,

J. acoust. Soc. Am.

,

71

(

6

),

1608

–

1609

.

10.1121/1.387816

Weaver

R.L.

,

1985

.

Diffuse elastic waves at a free surface

,

J. acoust. Soc. Am.

,

78

(

1

),

131

–

136

.

10.1121/1.392576

Weaver

R.L.

,

Yoritomo

J.Y.

,

2018

.

Temporally weighting a time varying noise field to improve Green function retrieval

,

J. acoust. Soc. Am.

,

143

(

6

),

3706

–

3719

.

10.1121/1.5043406

Webb

S.C.

,

1998

.

Broadband seismology and noise under the ocean

,

Rev. Geophys.

,

36

(

1

),

105

–

142

.

10.1029/97RG02287

Wu

R.-S.

,

Xu

Z.

,

Li

X.-P.

,

1994

.

Heterogeneity spectrum and scale-anisotropy in the upper crust revealed by the German Continental Deep-Drilling (KTB) Holes

,

Geophys. Res. Lett.

,

21

(

10

),

911

–

914

.

10.1029/94GL00772

10.1111/j.1365-246X.2006.03028.x

Yao

H.

,

van Der Hilst

R.D.

,

De Hoop

M.V.

,

2006

.

Surface-wave array tomography in SE Tibet from ambient seismic noise and two-station analysis—I. Phase velocity maps

,

J. geophys. Int.

,

166

(

2

),

732

–

744

.

Zhan

Z.

,

Tsai

V.C.

,

Clayton

R.W.

,

2013

.

Spurious velocity changes caused by temporal variations in ambient noise frequency content

,

J. geophys. Int.

,

194

(

3

),

1574

–

1581

.

10.1093/gji/ggt170

APPENDIX A: DISPERSION CURVES

In Fig. A1, we show the dispersion curves for models with two different layering types, as discussed in Section 5.3.

Figure A1.

Dispersion curves for the transverse (left-hand) and the vertical curl (rigth-hand) for the single layering type (top row), the multiple layering type (middle row) and the gradient layering type (bottom row). The correlation length a is 200 m and the scattering strength σ is 15%.