A scalar radiative transfer model including the coupling between surface and body waves

TL;DR

This paper develops a scalar radiative transfer model that captures the coupling between surface and body seismic waves, providing insights into energy transport and coda wave formation in Earth's subsurface.

Contribution

It introduces a coupled transport equation framework for surface and body waves, including depth-dependent scattering and wave coupling effects, which is novel in seismic modeling.

Findings

Equipartition between surface and body waves at long times

Depth-dependent scattering mean free paths for body waves

Coupling influences coda wave characteristics

Abstract

To describe the energy transport in the seismic coda, we introduce a system of radiative transfer equations for coupled surface and body waves in a scalar approximation. Our model is based on the Helmholtz equation in a half-space geometry with mixed boundary conditions. In this model, Green's function can be represented as a sum of body waves and surface waves, which mimics the situation on Earth. In a first step, we study the single-scattering problem for point-like objects in the Born approximation. Using the assumption that the phase of body waves is randomized by surface reflection or by interaction with the scatterers, we show that it becomes possible to define, in the usual manner, the cross-sections for surface-to-body and body-to-surface scattering. Adopting the independent scattering approximation, we then define the scattering mean free paths of body and surface waves including the coupling between the two types of waves. Using a phenomenological approach, we then derive a set of coupled transport equations satisfied by the specific energy density of surface and body waves in a medium containing a homogeneous distribution of point scatterers. In our model, the scattering mean free path of body waves is depth dependent as a consequence of the body-to-surface coupling. We demonstrate that an equipartition between surface and body waves is established at long lapse-time, with a ratio which is predicted by usual mode counting arguments. We derive a diffusion approximation from the set of transport equations and show that the diffusivity is both anisotropic and depth dependent. The physical origin of the two properties is discussed. Finally, we present Monte-Carlo solutions of the transport equations which illustrate the convergence towards equipartition at long lapse-time as well as the importance of the coupling between surface and body waves in the generation of coda waves.

Figures

Equations

$$\left(\rho {\partial }_{tt}-T\Delta \right)u(t,\mathbf{R})=0$$

$$\begin{array}{r}{({\partial }_z+\alpha )u(t,\mathbf{R})|}_{z=0}=0\\ +\text{radiation condition at }\infty \end{array}$$

$$L=\frac{1}{2}\left[\rho ({\partial }_{t}u(t,\mathbf{R}){)}^2-T(\nabla u(t,\mathbf{R}){)}^2+\alpha Tu(t,\mathbf{R}{)}^2\delta (z)\right]$$

$$\Delta u(\mathbf{R})+\frac{{\omega }^2}{c^2}u(\mathbf{R})=0$$

$$\mathbf{J}=$$

$$w=$$

$$u_s(\mathbf{r},z)=\sqrt{2\alpha }{e}^{-\alpha z}\frac{e^{i{\mathbf{k}}_{\parallel }\cdot \mathbf{r}}}{2\pi }$$

$$u_b(\mathbf{r},z)=\frac{1}{(2\pi {)}^{3/2}}(e^{-iqz}+r(q)e^{iqz})e^{i{\mathbf{k}}_{\parallel }\cdot \mathbf{r}},q\ge 0$$

$$r(q)=\frac{q+i\alpha }{q-i\alpha }$$

$${\psi }_b(\hat{\mathbf{n}},z)=(e^{-ikz\hat{\mathbf{n}}\cdot \hat{\mathbf{z}}}+r(k\hat{\mathbf{n}}\cdot \hat{\mathbf{z}})e^{ikz\hat{\mathbf{n}}\cdot \hat{\mathbf{z}}})$$

$$c_{\varphi }=\frac{c}{\sqrt{1+\frac{c^2{\alpha }^2}{{\omega }^2}}}$$

$$v_g=\frac{d\omega }{dk}=\frac{c^2}{c_{\varphi }}=c\sqrt{1+\frac{c^2{\alpha }^2}{{\omega }^2}}$$

$${\mathbf{v}}_s^E=\frac{⟨\mathbf{J}⟩}{⟨w⟩}=v_g{\hat{\mathbf{k}}}_{\parallel }$$

$${\mathbf{v}}_b^E=\underset{h\to \infty }{\lim }\frac{⟨\mathbf{J}{⟩}_h}{⟨w{⟩}_h}=c\sin j{\hat{\mathbf{k}}}_{\parallel },$$

$$\begin{array}{rl}G(\mathbf{r},z,z_0)=& \frac{1}{(2\pi {)}^3}{\int }_0^{+\infty }dq{\int }_{{\mathbb{R}}^2}\frac{e^{i{\mathbf{k}}_{\parallel }\cdot \mathbf{r}}(e^{-iqz}+r(q)e^{iqz})(e^{-iq{z}_0}+r(q)e^{iq{z}_0}{)}^{\ast }}{k^2-k_{\parallel }^2-q^2+iϵ}{d}^2{k}_{\parallel }\\ & +\frac{2\alpha e^{-\alpha (z+z_0)}}{(2\pi {)}^2}{\int }_{{\mathbb{R}}^2}\frac{e^{i{\mathbf{k}}_{\parallel }\cdot \mathbf{r}}}{k^2+{\alpha }^2-k_{\parallel }^2+iϵ}{d}^2{k}_{\parallel }\end{array},$$

$$G(\mathbf{r},z,z_0)=-\frac{e^{ikR}}{4\pi R}{\psi }_b(\hat{\mathbf{R}},z_0)-\frac{\alpha e^{-\alpha (z+z_0)+i{k}_{s}r+i\pi /4}}{\sqrt{2\pi k_{s}r}}$$

$$\begin{array}{rl}{\mathcal{E}}_b(z_0)=& \frac{\rho {\omega }^{2}c}{2}{\int }_H|G_b(\mathbf{r},z,z_0){|}^2{R}^{2}d\hat{R}\\ =& \frac{\rho {\omega }^{2}c}{32{\pi }^2}{\int }_{2\pi }|{\psi }_b(\hat{\mathbf{R}},z_0){|}^{2}d\hat{R},\end{array}$$

$$⟨{\mathcal{E}}_b⟩=\frac{\rho {\omega }^{2}c}{8\pi }$$

$$\begin{array}{rl}{\mathcal{E}}_s(z_0)=& \frac{\rho {\omega }^2{v}_g}{2}{\int }_C|G_s(\mathbf{r},z,z_0){|}^{2}rd\varphi dz\\ =& \frac{\rho {\omega }^2{v}_g{\alpha }^2{e}^{-2\alpha z_0}}{4\pi k_s}{\int }_0^{2\pi }d\varphi {\int }_0^h{e}^{-2\alpha z}dz\\ =& \frac{\rho \omega c^2\alpha }{4}{e}^{-2\alpha z_0}\end{array}$$

$$R(z_0)=\frac{{\mathcal{E}}_s(z_0)}{⟨{\mathcal{E}}_b⟩}=\frac{2\pi c\alpha }{\omega }{e}^{-2\alpha z_0}$$

$$n_{s,b}(z_0)=-{\frac{\text{Im}{G}_{s,b}(\mathbf{r},z_0,z_0)}{\pi }|}_{\mathbf{r}=0,z=z_0}\times \frac{d{k}_{s,b}^2(\omega )}{d\omega },$$

$$n_b(z_0)=$$

$$⟨n_b⟩=$$

$$n_s(z_0)=$$

$$\mathcal{R}={\int }_0^{\infty }\frac{n_s(z)}{⟨n_b⟩}dz=\frac{\pi c}{\omega }$$

$$\Delta u(\mathbf{r},z)+k^2(1+ϵa^3\delta (\mathbf{r})\delta (z-z_s))u(\mathbf{r},z)=0$$

$$u(\mathbf{r},z)=u_0(\mathbf{r},z)-k^2ϵa^{3}G(\mathbf{r},z,z_s)u_0(0,z_s)$$

$$U^{s\to b}=\frac{\rho {\omega }^2{S}^{2}c}{2(4\pi {)}^2{R}^2}{\int }_{2\pi }|{\psi }_b(\hat{\mathbf{R}},z_s){|}^2{R}^{2}d\hat{R}.$$

$$|{\mathbf{J}}_s|=\frac{\rho {\omega }^2{v}_g}{4\alpha }.$$

$${\sigma }^{s\to b}(z_s)=\frac{\alpha c{k}^4{a}^6{ϵ}^2{e}^{-2\alpha z_s}}{8{\pi }^2{v}_g}{\int }_{2\pi }|{\psi }_b(\hat{\mathbf{R}},z_s){|}^{2}d\hat{R}$$

$$⟨{\sigma }^{s\to b}(z_s)⟩=\frac{c\alpha k^4{a}^6{ϵ}^2{e}^{-2\alpha z_s}}{2\pi v_g}$$

$$\begin{array}{rl}{U}^{s\to s}=& \frac{\rho \omega c^2{\alpha }^2{S}^2}{4\pi }{\int }_{2\pi }{\int }_0^{+\infty }{e}^{-2\alpha (z+z_s)}d\varphi dz\\ =& \frac{\rho \omega c^2\alpha S^2{e}^{-2\alpha z_s}}{4}\end{array}$$

Full text

1 Introduction

In seismology, Radiative Transfer (RT) has been used for more than three decades to characterize the scattering and absorption properties of Earth’s crust (see for instance Fehler et al., 1992; Hoshiba, 1993; Carcolé & Sato, 2010; Eulenfeld & Wegler, 2017). Since its introduction by Wu (1985) for scalar waves in the stationary regime, RT has been considerably improved to bring it in closer agreement with real-world applications. In particular, the coupling between shear and compressional waves was developed in a series of papers by Weaver (1990); Turner & Weaver (1994); Zeng (1993); Sato (1994); Ryzhik et al. (1996). The model of Sato (1994) was applied to data from an active experiment by Yamamoto & Sato (2010) and showed impressive agreement between observations and elastic RT theory. For comprehensive introductions to RT, the reader is referred to the review chapter by Margerin (2005) or the monograph of Sato et al. (2012).

Parallel to the physical and mathematical developments of the theory, more and more realistic Monte-Carlo simulations of the transport process were developed over the years. This includes, for example, the treatment of non-isotropic scattering (Abubakirov & Gusev, 1990; Hoshiba, 1995; Gusev & Abubakirov, 1996; Jing et al., 2014; Sato & Emoto, 2018), velocity and heterogeneity stratification (Hoshiba, 1997; Margerin et al., 1998; Yoshimoto, 2000), coupling between shear and compressional waves (Margerin et al., 2000; Przybilla et al., 2009), laterally varying velocity and scattering structures (Sanborn et al., 2017). With the growth of computational power, Monte-Carlo simulations opened up new venues for the application of RT in seismology: imaging of deep Earth heterogeneity (Margerin & Nolet, 2003; Shearer & Earle, 2004; Mancinelli & Shearer, 2013; Mancinelli et al., 2016b), mapping of the depth-dependent scattering and absorption structure of the lithosphere (Mancinelli et al., 2016a; Takeuchi et al., 2017), modeling of propagation anomalies in the crust (Sens‐-Schönfelder et al., 2009; Sanborn & Cormier, 2018), P-to-S conversions in the teleseismic coda (Gaebler et al., 2015), to cite a few examples only.

Recently, RT has also been applied to the computation of sensitivity kernels for time-lapse imaging methods such as coda wave interferometry (see e.g. Poupinet et al., 1984; Snieder, 2006; Poupinet et al., 2008). Coda Wave Interferometry (CWI) exploits tiny changes of waveforms in the coda to map the temporal variations of seismic properties in 3-D. The mapping relies on the key concept of sensitivity kernels, which, in the framework of CWI, were introduced by Pacheco & Snieder (2005) in the diffusion regime and Pacheco & Snieder (2006) in the single-scattering regime. These kernels take the form of spatio-temporal convolutions of the mean intensity in the coda. It was later pointed out by Margerin et al. (2016) that an accurate computation of traveltime sensitivity kernels, valid for an arbitrary scattering order and an arbitrary spatial position, requires the knowledge of the angular distribution of energy fluxes in the coda. These fluxes, or specific intensities, are directly predicted by the radiative transfer model, which makes it attractive for imaging applications.

In noise-based monitoring (Wegler & Sens-Schoenfelder, 2007) -also known as Passive Image Interferometry (PII)- the virtual sources and receivers are located at the surface of the medium so that the early coda is expected to contain a significant proportion of Rayleigh waves. At longer lapse-time, the surface waves couple with body waves and the coda eventually reaches an equipartition regime when all the propagative surface and body wave modes are excited to equal energy (Weaver, 1982; Hennino et al., 2001). Because the volumes explored by surface and body waves are significantly different, the knowledge of the composition of the coda wavefield at a given lapse-time in the coda is key to locate accurately the changes at depth in the crust.

Obermann et al. (2013, 2016) proposed to express the sensitivity of coda waves as a linear combination of the sensitivity of surface and body waves, whose kernels are computed from scalar RT theory in 2-D and 3-D, respectively. The relative contribution of the 2-D and 3-D sensitivity kernels at a given lapse-time in the coda is determined by fitting the traveltime shift predicted by the theory against full wavefield numerical simulations in scattering media, where the background seismic velocity is perturbed in a fine layer at a given depth. This method has been validated in the case of 1-D perturbations through numerical tests and has the advantage of modeling exactly the complex coupling between surface and body waves in heterogeneous media. Furthermore, it can easily incorporate realistic topographies, which is important for the monitoring of volcanoes. The main drawbacks of the approach of Obermann et al. (2016) are the numerical cost and the fact that it requires a good knowledge of the scattering properties of the medium, which have to be determined by other methods such as MLTWA (Fehler et al., 1992; Hoshiba, 1993).

This brief overview illustrates that CWI and PII would benefit from a formulation of RT theory which incorporates the coupling between surface and body waves in a self-consistent way. In the case of a slab bounded by two free surfaces, Trégourès & Van Tiggelen (2002) derived from first principles a quasi 2-D RT equation where the wavefield is expanded onto a basis of Rayleigh, Lamb and Love eigenmodes. Thanks to the normal mode decomposition, this model incorporates the boundary conditions at the level of the wave equation. The energy exchange between surface and body waves is treated by normal mode coupling in the Born approximation. A notable advantage of this formulation is the capacity to predict directly the energy decay in the coda and its parttioning onto different components. The two main limitations for seismological applications are the slab geometry, which may not always be realistic and the fact that the disorder should be weak, i.e., the mean free time should be large compared to the vertical transit time of the waves through the slab.

Zeng (2006) proposed a system of coupled integral equations to describe the exchange of energy between surface waves and body waves in the seismic coda. The formalism used by the author is interesting and bears some similarities with the one developed in this work, although we formulate the theory in integro-differential form. A blind spot in the work of Zeng (2006) is the coupling between surface and body waves, which is introduced in an entirely phenomenological way and differs significantly from our findings. A very promising investigation of the energy exchange between surface and body waves on the basis of the elastodynamic equations in a half-space geometry was performed by Maeda et al. (2008). Using the Born approximation, these authors calculated the scattering coefficients between all possible modes of propagation in a medium containing random inhomogeneities. The main limitation of their theory comes from the fact that the conversion from body to surface waves is quantified by a non-dimensional coefficient, which makes it difficult to extend their results beyond the single-scattering regime. The authors argue that the absence of a characteristic scale-length for body-to-surface attenuation is a consequence of the fact that all conversions occur in approximately one Rayleigh wavelength in the vicinity of the surface.

In this work, we revisit the problem of coupling surface and body waves in a RT framework using an approach similar to the one of Maeda et al. (2008). For simplicity, we limit our investigations to a scalar model based on the Helmholtz equation with an impedance (or mixed) boundary condition in a half-space geometry. To make the presentation self-contained, we review the most important features of this particular wave equation. Specifically, we recall that the modes of propagation are composed of body waves, and surface waves whose penetration depth depends on the impedance condition only. Hence, our model mimics the situation on Earth while minimizing the mathematical complexity. We then introduce a simple point-scattering model and study its properties in the Born approximation. Using the additional assumption that the surface reflexion randomizes the phase of the reflected wave, we are able to derive simple expressions for the scattering mean free path of both body and surface waves including the coupling between the two. We elaborate on this result to establish a set of two coupled RT equations satisfied by the specific energy density of surface and body waves using a phenomenological approach. Some consequences of our simple theory are explored, in particular the establishment of a diffusion and equipartition regime. Monte-Carlo simulations show the potential of the approach to model the transport of energy in the seismic coda from single-scattering to diffusion.

2 Scalar wave equation model with surface and body waves

In this section, we present the basic ingredients of our scalar model based on the Helmholtz equation. We describe how an appropriate modification of boundary conditions at the surface of a half-space gives rise to the presence of a surface wave. We subsequently present an expression of the Green’s function and its asymptotic approximation. The concept of density of states, important for later developments, is recalled. For a thorough treatment of the mathematical foundations of our model, the interested reader is refered to the monograph of Hein Hoernig (2010).

2.1 Equation of motion

We consider a 3-D version of the membrane vibration equation in a half-space geometry:

[TABLE]

where $t$ is time and $\mathbf{R}$ is the position vector. It may be further decomposed as $\mathbf{R}=\mathbf{r}+z\mathbf{\hat{z}}$ ($z\geq 0$) with $\mathbf{r}=x\mathbf{\hat{x}}+y\mathbf{\hat{y}}$ and $(\mathbf{\hat{x}},\mathbf{\hat{y}},\mathbf{\hat{z}})$ denotes a Cartesian system. In Eq. (1) $u$, $\rho$ and $T$ may be thought of as the displacement, the density and the elastic constant of the medium, respectively. The wave Eq. (1) is supplemented with the boundary conditions:

[TABLE]

The case of interest to us corresponds to $\alpha>0$, i.e. when, as recalled below, the boundary can support a surface wave. Equation (1) can be derived by applying Hamilton’s principle to the following Lagrangian density:

[TABLE]

Thanks to the last term of the Lagrangian (3), which corresponds to a negative elastic potential energy stored at the surface $z=0$, the first B.C. in Eq. (2) becomes natural in the sense of variational principles. To make the presentation self-contained, we explain in Appendix A the origin of the delta function in Eq. (3) in the simple case of a finite string with mixed boundary conditions at one end and free boundary conditions at the other end.

In the case of a harmonic time dependence $u\propto e^{-i\omega t}$ ($\omega>0$), the vibrations are governed by Helmholtz Eq.:

[TABLE]

with $c=\sqrt{T/\rho}$ the speed of propagation of the waves in the bulk of the medium. Eq. (4) is complemented with the mixed boundary condition $\partial_{z}u+\alpha u=0$ at $z=0$ and an outgoing wave condition at infinity. From (3), we can deduce the energy flux density vector $\mathbf{J}$ and the energy density $w$ using the concept of stress-energy tensor (Morse & Ingard, 1986). For harmonic motions, their average value over a period can be expressed as:

[TABLE]

In the following section, we recall the consequences of mixed boundary conditions on the Helmoltz Eq., in particular the fact that it gives rise to a surface wave mode.

2.2 Eigenfunctions and Green’s function

Due to the translational invariance of the medium, we look for eigen-solutions of Eq. (4) in the form $u=\psi(z)e^{i\mathbf{k}_{\parallel}\cdot\mathbf{r}}$ with $\mathbf{k}_{\parallel}=(k_{x},k_{y},0)$. This leads to a self-adjoint eigenvalue problem in the $z$ variable only. For $\alpha>0$, part of the spectrum is discrete with eigenfunction :

[TABLE]

with $\mathbf{k}_{\parallel}\cdot\mathbf{k}_{\parallel}-\alpha^{2}=\dfrac{\omega^{2}}{c^{2}}$. The rest of the spectrum forms a continuum of body waves with normalized eigenfunctions:

[TABLE]

with $q^{2}+\mathbf{k}_{\parallel}\cdot\mathbf{k}_{\parallel}=\dfrac{\omega^{2}}{c^{2}}$ and:

[TABLE]

Note the relations: (1) $q=(\omega\cos j)/c=k\cos j$ with $j$ the incidence angle of the body wave and (2) $|r(q)|^{2}=1$, i.e., there is total reflection at the surface. For later reference, we introduce a specific notation for the vertical eigenfunction of body waves:

[TABLE]

with $\mathbf{\hat{n}}\cdot\mathbf{\hat{z}}=\cos j$. Note that throughout the paper, we use a hat to denote a unit vector. The surface waves (7) and body waves (8) are normalized and orthogonal in the sense of the scalar product $\langle u|v\rangle=\int_{\mathbb{R}^{3}_{+}}u(\mathbf{R})^{*}v(\mathbf{R})d^{3}R$, where $u$ and $v$ are arbitrary square integrable functions. The surface wave phase velocity $c_{\phi}$ is given by:

[TABLE]

and is always smaller than the bulk wave velocity $c$. The group velocity may be obtained in two different manners, namely, (1) using the classical formula based on the interference of a wave packet:

[TABLE]

and (2) using the principle of energy conservation:

[TABLE]

where the brackets indicate an integration over the whole depth range.

The dispersion properties of the surface wave in our scalar model are illustrated in Figure 1. As is evident from Eqs (11)-(12), the group velocity is always faster than both the phase and body wave velocity. In the high-frequency limit, the phase and group velocity tend to the common value $c$. Using definition (13) it is possible to define the energy velocity of a body wave eigenmode (see Eq. 8):

[TABLE]

where $\langle\rangle_{h}$ denotes an integration from the surface to depth $h$. The depth averaging smoothes out the oscillations of $\mathbf{J}$ and $w$ caused by the interference between the incident and reflected amplitudes. The passage to the limit is necessary because the integrals over depth diverge. Eq. (14) can be interpreted as follows. In the full-space case, the current vector of a single unit-amplitude plane wave with wavevector $\mathbf{k}=k(\cos j\mathbf{\hat{z}}+\sin j\mathbf{\hat{k}}_{\parallel})$ is given by $\mathbf{J}=\rho\omega c^{2}\mathbf{k}/2$ and carries an energy density $w=\rho\omega^{2}/2$. If we define $\mathbf{k}_{r}$ as the mirror image of $\mathbf{k}$ across the plane $z=0$ and consider the sum of the current vector of two plane waves with wavevectors $\mathbf{k}$ and $\mathbf{k}_{r}$ we obtain $\mathbf{J}+\mathbf{J}_{r}=\rho\omega^{2}c\sin j\mathbf{\hat{k}}_{\parallel}$. After normalization by the sum of energy densities, the result (14) is recovered. In other words, on average, the energy transported by a body wave mode is simply the sum of the energies transported by the incident and reflected waves, as if the two were independent.

Using the eigenmodes (7) and (8), one may obtain an exact representation of the Green’s function of Helmholtz Eq. with mixed BC in the form:

[TABLE]

where $z_{0}$ denotes the source depth, $\epsilon$ is a small positive number which guarantees the convergence of the integrals and the star ∗ denotes complex conjugation. In Eq. (15) the first (resp. second) line represents the body wave (resp. surface wave) contribution. As shown in Appendix B, the surface wave term can be computed analytically in terms of Hankel functions. The following far-field approximation of the Green’s function of the Helmholtz Eq. (4) can be obtained using the stationary phase approximation for the body wave term:

[TABLE]

with $k_{s}=\omega/c_{\phi}$, $\mathbf{\hat{R}}=\mathbf{R}/R$ and $\mathbf{R}=\mathbf{r}+z\mathbf{\hat{z}}$. The expansion (16) is performed with respect to the midpoint of the source point and its mirror image by the surface $z=0$. The $z$ dependence of the first term is simply given by the body wave eigenfunction (8). For further computational details, the reader may consult Appendix B.

2.3 Source radiation and density of states

We now compute the energy radiated by a point source located at $(0,0,z_{0})$. To do so, we introduce a cylindrical surface $C$ of radius $R_{c}$ extending from the free surface to a depth $h$ greater than $z_{0}$ and large compared to $1/\alpha$. We close this surface with a hemispherical cap $H$ of radius $R$ centered at the surface point $(0,0,0)$. The geometry is schematically depicted in Figure 2. The energy flux vector (5) of the radiated field contains terms that are purely surface, purely bulk and cross-terms. The contribution of surface waves to the flux across the hemispherical surface is negligible (the error made is exponentially small). The contribution of body waves to the flux across the lateral cylindrical surface is also negligible because this surface subtends a solid angle which goes to [math] as $R_{c}$ increases. The cross-terms are negligible across the whole surface because the coupled surface/body wave term decays algebraically faster than the surface wave term on the cylindrical surface and exponentially faster than the body wave term on the hemispherical surface. Hence, we may split the flux of radiated waves into a contribution of surface and body waves, respectively.

The energy transported per unit time by body waves through the hemispherical cap $H$ is given by:

[TABLE]

with $G_{b}$ the body wave part of Green’s function. In the second line of Eq. (17), the integral is over the space directions subtended by the hemispherical cap. (N.B.: strictly speaking, the total solid angle is not equal to $2\pi$ because one should remove the directions corresponding to the cylinder. But as noted before, the measure of this set of directions goes to zero as $R_{c}$ goes to infinity.) The function defined in Eq. (17) oscillates with depth around the following mean value:

[TABLE]

Here the brackets may have at least 2 different meanings. The most obvious is an average over depth, as was done in the calculation of the group velocity. But we may also assume that the surface “scrambles” the phase of the reflected wave $\phi_{r}$ so that it becomes a random variable. In this scenario, the brackets would mean an average over all realizations of the random reflection process. Upon averaging over phase or depth, the interference pattern between the incident and reflected wave is smoothed out, so that the two approaches yield the same result. Note that the randomization of the phase does not affect energy conservation because the incident flux is still totally reflected. In particular, the discussion following the interpretation of Eq. (14) would still be valid. In practice, the assumption that the phase of the reflected wave is randomized by the surface may not be as unrealistic as it seems. Observations of reflected $SH$ waves by Kinoshita (1993) at borehole stations in Japan indeed suggest that the reflected field is a distorted version of the incident one. This concurs with the general view that the subsurface of the Earth is highly heterogeneous at scales that can be much smaller than the wavelength and brings support to the idea that aberrating fine layers could indeed randomize the phase of the reflected wave as we hypothesize. In what follows, the scrambled-phase assumption will be adopted to simplify the treatment of the reflection of body waves at the surface.

The energy transported per unit time by surface waves through the lateral cylindrical surface $C$ is given by:

[TABLE]

with $G_{s}$ the surface wave part of Green’s function. Because $h$ is large compared to $1/\alpha$ the integral over depth may be performed from [math] to $+\infty$ (the error incurred is exponentially small). We find the depth dependent surface-to-body energy ratio:

[TABLE]

Formula (20) can also be understood in the light of the local density of states $n_{s,b}$ defined as (Sheng, 2006):

[TABLE]

where $k_{s,b}(\omega)$ stands for the wavenumber of surface or body waves. Using the spectral representation (15), one obtains the (exact) formulas:

[TABLE]

which show that the partitioning of the energy radiated into surface and body waves by the source $R(z_{0})$ is given by the ratio of their local density of states $n_{s}(z_{0})/\langle n_{b}\rangle$. Finally, we may compute the partitioning ratio $\mathcal{R}$ between the modal density of surface and body waves by integrating Eq.(24) over $z$ and taking the ratio with (23). This yields the simple result:

[TABLE]

where it is to be noted that the modal density ratio $\mathcal{R}$ is independent of the scale length $\alpha$ appearing in the mixed boundary condition of the Helmholtz equation. This result could have been deduced directly from the dispersion relations of body and surface waves using classical mode counting arguments (Kittel et al., 1976). It is worth noting that the local density of states (23) is exactly the same as in the case of the Helmholtz equation in full 3-D space. Although the integral in (22) is carried over one hemisphere only, each eigenmode $\psi_{b}$ is composed of an incident and a reflected wave, which -on average- doubles its contribution compared to a single plane wave state. In the next section, we use our knowledge of the Green’s function to derive the scattering properties of surface and body waves including the coupling between the two modes of propagation.

3 Single scattering by a point scatterer

In this section, we calculate the energy radiated by a single scatterer in a half-space geometry for incident surface or body waves. For simplicity, we restrict our investigations to point scatterers and employ Born’s approximation. The resulting expressions are simplified following the scrambled phase approximation and interpreted in terms of scattering cross-sections.

3.1 Scattering of a surface wave

We now consider the following perturbed Helmholtz Eq.:

[TABLE]

Here $a$ represents the typical linear dimension of the scatterer located at $(0,0,z_{s})$ and it is understood that $ka\ll 1$. $\epsilon$ is the local perturbation of inverse squared velocity. Following the standard procedure (Snieder, 1986), we look for solutions of Eq. (26) of the form: $u=u_{0}+u_{s}$, where $u_{0}(\mathbf{r},z)=e^{-\alpha z+i\sqrt{\alpha^{2}+\omega^{2}/c^{2}}x}$ is a surface wave eigenmode of the Helmholtz equation (the incident field) and $u_{s}$ is the scattered field. Using the Born approximation, one obtains:

[TABLE]

Introducing the coupling strength $S=k^{2}a^{3}\epsilon e^{-\alpha z_{s}}$, one may express the energy radiated by the body waves through the hemispherical cap $H$ (per unit time) as:

[TABLE]

Note that the coupling term $S$ depends on both intrinsic properties of the scatterer -size and strength of perturbations- as well as on the properties of the incident wave -depth dependence of eigenfunction and frequency-.

For a unit amplitude surface wave, the vertically-integrated time-averaged energy flux density is given by

[TABLE]

The ratio of (28) and (29) gives the surface-to-body scattering cross-section:

[TABLE]

This cross-section has unit of length. In the case where the surface scrambles the phase of the reflected wave, we may compute the mean conversion scattering cross-section by taking the average over the random phase $\phi_{r}$.

[TABLE]

Let us remark again that the averaging smoothes out the interference pattern of the body wave eigenfunction $\psi_{b}$ but does not affect the conservation of energy. Furthermore, the averaging procedure makes the scattering pattern isotropic since $\langle|\psi_{b}(\mathbf{\hat{R}},z_{s})|^{2}\rangle=2$ with an equal contribution of upgoing and downgoing waves. The computation of the surface-to-surface scattering cross-section proceeds in a similar way. The surface-wave energy radiated through the lateral surface is given by:

[TABLE]

Normalizing the result (32) by the incident flux yields the surface-to-surface scattering cross section:

[TABLE]

again with unit of length. If we have a collection of point scatterers with volume density $n$, we may define a surface wave scattering mean free path using the independent scattering approximation as (Lagendijk & Van Tiggelen, 1996; Trégourès & Van Tiggelen, 2002; Maeda et al., 2008):

[TABLE]

The approximate formula (34) neglects all recurrent interactions between the scatterers, which is valid for sufficiently low concentrations of inclusions.

3.2 Scattering of body waves

In the case of an incident body wave mode (8), the computation of the energy radiated in the form of body or surface waves can be performed as in the previous section. The definition of the scattering cross-section, however, is problematic if we stick to the modal description. Indeed, the vertically integrated energy flux density of a body wave mode, as defined in Eq. (8), diverges. We must therefore come back to a conventional plane wave description. We first consider the situation where the scatterer is located at a large depth in the half-space. In this case, it appears reasonable to think that the scattering of a body wave mode should be equivalent to the scattering of a plane wave in the full-space case, at least in a sense to be specified below. To verify the correctness of this assertion, we start by computing the body wave cross-section in absence of a boundary. Assuming a unit amplitude incident plane wave and using again the Born approximation, one may express the scattered energy as:

[TABLE]

Normalizing the result by the energy flux density of the incident wave:

[TABLE]

we obtain:

[TABLE]

In the half-space geometry, the incident wave has the form (8). The total energy scattered in the form of body waves is still given by Eq. (28) provided one redefines the coupling constant as $S=k^{2}a^{3}\epsilon|\psi_{b}(\mathbf{\hat{R}}_{i},z_{s})|^{2}$, where $\mathbf{\hat{R}}_{i}$ refers to the incidence direction of the body wave. Eq. (28) differs from formula (35) as a consequence of the interference between the incident and reflected waves. In the case of a scattering medium, we may expect these interferences to be blurred due to the randomization of the phase by the scattering events. In this scenario, the phase of the incident and reflected waves may be expected to be uncorrelated. Upon averaging the result (28) over the random phase of the reflected wave, one obtains:

[TABLE]

which is exactly the double of the full-space result (35). Keeping in mind that the energy flux of the incident plane wave interacts twice with the scatterer (direct interaction + interaction after reflexion) and using the assumption that the energy fluxes of incident and reflected waves do not interfere and may therefore be added, we obtain the result:

[TABLE]

As announced, the average scattering cross-section of body waves in the half space is the same as in the full space for a scatterer located far away from the boundary. Conceptually, the “scrambled phase” approximation allows us to extend the result (39) to a scatterer located at an arbitrary depth in the medium thanks again to the assumption that the incident and reflected waves are statistically independent. In this scenario, on average, the surface does not modify anything to the scattering of body waves into body waves as compared to the full-space case.

Following the same approach and approximation, we can calculate the body-to-surface scattering cross-section. The energy radiated in the form of surface waves is given by:

[TABLE]

After averaging and normalization by the total energy flux (incident + reflected), one finds:

[TABLE]

Note that all the remarks pertaining to the mean scattering pattern made after Eq. (30) also apply to the derivation of Eq. (39) and (41). The scattering cross-sections $\sigma^{b\rightarrow b}$ and $\sigma^{b\rightarrow s}$ have unit of surface. Using the independent scattering approximation again, we conclude that the inverse scattering mean free path of body waves defined as:

[TABLE]

depends on the depth in the medium, as a consequence of the coupling with surface waves. A simple but fundamental reciprocity relation may be established between the surface-to-body and body-to-surface scattering mean free time. The later may be expressed as:

[TABLE]

while the former may be obtained after performing the integral over depth in Eq. (34):

[TABLE]

The ratio between the two quantities is given by:

[TABLE]

where Eq. (20) has been used. Eq. (45) establishes a link between the surface-to-body versus body-to-surface conversion rates and the local density of states. In the case of elastic waves, a similar relation applies (Weaver, 1990; Ryzhik et al., 1996): the ratio between the $P$-to-$S$ and $S$-to-$P$ scattering mean free times is given by the ratio of the density of states of $P$ and $S$ waves. As shown in the next section, the relation (45) plays a key role in the establishment of an equipartion between surface and body waves.

4 Equation of radiative transfer

In this section, we employ standard energy balance arguments to derive a set of coupled equations of RT for surface and body waves in a half-space containing a uniform distribution of point-scatterers. A notable feature of our formulation is the appearance of the penetration depth of surface waves as a parameter in the Equations.

4.1 Phenomenological derivation

Before we establish the transport equation, a few remarks are in order. It is clear that the transport of surface wave energy is naturally described by a specific surface energy density $\epsilon_{s}(t,\mathbf{r},\mathbf{\hat{n}})$ where $\mathbf{\hat{n}}$ is a unit vector in the horizontal plane. The surface energy density of surface waves may in turn be defined as:

[TABLE]

where the integral is carried over all propagation directions in the horizontal plane. Note that the $\epsilon_{s}$ symbol should not be confused with the strength of fluctuations in Eq. (26). In contrast with surface waves, the transport of body waves is described by a specific volumetric energy density $e_{b}(t,\mathbf{r},\mathbf{\hat{k}},z)$ where $\mathbf{\hat{k}}$ is a vector on the unit sphere in 3-D. The volumetric energy density of body waves is again obtained by integration of the specific energy density over all propagation directions in 3-D:

[TABLE]

In the usual formulation of transport equations, energy densities have the same unit (either surfacic or volumetric). In order to treat on the same footing surface and body waves, we introduce the following volumetric energy density of surface waves as:

[TABLE]

It is clear that upon integration of $e_{s}$ over depth, one recovers the surface density $\epsilon_{s}$. The exponential decay of the surface wave energy density is directly inherited from the modal shape and implies that the coupling between surface and body waves mostly occurs within a skin layer of typical thickness $1/2\alpha$. For future reference, we introduce the following notation:

[TABLE]

to represent the volumetric energy density of surface waves, consistent with Eq. (47). Note that the total energy density at a given point will be defined as the sum $(E_{b}(t,\mathbf{r},z)+E_{s}(t,\mathbf{r},z))$, thereby implying the incoherence between the two types of waves.

With this definition of the surface wave energy density, the phenomenological derivation of the radiative transport equation follows exactly the same procedure as in the multi-modal case (Turner & Weaver, 1994). A beam of energy followed along its path around direction $\mathbf{\hat{n}}^{i}$ is affected by (1) conversion of energy into other propagating modes and/or deflection into other propagation directions $\mathbf{\hat{n}}^{o}\neq\mathbf{\hat{n}}^{i}$; (2) a gain of energy thanks to the reciprocal process: a wave with mode $i$ propagating in direction $\mathbf{\hat{n}}^{i}$ can be converted into a wave with mode $o$ propagating in direction $\mathbf{\hat{n}}^{o}$ by scattering. In the general case of finite size scatterers we may anticipate scattering to be anisotropic. To describe such an angular dependence of the scattering process, we may introduce normalized phase functions $p^{i\rightarrow o}(\mathbf{\hat{n}}^{o},\mathbf{\hat{n}}^{i})$. The phase function may be understood as the probability that a wave of mode $i$ propagating in direction $\mathbf{\hat{n}}^{i}$ be converted into a wave of mode $o$ propagating in direction $\mathbf{\hat{n}}^{o}$. To be interpretable probabilistically, its integral over all outgoing directions ($\mathbf{\hat{n}}^{o}$) should equal 1. Note that in addition to the incoming and outgoing propagation directions, the phase function may also depend on the depth in the medium as a consequence of the presence of the surface. In the case of point scatterers, this complexity disappears thanks to the scrambled phase assumption and will therefore not be considered in our formalism.

A detailed local balance of energy yields the following system of coupled transport equations:

[TABLE]

where the terms $s_{s,b}$ represent sources of surface and body waves. To take into account the reflection of body waves at the surface, the sytem (50) is supplemented with the following boundary condition for the energy density $e_{b}$:

[TABLE]

where $\mathbf{\hat{k}}_{r}$ is the mirror image of the incident direction $\mathbf{\hat{k}}$ $(\mathbf{\hat{k}}\cdot\mathbf{\hat{z}}<0)$ across the horizontal plane $z=0$. The boundary condition (51) is compatible with the assumption that the incident and reflected waves are statistically independent.

In the simple case of a unit point-like source at depth $z_{0}$, the terms $s_{s}$ and $s_{b}$ are given respectively by:

[TABLE]

where $R(z_{0})$ is the energy partitioning ratio defined in Eq. (20). Note that $s_{s}$ follows the same exponential decay as $e_{s}$ with depth $z$ (see Eq. 48), which takes into account the vertical dependence of the surface wave eigenfunction. In Eq. (52), the complex dependence of the body wave radiation with depth has been simplified by averaging the exact result (17) over the random phase of the reflected wave. As a consequence, the energy radiation at the source covers uniformly the whole sphere of space directions in 3-D.

A similar remark applies to the scattering from body waves to body waves, which have been treated as if the surface was absent in Eq. (50). Again, this approximation is admissible if the scattered upgoing and downgoing energy fluxes are statistically independent, a condition which is guaranteed by the randomization of the phase of the waves upon reflection at the surface in our model. As discussed in the previous section, in the case of point scatterers the phase averaging procedure makes all scattering processes isotropic which allows us to simplify the system of Eqs (50) by evaluating the scattering integrals on the right-hand side of Eq. (50):

[TABLE]

where we have used the normalization condition of the phase functions. The decreasing efficacy of scattering conversions with depth is guaranteed by the exponential decay of $\tau^{b\rightarrow s}(z)^{-1}$ and $E_{s}(t,\mathbf{r},z)$ in the first and second Eq. of the coupled sytem (53), respectively. It is worth recalling that the depth dependence of the scattering mean free times of body waves $\tau^{b}(z)$ and $\tau^{b\rightarrow s}(z)$ is caused by the coupling with surface waves and not by a stratification of heterogeneity.

4.2 Energy conservation and equipartition

A self-consistent formulation of coupled transport equations should verify two elementary principles: energy conservation and equipartition. To demonstrate these properties from the basic set of equations (50) or (53), we proceed in the usual fashion (Turner & Weaver, 1994). An integration of each equation over all possible propagation directions $\mathbf{\hat{n}}$ (in 2-D) or $\mathbf{\hat{k}}$ (in 3-D) yields:

[TABLE]

where we have used the definition of the mean free times on the RHS of (54). Note that we have dropped the source terms since they are not essential to our argumentation. In Eq. (54), we have introduced the energy flux density vector of surface and body waves:

[TABLE]

Note that $\mathbf{J}_{s}$ is contained in a horizontal plane. Upon integration over the whole half-space, the terms which contain the current density vectors can be converted into surface integrals that vanish. Denoting by a double bar an integration over $\mathbf{r}$ and $z$, and summing the two Eqs of the system (54) leaves us with:

[TABLE]

which demonstrates the conservation of energy.

In order to prove the existence of equipartition, we integrate each Eq. of the system (54) over the horizontal plane and from the surface to a finite depth $h$, typically large compared to the surface wave penetration depth. To simplify the derivation, we assume that at $z=h$ lies a perfectly reflecting surface through which no energy can flow. Because $\alpha h\gg 1$, the medium may be still be considered as a half-space in the treatment of surface waves. This leads us to:

[TABLE]

where the single bar denotes an integration over $\mathbf{r}$ only. Note that $\bar{\bar{E}}_{s}$ is the total energy of surface waves (up to an exponentially small correction term), in contrast with $\bar{E}_{b}$ which represents the body wave energy per unit depth. Our goal is not to solve the system of Eq. (57) in its full generality but rather to exhibit an equipartition solution. In this regime, we expect the distribution of body wave energy to become independent of depth. Hence we look for solutions of the system (57) of the form:

[TABLE]

Reporting the ansatz (58) into the integro-differential system (57), one arrives at a linear and homogeneous system of algebraic equations. A non-zero solution is obtained only if the determinant of the sytem vanishes which implies:

[TABLE]

One of the solutions is given by $\lambda=0$ which corresponds to the asymptotic equipartition state such as:

[TABLE]

where Eq. (25) and (45) have been used, and the depth integral is again extended to $+\infty$ thanks to the assumption $\alpha h\gg 1$. Eq. (60) illustrates that the ratio between the energy density of surface and body waves approaches the ratio of their density of states $\mathcal{R}$ at long lapse-time. Note that the energy density is not perfectly homogenized spatially, even at equipartition, as a consequence of the decay of the surface wave eigenfunction with depth. This may be related to the depth-dependence of the density of states near the boundary (see e.g. Hennino et al., 2001). The result (60) could have been predicted using the usual concept of equipartition which states that when filtered around a narrow frequency band, all the propagating modes of a diffuse field should be excited to equal energy (Weaver, 1982). In the case where the medium is unbounded at depth, there will be a flux of energy across the lower boundary $z=h$ which vanishes as the lapse-time increases. As demonstrated through numerical simulations later in this paper, equipartition also sets in in this configuration, though probably more slowly than in the slab geometry.

4.3 Diffusion Approximation

Having established the existence of an equipartition state, we now derive a diffusion approximation for the transport process. At the outset, it should be clear that the volumetric energy density $E_{s}$ may not be the solution of a diffusion equation because it exhibits a decay with depth which is inherited from the modal shape and therefore independent of the scattering properties. To circumvent the difficulty, we derive a closed diffusion equation for $E_{b}$ from which we subsequently deduce the energy density of surface waves. To simplify the calculations, we assume that the scattering is isotropic. Proceeding in the usual fashion, we expand the specific energy density into its first two angular moments (Akkermans & Montambaux, 2007):

[TABLE]

Multiplying, respectively, the first and second line of Eq. (53) by the unit vectors $\mathbf{\hat{n}}$ and $\mathbf{\hat{k}}$, integrating over all possible directions and employing the moment expansion (61), we obtain the following set of Equations:

[TABLE]

where $\nabla_{\parallel}$ denotes the gradient operator in the horizontal plane. Note that the expansions (61) are used only to evaluate the integral of the gradient on the left-hand side of the RT Equation. All other terms follow directly either from the definition of the energy flux density vector or the assumption of isotropic scattering. Our final approximation consists in neglecting the derivative of the current vector with respect to time which yields the equivalent of Ohm’s law for the mutiply-scattered waves:

[TABLE]

where the following notations have been introduced:

[TABLE]

Note that the total reflection condition (51) imposes that there is no net flux of body waves across the surface $z=0$. To make progress, we now invoke the equipartition principle to fix the ratio between the energy densities of surface and body waves:

[TABLE]

in agreement with Eq. (60). This allows us to express the total energy flux as:

[TABLE]

which may be interpreted as a generalization of Ohm’s law for diffuse waves. Eq (66) demonstrates that the coupling between surface and body waves renders the energy transport both depth-dependent and anisotropic. It is worth emphasizing that depth-dependence and anisotropy are caused neither by specific orientations/shapes of the scatterers nor by the non-homogeneity of the statistical properties. As further discussed below, these properties stem from the coupling between body and surface wave modes.

The conservation equation for the total energy $E=E_{s}+E_{b}$ excited by a point source at $t=0$ is obtained by taking the sum of the set of Eqs (54):

[TABLE]

where the right-hand side now contains the source term with $z_{0}$ the source depth. Making use of Eqs (65) and (63), we obtain the following diffusion-like equation verified by the body wave energy density:

[TABLE]

The last term on the left-hand side of Eq. (68) differs from the traditional form for the diffusion model due to the $(1+R(z))^{-1}$ factor in front of the derivative operators. Actually, this difference is purely formal as may be shown by the change of variable $z\rightarrow z^{\prime}$ where $z^{\prime}$ is defined as:

[TABLE]

where Eq. (20) has been used. In the new variables, Eq. (68) may be rewritten as:

[TABLE]

In Eq. (70), we have introduced the notations $E^{\prime}_{b}(t,\mathbf{r},z^{\prime})=E_{b}(t,\mathbf{r},z)$, $z^{\prime}_{0}=\int_{0}^{z_{0}}(1+R(x))dx$, as well as the following definitions of the horizontal and vertical diffusivities:

[TABLE]

Hence, a simple change of scale in the vertical direction reduces Eq.(68) to the conventional diffusion Eq. (70). Note that the $(1+R(z_{0}))^{-1}$ factor on the right-hand side has been absorbed by the change of variable (69).

We explore the consequences of Eq. (71) by first considering the case $z^{\prime}\to\infty$. According to Eq. (69), this implies $z^{\prime}\approx z$. Since the partitioning ratio $R(z)$ goes to zero at large depth in the medium, Eq. (71) indicates that the vertical and horizontal diffusivities become equal and the diffusion tensor isotropic. Furthermore, because the coupling between surface and body waves is negligible at large depth, its magnitude tends to the constant value $D_{b}^{bulk}=c^{2}\tau^{b\rightarrow b}/3$, as expected on physical grounds. In other words, the diffusion process at depth is governed by a simple 3-D diffusion equation for body waves with diffusion constant $D_{b}^{bulk}$. This in turn suggests that at long lapse-time, the coda should decay as $t^{-3/2}$ in a non-absorbing medium. This point will be further substantiated by numerical simulations.

In the vicinity of the surface $z=O(\alpha^{-1})$, Eq. (71) shows that the diffusivity of coupled body and surface waves is both depth dependent and non-isotropic. The origin of the $z$-dependence is clear since the efficacy of the coupling between surface and body waves decays exponentially with depth. In the vicinity of the surface, the anisotropy stems from the transport of a fraction of the energy by surface waves whose velocity and scattering mean free time differ from the one of body waves. In Eq. (71) the transverse diffusivity is recognized as a weighted average of the surface and body wave diffusivities with coefficients dictated by the equipartition principle. The vertical diffusivity is -up to the $(1+R(z))$ pre-factor inherited from the change of scale in the vertical direction- equal to the diffusivity of body waves. In the next section, we illustrate the transport process of coupled body and surface waves by numerically simulating the system of Eq. (50).

5 Monte-Carlo Simulations

In this section, we explore some of the key features of our model with the aid of numerical simulations. The approach to equipartition as well as the role of mode coupling in the coda excitation are illustrated.

5.1 Overview of the method

As outlined in introduction, Monte-Carlo simulations have been used for more than thirty years in seismology to simulate the transport of seismic energy in heterogeneous media. Our approach to the solution of the coupled set of transport equations (53) for surface and body waves follows closely the approach of Margerin et al. (2000), with some appropriate modifications which we outline briefly.

Energy transport is modeled by the simulation of a large number of random trajectories of particles or seismic phonons (Shearer & Earle, 2004). Each particle is described by its mode, position, propagation direction and time. The initial mode of propagation is randomly selected, following the source energy partitioning ratio (20), and the initial propagation direction is a uniformly distributed random vector in 2-D (resp. 3-D) for surface (resp. body) waves. Note that when the particle is of surface type, the particle propagates in a horizontal plane and its exact depth is immaterial. In fact, we may say that a particle of surface type is present at all depth with a probability distribution given by $p_{s}(z)=2\alpha\exp(-2\alpha z)$ inherited from the modal shape. The lapse-time to the first scattering event is randomly determined and obeys a simple exponential distribution when the particle represents surface waves. In the case of body waves, the selection process is more complicated because their scattering mean free time depends on the depth in the medium. To address this difficulty, we employ the method of delta collisions, which simulates in a simple and exact way a completely general distribution of scattering mean free time. We will not detail the method here and refer the interested reader to the pedagogical treatment by Lux & Koblinger (1991). At each scattering event, the mode of the particle is randomly selected by interpreting probabilistically Eqs (34) and (42) defining the scattering attenuations. As an example, $(1/l^{s\rightarrow b})/(1/l^{s})$ may be interpreted as the transition probability from a surface to a body wave mode. Note that when such an event occurs, the particle is reinjected at a random depth in the medium following the probability distribution $p_{s}(z)$. To obtain energy envelopes, the position and mode of the particle is monitored on a cylindrical grid at regular time intervals. The local energy density is estimated by averaging the number of particles per cell over a sufficiently large number of random walks. For accuracy, it is important that the cells be relatively small compared to the shortest mean free path in the medium.

5.2 Numerical results

Figure 3 illustrates the striking difference between the global and local partitioning of the seismic energy into surface and body waves. The following parameters have been employed in the simulation: $\alpha=1$km*-1*, $c=3$km/s and $\tau^{s\rightarrow s}=20$s, $\tau^{s\rightarrow b}=30$s, $\tau^{b\rightarrow b}=30$s, $\tau^{b\rightarrow s}(z)=10\exp(2z)$s. Note that in our model the group velocity of surface waves $v_{g}\approx 3.32$km/s is slightly faster than the speed of propagation of body waves. Two source depths are considered: a relatively shallow one ($z_{0}=1$km) and a deep one ($z_{0}=5$km), which radiate approximately $29\%$ and $0.01\%$ energy as surface waves, respectively.

On the left, we show the temporal evolution of the ratio between the total energy of surface waves $\bar{\bar{E}}_{s}$ (see the remarks before Eq. 56 for a reminder of the notations), and the horizontally-integrated energy density of body waves $\bar{E}_{b}$ at the surface $z=0$, averaged over a depth range $\Delta z=1$km. Hence, the ratio $\bar{\bar{E}}_{s}/\bar{E}_{b}$ has unit of inverse length. Independent of the source depth, we find that the partitioning of the energy density at the surface -into surfacic energy of surface waves and volumetric energy of body waves- converges toward the prediction of equipartition theory, at long lapse-time in the coda (see Eq. 25 and 60). This numerical result confirms that the analysis of equipartition given in the previous section in slab geometry extends to the half-space geometry. Furthermore, we find that the surface-to-body energy ratio overshoots the prediction of equipartition theory for the two sources at short lapse-time, by a factor which decreases with the source depth $z_{0}$.

The stabilization of the local energy density ratio of surface and body waves at the surface of the half-space is to be contrasted with the evolution of the global partitioning of the energy into surface and body wave modes. Figure 3 (right) shows that after a few mean free times, most of the energy is carried in the form of body waves in the medium. The Figure also suggests that the global transfer of energy from surface waves to body waves occurs at a rate proportional to $t^{-1/2}$ at long lapse-time.

Further insight into the equipartition process is offered in Figure 4, where we show the depth dependence of the horizontally-averaged body wave energy density at different lapse-time in the coda. All the parameters of the simulation are the same as in Figure 3, except for the much finer spatial resolution $\Delta z=\alpha^{-1}/5=0.2$ km, which allows us to track processes that occur in the skin layer where the coupling between surface and body waves occurs. We observe that after roughly 10 mean free times, the depth distribution of body wave energy becomes homogeneous over a depth range at least as large as $10\alpha^{-1}$, independent of the source depth. This simulation therefore confirms the theoretical analysis performed in slab geometry. The homogenization of the energy of body waves is a dynamic process: the energy density of surface waves increases exponentially near the surface, thereby generating a larger amount of body-wave converted energy; this process is compensated by the exponential increase of the conversion rate from body to surface waves, which eventually yields an equilibrium. Note that the total energy density does not homogenize with depth, due to the exponential decay of the surface wave eigenfunction with depth.

In Figure 5 and 7, we illustrate in greater details the multiple-scattering process by showing snapshots of the surfacic and volumetric energy densities $\varepsilon_{s}(t,\mathbf{r})$ and $E_{b}(t,\mathbf{r},z)|_{z=0}$ at regular time intervals $\Delta t=1\tau^{s}$ starting at a lapse time $t=0.8\tau^{s}$ for a shallow ($z_{0}=1$km) and a deep ($z_{0}=5$km) source, respectively. The scattering parameters are the same as in Figure 3 and the energy is averaged over a range of epicentral distance $\Delta r=5$km and depth $\Delta z=5$km. We use a double horizontal axis on Figures 5-7 to show simultaneously the epicentral distance in kms and in units of mean free path. Note that in the case of body waves, we take the value of the mean free path in the bulk of the medium $\tau^{b\rightarrow b}$. For comparison, we show in Figure 6, snapshots of energy density of surface waves and body waves when the coupling between the two is deactivated. In the case of surface waves, this amounts to computing the solution of a conventional 2-D multiple-scattering process with the mean free time $\tau^{s}$. In the case of body waves, we consider a conventional 3-D multiple-scattering process in a half-space with a constant mean free time $\tau^{b\rightarrow b}$, i.e., we remove the boundary layer where the coupling with surface waves occurs. To facilitate the comparison between Figure 5 and 6, we have adjusted the strength of the source term in the conventional multiple-scatttering simulations so that they match exactly the energy released at the source in the form of body and shear waves in the coupled case.

We first analyze the transport of surface waves in the case of a shallow source. As compared to the conventional 2-D case, mode coupling has at least two visible effects on the spatial distribution of the surface wave energy. First, it lowers the energy level in the coda. As an illustration, we observe that after 10 mean free times the coda intensity is reduced by a factor at least equal to 10 in Figure 5 compared to Figure 6. Second, mode coupling appears to enhance the visibility of ballistic surface waves. In Figure 6, the ballistic term is completely masked by the diffuse contribution at roughly 6 mean free path from the source, whereas a small ballistic peak is still visible at roughly 10 mean free paths from the source in Figure 5. It is worth noting that the ballistic contribution is exactly identical in Figures 5 (left) and 6 (left). Again, this is the strong decrease of the energy of scattered coda waves which explains the difference between the two Figures. Figure 8, which displays the spatial distribution of the mean order of scattering in the coda, reveals that the coda of coupled surface waves is depleted in high-order multiply-scattered waves compared to a conventional 2-D transport process. In other words, mode conversions entail a strong conversion of multiply-scattered surface waves into body waves which decreases the energy level of surface wave coda and, by comparison, enhances the ballistic contribution. Examination of Figures 5 (right) and 6 (right) reveals that the effects of mode coupling on body waves are opposite to the ones just described for surface waves. Thus, we observe that the energy level in the coda is slightly increased by the transfer energy from surface wave to body waves. An additional contribution comes from the increase of the scattering strength of body waves near the surface which attenuates the ballistic waves and transfers their energy into the coda. Examination of the decay of the ballistic peak of body waves with epicentral distance in Figures 5 and 6 confirms the increased attenuation entailed by the coupling with surface waves. Other more exotic phenomena are also visible in Figure 5 such as some precursory body waves arrival due to the coupling from surface waves to body waves. However this process is a very peculiar feature of our model, due to the higher wavespeed of surface waves compared to the one of body waves.

Further differences between our coupled model for surface and body waves and conventional transport theory is illustrated in Figure 7 where we show snapshots of the energy distribution of surface and body waves in the case of a deep source. Note that in that case, surface waves can only be generated by mode conversions so that ballistic arrivals are absent in Figure 7 (left). Interestingly, our numerical simulations indicate that surface waves are rapidly excited to a non-negligible level in the coda. Examination of Figure 8 (right) further indicates that multiple-scattering is at the origin of the generation of surface waves in the coda when the source is located at large depth. These observations agree with our theoretical analysis of equipartition, which implies that, independent of the source depth, the coda at the surface of a half-space always appears as a mixture of surface and body waves.

In Figure 9, we show envelopes of energy densities for surface and body waves in the case of a shallow source ($z_{0}=1$km) and a deep source ($z_{0}=5$km) at an epicentral distance of 50km. The scattering parameters are the same as in all previous Figures and the spatial resolution of the computation is 5km again. The impact of the depth of the source on the excitation of ballistic waves is obvious in Figure 9 and confirms the analysis of Figures 5 and 7. In particular, it is apparent that the direct body waves are less attenuated in the case of the deep source, as a consequence of the exponential decay of the scattering conversions from body to surface waves with depth. To facilitate the identification of different propagation regimes in the coda, we have superposed on the graphs some typical algebraic decays: $t^{-1}$ for scattering in 2-D (either single or multiple, see e.g. Paasschens, 1997), $t^{-3/2}$ for multiple scattering in 3-D, and $t^{-2}$ for single-scattering in 3-D. In Figure 9, we observe that for both body and surface waves the coda obeys a $t^{-3/2}$ decay law at long lapse-time, independent of the source depth, which is characteristic of a 3-D diffusion process. This supports the predictions of the diffusion model and confirms the dominance of body waves in the transport process at large lapse-time. At short lapse-time, we observe a distinct behavior between the two kinds of waves, particularly in the case of a shallow source. After the passage of the ballistic waves, body waves appear to decay slightly more slowly than the asymptotic $t^{-3/2}$ behavior. This may reflect the conversion of surface waves to body waves as discussed in the analysis of Figure 5. Two propagation regimes show up clearly on the surface wave energy envelope, with a transition between the two around a lapse-time of 10 mean free times. At short time, the decay of surface waves appears to be faster than the one of body waves, probably as a consequence of the transfer of surface wave energy into the volume as discussed in relation with Figure 5. Taken together, Figures 5-9 illustrate the much richer behavior of the coupled system of transport Eq. (50), compared to the conventional transport process without coupling between surface and body waves.

6 Conclusions

This work represents a first attempt at formulating a self-consistent theory of RT of seismic waves in a half-space geometry including the coupling between surface and body waves. The main approximation underlying our work is that, upon reflection at the surface, the phase of body waves is randomized so that upgoing and downgoing fluxes may be considered as statistically independent. Our approach distinguishes itself from the standard Eqs of RT for scalar waves found in the literature in one important way: it keeps track - to some extent - of the wave behavior in the vicinity of the surface. This has a number of consequences: (1) surface and body waves are coupled by conversion scattering (2) even in a statistically homogeneous medium it requires that the scattering properties of body waves depend on the depth in the medium. Furthermore, the reciprocity relation between the surface-to-body vs body-to-surface mean free times plays a prominent role in the establishment of an equipartition regime with a ratio that conforms to the predictions of standard mode counting arguments. Besides equipartition, a notable outcome of our RT equations is the anisotropy of the diffusivity of seismic waves, due to the difference in scattering properties and wave velocities of body and surface waves. We also show that our RT Eqs are operational, in the sense that they are readily amenable to numerical solutions by Monte-Carlo simulations. These simulations could be used in the future to study in more details the dynamics of equipartition, in particular, how the equipartition time varies as a function of the ratio between the penetration depth of surface waves and the scattering mean free path for body-to-surface wave coupling.

Before becoming a viable alternative to current approaches, our theory needs to be tested and improved. In the future, we plan to address the following issues: (1) Evaluate the impact of neglecting the interference between upgoing and downgoing body waves on the scattering cross-section and, if possible, go beyond this approximation. (2) Extend the theory to more realistic finite size scatterers and more general spatial distributions of scatterers. (3) Incorporate polarization effects for elastic waves at a free surface. (4) Absorption of energy is also a very important mechanism of attenuation, which has been entirely neglected in this work for simplicity. Because the sub-surface of the Earth is thought to be very strongly attenuating due to the widespread presence of fluids, we may expect dissipation to affect more severely surface waves than body waves. In turn, this may modify the partitioning of the energy in the coda as was previously shown by Margerin et al. (2001) in the case of coupled $S$ and $P$ waves. Special efforts should be devoted to this important topic before our formalism can be applied to real seismic data.

Acknowledgements.

The authors wish to thank the Associate Editor S. Ni and an anonymous referee for their suggestions to clarify the presentation of the results. The careful comments and constructive criticisms of H. Sato contributed to significant changes and improvements in the content of the manuscript. The authors acknowledge the European Research Council under the European Union Horizon 2020 research and innovation program (grant agreement no. 742335 - F-IMAGE).

Appendix A Variational formulation for mixed boundary conditions

Here, we recall briefly on a simple one-dimensional example how mixed boundary conditions of the type used in Eq. (2) can be incorporated into a variational formulation. The interested reader will find further details and more examples in the classic book by Gelfand & Fomin (1963), after which our treatment is modeled. For simplicity we consider a vibrating string of density $\rho$, tension $T$ and length $L$. For the moment, we do not specify the boundary conditions. The total kinetic energy stored in the string at time $t$ is given by:

[TABLE]

where $u$ denotes the displacement field. The instantaneous potential energy stored in the string may be expressed as

[TABLE]

According to Hamilton’s principle, among all possible displacement fields, the one that satisfies the actual equations of motion should make the following action integral:

[TABLE]

stationary. Mathematically, this principle of stationary action may be expressed as:

[TABLE]

where $\psi$ is an arbitrary function. This is sometimes written as $\delta I=0$, where $\delta I$ is known as the first variation of the action integral. Using integration by parts, Gelfand & Fomin (1963) establish that:

[TABLE]

The arbitrariness of the function $\psi$ in Eq. (LABEL:deltai) implies both the governing wave equation for the vibrating string:

[TABLE]

as well as the so-called natural boundary conditions:

[TABLE]

which correspond to a string with free ends. In order to obtain mixed boundary conditions, it suffices to add to the potential energy (73), a term of the form $\chi u(0,t)^{2}$ where $\chi$ is a constant. Eq. (LABEL:deltai) must be modified accordingly by adding the new contribution $-\epsilon\chi\int_{t_{0}}^{t_{1}}u(0,t)\psi(0,t)dt$ which, in turn, implies a natural boundary condition of the mixed type at $x=0$:

[TABLE]

The total potential energy may be rewritten in integral form as follows:

[TABLE]

which justifies the appearance of the delta function in Eq. (3) and Eq. (6) in a simplified context.

Appendix B Far-field expression of the Green’s function for scalar waves in a half-space with mixed B.C.

In this Appendix, we summarize the key steps to the derivation of Eq. (16 ) from Eq. (15). We split the computation into two parts and begin with the surface wave contribution:

[TABLE]

Introducing cylindrical coordinates $(k_{\parallel},\phi)$ and integrating over angle yields:

[TABLE]

where $J_{0}$ denotes the standard Bessel function of order 0. Using the same trick as in Aki & Richards (2002, Chapter 6), we extend the wavenumber integral over the whole $k_{\parallel}$ axis using the Hankel function of the first kind instead of the Bessel function:

[TABLE]

In the last step, we employ the residue theorem by closing the contour in the upper half of the complex plane with a semi-circle of radius $R\to+\infty$ and note the presence of pole at $k_{\parallel}=\sqrt{k^{2}+\alpha^{2}}+i\eta$, $(\eta\to 0^{+})$. Thanks to the exponential decay of the integrand, the integral on the semi-circle vanishes which yields:

[TABLE]

The result (84) is exact. The far-field approximation (16) follows by application of standard asymptotic expansions to the Hankel function.

The computation of the body wave contribution can also be split into two parts:

[TABLE]

where the unitarity of the reflection coefficient has been used and the $q$ integral has been extended from $-\infty$ to $+\infty$ thanks to the relation $r(q)^{*}=r(-q)$. The first term in the second equality of (85) may be recognized as the full-space solution to the Helmholtz Eq.:

[TABLE]

where $R_{0}=\sqrt{r^{2}+(z-z_{0})^{2}}$. The second term in the second equality of (85) represents the waves reflected at the surface:

[TABLE]

The computation of this integral may be attacked in exactly the same way as we did for the surface wave term $G_{s}$ to obtain:

[TABLE]

To approximate this last integral in the far-field of the source, we first remark that for $|q|>k$ the cylindrical waves are evanescent so that we may legitimately take $-k$ and $+k$ as integration limits. We next make use of the far-field expansion of the Hankel function to obtain the following oscillatory integral representation:

[TABLE]

Further noting that the derivative of the phase term:

[TABLE]

vanishes at :

[TABLE]

with $R^{\prime}_{0}=\sqrt{r^{2}+(z+z_{0})^{2}}$, we apply the stationary phase formula to obtain after some straightforward algebra:

[TABLE]

This term may be interpreted as the contribution of the image point of the source with a strength given by the reflection coefficient evaluated at an incidence angle corresponding to the specularly reflected ray connecting the source to the detection point (see Eq. 91). To complete the far-field approximation, we first note the following expansions: $R_{0}=R-z_{0}z/R+o(1/R)$, $R^{\prime}_{0}=R+z_{0}z/R+o(1/R)$ where $R=\sqrt{r^{2}+z^{2}}$. Neglecting all terms smaller than $1/R$ for the amplitude, all terms smaller than $z_{0}/R$ for the phase and further approximating $q_{0}$ as $kz/R$, formula (16) is recovered.