A scalar radiative transfer model including the coupling between surface and body waves
GGeophys. J. Int. (0000) , 000–000
A scalar radiative transfer model including the couplingbetween surface and body waves
Ludovic Margerin , Andres Bajaras , Michel Campillo Institut de Recherche en Astrophysique et Plan´etologie, Observatoire Midi-Pyr´en´ees,Universit´e Paul Sabatier, C.N.R.S., C.N.E.S., 14 Avenue Edouard Belin, 31400 Toulouse, France. Institut des Sciences de la Terre, Observatoire des Sciences de l’Univers de Grenoble,Universit Grenoble Alpes, C.N.R.S., I.R.D., CS 40700, 38058 GRENOBLE Cedex 9, France
15 July 2019
SUMMARY
To describe the energy transport in the seismic coda, we introduce a system of radiativetransfer equations for coupled surface and body waves in a scalar approximation. Ourmodel is based on the Helmholtz equation in a half-space geometry with mixed boundaryconditions. In this model, Green’s function can be represented as a sum of body wavesand surface waves, which mimics the situation on Earth. In a first step, we studythe single-scattering problem for point-like objects in the Born approximation. Usingthe assumption that the phase of body waves is randomized by surface reflection orby interaction with the scatterers, we show that it becomes possible to define, in theusual manner, the cross-sections for surface-to-body and body-to-surface scattering.Adopting the independent scattering approximation, we then define the scattering meanfree paths of body and surface waves including the coupling between the two types ofwaves. Using a phenomenological approach, we then derive a set of coupled transportequations satisfied by the specific energy density of surface and body waves in a mediumcontaining a homogeneous distribution of point scatterers. In our model, the scatteringmean 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 wavesis established at long lapse-time, with a ratio which is predicted by usual mode countingarguments. We derive a diffusion approximation from the set of transport equationsand show that the diffusivity is both anisotropic and depth dependent. The physicalorigin of the two properties is discussed. Finally, we present Monte-Carlo solutions ofthe transport equations which illustrate the convergence towards equipartition at longlapse-time as well as the importance of the coupling between surface and body wavesin the generation of coda waves.
Key words:
Coda waves, Wave scattering and diffraction, Theoretical seismology
In seismology, Radiative Transfer (RT) has been used for more than three decades to characterize the scattering and absorptionproperties of Earth’s crust (see for instance Fehler et al. 1992; Hoshiba 1993; Carcol´e & Sato 2010; Eulenfeld & Wegler 2017).Since its introduction by Wu (1985) for scalar waves in the stationary regime, RT has been considerably improved to bringit in closer agreement with real-world applications. In particular, the coupling between shear and compressional waves wasdeveloped 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 impressiveagreement between observations and elastic RT theory. For comprehensive introductions to RT, the reader is referred to thereview 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 simulationsof 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 a r X i v : . [ phy s i c s . g e o - ph ] J u l L. Margerin et al. heterogeneity stratification (Hoshiba 1997; Margerin et al. 1998; Yoshimoto 2000), coupling between shear and compressionalwaves (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 inseismology: 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 (Mancinelliet al. 2016a; Takeuchi et al. 2017), modeling of propagation anomalies in the crust (Sens-Sch¨onfelder 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 codawave 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 mappingrelies 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 ofspatio-temporal convolutions of the mean intensity in the coda. It was later pointed out by Margerin et al. (2016) that anaccurate 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 directlypredicted 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)- thevirtual sources and receivers are located at the surface of the medium so that the early coda is expected to contain a significantproportion of Rayleigh waves. At longer lapse-time, the surface waves couple with body waves and the coda eventually reachesan 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 ofthe composition of the coda wavefield at a given lapse-time in the coda is key to locate accurately the changes at depth inthe crust.Obermann et al. (2013, 2016) proposed to express the sensitivity of coda waves as a linear combination of the sensitivityof surface and body waves, whose kernels are computed from scalar RT theory in 2-D and 3-D, respectively. The relativecontribution of the 2-D and 3-D sensitivity kernels at a given lapse-time in the coda is determined by fitting the traveltimeshift predicted by the theory against full wavefield numerical simulations in scattering media, where the background seismicvelocity is perturbed in a fine layer at a given depth. This method has been validated in the case of 1-D perturbationsthrough numerical tests and has the advantage of modeling exactly the complex coupling between surface and body wavesin heterogeneous media. Furthermore, it can easily incorporate realistic topographies, which is important for the monitoringof volcanoes. The main drawbacks of the approach of Obermann et al. (2016) are the numerical cost and the fact that itrequires a good knowledge of the scattering properties of the medium, which have to be determined by other methods suchas 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 thecoupling between surface and body waves in a self-consistent way. In the case of a slab bounded by two free surfaces, Tr´egour`es& Van Tiggelen (2002) derived from first principles a quasi 2-D RT equation where the wavefield is expanded onto a basisof Rayleigh, Lamb and Love eigenmodes. Thanks to the normal mode decomposition, this model incorporates the boundaryconditions at the level of the wave equation. The energy exchange between surface and body waves is treated by normal modecoupling in the Born approximation. A notable advantage of this formulation is the capacity to predict directly the energydecay in the coda and its parttioning onto different components. The two main limitations for seismological applications arethe slab geometry, which may not always be realistic and the fact that the disorder should be weak, i.e., the mean free timeshould 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 wavesand body waves in the seismic coda. The formalism used by the author is interesting and bears some similarities with theone 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 differssignificantly from our findings. A very promising investigation of the energy exchange between surface and body waves onthe basis of the elastodynamic equations in a half-space geometry was performed by Maeda et al. (2008). Using the Bornapproximation, these authors calculated the scattering coefficients between all possible modes of propagation in a mediumcontaining random inhomogeneities. The main limitation of their theory comes from the fact that the conversion from bodyto surface waves is quantified by a non-dimensional coefficient, which makes it difficult to extend their results beyond thesingle-scattering regime. The authors argue that the absence of a characteristic scale-length for body-to-surface attenuation isa 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 tothe one of Maeda et al. (2008). For simplicity, we limit our investigations to a scalar model based on the Helmholtz equationwith an impedance (or mixed) boundary condition in a half-space geometry. To make the presentation self-contained, wereview the most important features of this particular wave equation. Specifically, we recall that the modes of propagation adiative transfer of body and surface waves 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 surfacereflexion randomizes the phase of the reflected wave, we are able to derive simple expressions for the scattering mean freepath of both body and surface waves including the coupling between the two. We elaborate on this result to establish aset of two coupled RT equations satisfied by the specific energy density of surface and body waves using a phenomenologicalapproach. Some consequences of our simple theory are explored, in particular the establishment of a diffusion and equipartitionregime. Monte-Carlo simulations show the potential of the approach to model the transport of energy in the seismic codafrom single-scattering to diffusion. In this section, we present the basic ingredients of our scalar model based on the Helmholtz equation. We describe how anappropriate 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 ofstates, 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).
We consider a 3-D version of the membrane vibration equation in a half-space geometry:( ρ∂ tt − T ∆) u ( t, R ) = 0 (1)where t is time and R is the position vector. It may be further decomposed as R = r + z ˆz ( z ≥
0) with r = x ˆx + y ˆy and( ˆx , ˆy , ˆz ) denotes a Cartesian system. In Eq. (1) u , ρ and T may be thought of as the displacement, the density and the elasticconstant of the medium, respectively. The wave Eq. (1) is supplemented with the boundary conditions:( ∂ z + α ) u ( t, R ) | z =0 = 0+radiation condition at ∞ (2)The case of interest to us corresponds to α >
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: L = 12 (cid:2) ρ ( ∂ t u ( t, R )) − T ( ∇ u ( t, R )) + αT u ( t, R ) δ ( z ) (cid:3) (3)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 boundaryconditions at one end and free boundary conditions at the other end.In the case of a harmonic time dependence u ∝ e − iωt ( ω > u ( R ) + ω c u ( R ) = 0 (4)with c = (cid:112) T /ρ the speed of propagation of the waves in the bulk of the medium. Eq. (4) is complemented with the mixedboundary condition ∂ z u + αu = 0 at z = 0 and an outgoing wave condition at infinity. From (3), we can deduce the energyflux density vector J and the energy density w using the concept of stress-energy tensor (Morse & Ingard 1986). For harmonicmotions, their average value over a period can be expressed as: J = − T { iωu ∗ ∇ u } (5) w = 14 (cid:8) ρω | u | + T |∇ u | − αT | u | δ ( z ) (cid:9) (6)In the following section, we recall the consequences of mixed boundary conditions on the Helmoltz Eq., in particular the factthat it gives rise to a surface wave mode. Due to the translational invariance of the medium, we look for eigen-solutions of Eq. (4) in the form u = ψ ( z ) e i k (cid:107) · r with k (cid:107) = ( k x , k y , z variable only. For α >
0, part of the spectrum is
L. Margerin et al. Normalized Frequency0123456 N o r m a li z e d V e l o c i t y PhaseBulkGroup
Figure 1.
Dispersion of surface waves in a half-space with mixed boundary conditions. On the horizontal axis, the normalized frequencyis defined as ωα/c . On the vertical axis, the velocity is normalized by the speed of body waves c . discrete with eigenfunction : u s ( r , z ) = √ αe − αz e i k (cid:107) · r π (7)with k (cid:107) · k (cid:107) − α = ω c . The rest of the spectrum forms a continuum of body waves with normalized eigenfunctions: u b ( r , z ) = 1(2 π ) / ( e − iqz + r ( q ) e iqz ) e i k (cid:107) · r , q ≥ q + k (cid:107) · k (cid:107) = ω c and: r ( q ) = q + iαq − iα (9)Note the relations: (1) q = ( ω cos j ) /c = k cos j with j the incidence angle of the body wave and (2) | r ( q ) | = 1, i.e., thereis total reflection at the surface. For later reference, we introduce a specific notation for the vertical eigenfunction of bodywaves: ψ b ( ˆn , z ) = ( e − ikz ˆn · ˆz + r ( k ˆn · ˆz ) e ikz ˆn · ˆz ) (10)with ˆn · ˆz = cos j . Note that throughout the paper, we use a hat to denote a unit vector. The surface waves (7) and bodywaves (8) are normalized and orthogonal in the sense of the scalar product (cid:104) u | v (cid:105) = (cid:82) R u ( R ) ∗ v ( R ) d R , where u and v arearbitrary square integrable functions. The surface wave phase velocity c φ is given by: c φ = c (cid:114) c α ω (11)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: v g = dωdk = c c φ = c (cid:114) c α ω (12)and (2) using the principle of energy conservation: v Es = (cid:104) J (cid:105)(cid:104) w (cid:105) = v g ˆk (cid:107) (13)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 adiative transfer of body and surface waves Figure 2.
Geometry of the surface employed to compute the energy radiation of a point source located at (0 , , z ). C is a cylindricalsurface of radius R c and height h . H is a hemispherical surface of radius R . phase and group velocity tend to the common value c . Using definition (13) it is possible to define the energy velocity of abody wave eigenmode (see Eq. 8): v Eb = lim h →∞ (cid:104) J (cid:105) h (cid:104) w (cid:105) h = c sin j ˆk (cid:107) , (14)where (cid:104)(cid:105) h denotes an integration from the surface to depth h . The depth averaging smoothes out the oscillations of J and w caused by the interference between the incident and reflected amplitudes. The passage to the limit is necessary becausethe integrals over depth diverge. Eq. (14) can be interpreted as follows. In the full-space case, the current vector of a singleunit-amplitude plane wave with wavevector k = k (cos j ˆz + sin j ˆk (cid:107) ) is given by J = ρωc k / w = ρω /
2. If we define k r as the mirror image of k across the plane z = 0 and consider the sum of the current vector of twoplane waves with wavevectors k and k r we obtain J + J r = ρω c sin j ˆk (cid:107) . 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 ofthe 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. withmixed BC in the form: G ( r , z, z ) = 1(2 π ) (cid:90) + ∞ dq (cid:90) R e i k (cid:107) · r ( e − iqz + r ( q ) e iqz )( e − iqz + r ( q ) e iqz ) ∗ k − k (cid:107) − q + i(cid:15) d k (cid:107) + 2 αe − α ( z + z ) (2 π ) (cid:90) R e i k (cid:107) · r k + α − k (cid:107) + i(cid:15) d k (cid:107) , (15)where z denotes the source depth, (cid:15) is a small positive number which guarantees the convergence of the integrals and thestar ∗ 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 stationaryphase approximation for the body wave term: G ( r , z, z ) = − e ikR πR ψ b ( ˆR , z ) − αe − α ( z + z )+ ik s r + iπ/ √ πk s r (16)with k s = ω/c φ , ˆR = R /R and R = r + z ˆz . The expansion (16) is performed with respect to the midpoint of the source pointand 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. We now compute the energy radiated by a point source located at (0 , , z ). To do so, we introduce a cylindrical surface C ofradius R c extending from the free surface to a depth h greater than z and large compared to 1 /α . We close this surface witha hemispherical cap H of radius R centered at the surface point (0 , , L. Margerin et al.
2. The energy flux vector (5) of the radiated field contains terms that are purely surface, purely bulk and cross-terms. Thecontribution 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 subtendsa solid angle which goes to 0 as R c increases. The cross-terms are negligible across the whole surface because the coupledsurface/body wave term decays algebraically faster than the surface wave term on the cylindrical surface and exponentiallyfaster than the body wave term on the hemispherical surface. Hence, we may split the flux of radiated waves into a contributionof surface and body waves, respectively.The energy transported per unit time by body waves through the hemispherical cap H is given by: E b ( z ) = ρω c (cid:90) H | G b ( r , z, z ) | R d ˆ R = ρω c π (cid:90) π | ψ b ( ˆR , z ) | d ˆ R, (17)with G b the body wave part of Green’s function. In the second line of Eq. (17), the integral is over the space directionssubtended by the hemispherical cap. (N.B.: strictly speaking, the total solid angle is not equal to 2 π because one shouldremove 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: (cid:104)E b (cid:105) = ρω c π (18)Here the brackets may have at least 2 different meanings. The most obvious is an average over depth, as was done in thecalculation of the group velocity. But we may also assume that the surface “scrambles” the phase of the reflected wave φ r sothat it becomes a random variable. In this scenario, the brackets would mean an average over all realizations of the randomreflection process. Upon averaging over phase or depth, the interference pattern between the incident and reflected wave issmoothed out, so that the two approaches yield the same result. Note that the randomization of the phase does not affectenergy conservation because the incident flux is still totally reflected. In particular, the discussion following the interpretationof Eq. (14) would still be valid. In practice, the assumption that the phase of the reflected wave is randomized by the surfacemay not be as unrealistic as it seems. Observations of reflected SH waves by Kinoshita (1993) at borehole stations in Japanindeed suggest that the reflected field is a distorted version of the incident one. This concurs with the general view that thesubsurface of the Earth is highly heterogeneous at scales that can be much smaller than the wavelength and brings support tothe 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: E s ( z ) = ρω v g (cid:90) C | G s ( r , z, z ) | rdφdz = ρω v g α e − αz πk s (cid:90) π dφ (cid:90) h e − αz dz = ρωc α e − αz (19)with G s the surface wave part of Green’s function. Because h is large compared to 1 /α the integral over depth may beperformed from 0 to + ∞ (the error incurred is exponentially small). We find the depth dependent surface-to-body energyratio: R ( z ) = E s ( z ) (cid:104)E b (cid:105) = 2 πcαω e − αz (20)Formula (20) can also be understood in the light of the local density of states n s,b defined as (Sheng 2006): n s,b ( z ) = − Im G s,b ( r , z , z ) π (cid:12)(cid:12)(cid:12)(cid:12) r = ,z = z × dk s,b ( ω ) dω , (21)where k s,b ( ω ) stands for the wavenumber of surface or body waves. Using the spectral representation (15), one obtains the(exact) formulas: n b ( z ) = ω π c (cid:90) π | ψ b ( ˆR , z ) | d ˆ R (22) (cid:104) n b (cid:105) = ω πc (23) n s ( z ) = αωc e − αz (24)which show that the partitioning of the energy radiated into surface and body waves by the source R ( z ) is given by the ratioof their local density of states n s ( z ) / (cid:104) n b (cid:105) . Finally, we may compute the partitioning ratio R between the modal density of adiative transfer of body and surface waves surface and body waves by integrating Eq.(24) over z and taking the ratio with (23). This yields the simple result: R = (cid:90) ∞ n s ( z ) (cid:104) n b (cid:105) dz = πcω (25)where it is to be noted that the modal density ratio R is independent of the scale length α appearing in the mixed boundarycondition of the Helmholtz equation. This result could have been deduced directly from the dispersion relations of body andsurface 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 carriedover one hemisphere only, each eigenmode ψ b is composed of an incident and a reflected wave, which -on average- doublesits contribution compared to a single plane wave state. In the next section, we use our knowledge of the Green’s function toderive the scattering properties of surface and body waves including the coupling between the two modes of propagation. In this section, we calculate the energy radiated by a single scatterer in a half-space geometry for incident surface or bodywaves. For simplicity, we restrict our investigations to point scatterers and employ Born’s approximation. The resultingexpressions are simplified following the scrambled phase approximation and interpreted in terms of scattering cross-sections.
We now consider the following perturbed Helmholtz Eq.:∆ u ( r , z ) + k (1 + (cid:15)a δ ( r ) δ ( z − z s )) u ( r , z ) = 0 (26)Here a represents the typical linear dimension of the scatterer located at (0 , , z s ) and it is understood that ka (cid:28) (cid:15) is thelocal perturbation of inverse squared velocity. Following the standard procedure (Snieder 1986), we look for solutions of Eq.(26) of the form: u = u + u s , where u ( r , z ) = e − αz + i √ α + ω /c x is a surface wave eigenmode of the Helmholtz equation (theincident field) and u s is the scattered field. Using the Born approximation, one obtains: u ( r , z ) = u ( r , z ) − k (cid:15)a G ( r , z, z s ) u ( , z s ) (27)Introducing the coupling strength S = k a (cid:15)e − αz s , one may express the energy radiated by the body waves through thehemispherical cap H (per unit time) as: U s → b = ρω S c π ) R (cid:90) π | ψ b ( ˆR , z s ) | R d ˆ R. (28)Note that the coupling term S depends on both intrinsic properties of the scatterer -size and strength of perturbations- aswell 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 | J s | = ρω v g α . (29)The ratio of (28) and (29) gives the surface-to-body scattering cross-section: σ s → b ( z s ) = αck a (cid:15) e − αz s π v g (cid:90) π | ψ b ( ˆR , z s ) | d ˆ R (30)This cross-section has unit of length. In the case where the surface scrambles the phase of the reflected wave, we may computethe mean conversion scattering cross-section by taking the average over the random phase φ r . (cid:104) σ s → b ( z s ) (cid:105) = cαk a (cid:15) e − αz s πv g (31)Let us remark again that the averaging smoothes out the interference pattern of the body wave eigenfunction ψ b but doesnot affect the conservation of energy. Furthermore, the averaging procedure makes the scattering pattern isotropic since (cid:104)| ψ b ( ˆR , z s ) | (cid:105) = 2 with an equal contribution of upgoing and downgoing waves. The computation of the surface-to-surfacescattering cross-section proceeds in a similar way. The surface-wave energy radiated through the lateral surface is given by: U s → s = ρωc α S π (cid:90) π (cid:90) + ∞ e − α ( z + z s ) dφdz = ρωc αS e − αz s σ s → s ( z s ) = cα k a (cid:15) e − αz s v g , (33) L. Margerin et al. again with unit of length. If we have a collection of point scatterers with volume density n , we may define a surface wavescattering mean free path using the independent scattering approximation as (Lagendijk & Van Tiggelen 1996; Tr´egour`es &Van Tiggelen 2002; Maeda et al. 2008): 1 l s = (cid:90) ∞ n (cid:16) σ s → s ( z ) + σ s → b ( z ) (cid:17) dz (34)The approximate formula (34) neglects all recurrent interactions between the scatterers, which is valid for sufficiently lowconcentrations of inclusions. In the case of an incident body wave mode (8), the computation of the energy radiated in the form of body or surface wavescan be performed as in the previous section. The definition of the scattering cross-section, however, is problematic if we stickto 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 thescatterer is located at a large depth in the half-space. In this case, it appears reasonable to think that the scattering ofa body wave mode should be equivalent to the scattering of a plane wave in the full-space case, at least in a sense to bespecified below. To verify the correctness of this assertion, we start by computing the body wave cross-section in absence ofa boundary. Assuming a unit amplitude incident plane wave and using again the Born approximation, one may express thescattered energy as: U full = ρω c ( k (cid:15)a ) π , (35)Normalizing the result by the energy flux density of the incident wave: | J b | = ρω c σ full = ( k (cid:15)a ) π (37)In the half-space geometry, the incident wave has the form (8). The total energy scattered in the form of body waves isstill given by Eq. (28) provided one redefines the coupling constant as S = k a (cid:15) | ψ b ( ˆR i , z s ) | , where ˆR i refers to the incidencedirection of the body wave. Eq. (28) differs from formula (35) as a consequence of the interference between the incident andreflected waves. In the case of a scattering medium, we may expect these interferences to be blurred due to the randomizationof the phase by the scattering events. In this scenario, the phase of the incident and reflected waves may be expected to beuncorrelated. Upon averaging the result (28) over the random phase of the reflected wave, one obtains: U b → b = ρω c ( k (cid:15)a ) π , (38)which is exactly the double of the full-space result (35). Keeping in mind that the energy flux of the incident plane waveinteracts twice with the scatterer (direct interaction + interaction after reflexion) and using the assumption that the energyfluxes of incident and reflected waves do not interfere and may therefore be added, we obtain the result: (cid:104) σ b → b (cid:105) = σ full (39)As announced, the average scattering cross-section of body waves in the half space is the same as in the full space for ascatterer located far away from the boundary. Conceptually, the “scrambled phase” approximation allows us to extend theresult (39) to a scatterer located at an arbitrary depth in the medium thanks again to the assumption that the incidentand reflected waves are statistically independent. In this scenario, on average, the surface does not modify anything to thescattering 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 energyradiated in the form of surface waves is given by: U b → s = ρωc α ( k (cid:15)a ) e − αz s | ψ b ( ˆR i , z s ) | (cid:104) σ b → s (cid:105) ( z ) = α ( k (cid:15)a ) e − αz s k (41)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 σ b → b and σ b → s have unit of surface. Using the independent scattering approximationagain, we conclude that the inverse scattering mean free path of body waves defined as:1 l b ( z ) = n ( σ b → s ( z ) + σ b → b ) (42) adiative transfer of body and surface waves depends on the depth in the medium, as a consequence of the coupling with surface waves. A simple but fundamentalreciprocity relation may be established between the surface-to-body and body-to-surface scattering mean free time. The latermay be expressed as: τ b → s ( z ) = 2 ke αz ncα ( k a (cid:15) ) (43)while the former may be obtained after performing the integral over depth in Eq. (34): τ s → b = 4 πnc ( k a (cid:15) ) (44)The ratio between the two quantities is given by: τ s → b τ b → s ( z ) = 2 παe − αz k = R ( z ) (45)where Eq. (20) has been used. Eq. (45) establishes a link between the surface-to-body versus body-to-surface conversion ratesand the local density of states. In the case of elastic waves, a similar relation applies (Weaver 1990; Ryzhik et al. 1996): theratio 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 surfaceand body waves. In this section, we employ standard energy balance arguments to derive a set of coupled equations of RT for surface andbody waves in a half-space containing a uniform distribution of point-scatterers. A notable feature of our formulation is theappearance of the penetration depth of surface waves as a parameter in the Equations.
Before we establish the transport equation, a few remarks are in order. It is clear that the transport of surface wave energy isnaturally described by a specific surface energy density (cid:15) s ( t, r , ˆn ) where ˆn is a unit vector in the horizontal plane. The surfaceenergy density of surface waves may in turn be defined as: ε s ( t, r ) = (cid:90) π (cid:15) s ( t, r , ˆn ) d ˆ n (46)where the integral is carried over all propagation directions in the horizontal plane. Note that the (cid:15) s symbol should not beconfused with the strength of fluctuations in Eq. (26). In contrast with surface waves, the transport of body waves is describedby a specific volumetric energy density e b ( t, r , ˆk , z ) where ˆk is a vector on the unit sphere in 3-D. The volumetric energydensity of body waves is again obtained by integration of the specific energy density over all propagation directions in 3-D: E b ( t, r , z ) = (cid:90) π e b ( t, r , z, ˆk ) d ˆ k (47)In the usual formulation of transport equations, energy densities have the same unit (either surfacic or volumetric). In orderto treat on the same footing surface and body waves, we introduce the following volumetric energy density of surface wavesas: e s ( t, r , z, ˆn ) = 2 α(cid:15) s ( t, r , ˆn ) e − αz (48)It is clear that upon integration of e s over depth, one recovers the surface density (cid:15) s . The exponential decay of the surfacewave energy density is directly inherited from the modal shape and implies that the coupling between surface and body wavesmostly occurs within a skin layer of typical thickness 1 / α . For future reference, we introduce the following notation: E s ( t, r , z ) = (cid:90) π e s ( t, r , z, ˆn ) d ˆ n (49)to represent the volumetric energy density of surface waves, consistent with Eq. (47). Note that the total energy density ata given point will be defined as the sum ( E b ( t, r , z ) + E s ( t, r , z )), thereby implying the incoherence between the two types ofwaves.With this definition of the surface wave energy density, the phenomenological derivation of the radiative transport equationfollows exactly the same procedure as in the multi-modal case (Turner & Weaver 1994). A beam of energy followed along itspath around direction ˆn i is affected by (1) conversion of energy into other propagating modes and/or deflection into otherpropagation directions ˆn o (cid:54) = ˆn i ; (2) a gain of energy thanks to the reciprocal process: a wave with mode i propagating indirection ˆn i can be converted into a wave with mode o propagating in direction ˆn o by scattering. In the general case of finite sizescatterers we may anticipate scattering to be anisotropic. To describe such an angular dependence of the scattering process, we L. Margerin et al. may introduce normalized phase functions p i → o ( ˆn o , ˆn i ). The phase function may be understood as the probability that a waveof mode i propagating in direction ˆn i be converted into a wave of mode o propagating in direction ˆn o . To be interpretableprobabilistically, its integral over all outgoing directions ( ˆn o ) should equal 1. Note that in addition to the incoming andoutgoing propagation directions, the phase function may also depend on the depth in the medium as a consequence of thepresence of the surface. In the case of point scatterers, this complexity disappears thanks to the scrambled phase assumptionand will therefore not be considered in our formalism.A detailed local balance of energy yields the following system of coupled transport equations:( ∂ t + v g ˆn · ∇ ) e s ( t, r , z, ˆn ) = − e s ( t, r , z, ˆn ) τ s + 1 τ s → s (cid:90) π p s → s ( ˆn , ˆn (cid:48) ) e s ( t, r , z, ˆn (cid:48) ) d ˆ n (cid:48) + 1 τ b → s ( z ) (cid:90) π p b → s ( ˆn , ˆk (cid:48) ; z ) e b ( t, r , z, ˆk (cid:48) ) d ˆ k (cid:48) + s s ( t, r , z, ˆn ) (cid:16) ∂ t + c ˆk · ∇ (cid:17) e b ( t, r , z, ˆk ) = − e b ( t, r , z, ˆk ) τ b ( z ) + 1 τ b → b (cid:90) π p b → b ( ˆk , ˆk (cid:48) ) e b ( t, r , z, ˆk (cid:48) ) d ˆ k (cid:48) + 1 τ s → b (cid:90) π p s → b ( ˆk , ˆn (cid:48) ) e s ( t, r , z, ˆn (cid:48) ) d ˆ n (cid:48) + s b ( t, r , z, ˆn ) (50)where the terms s s,b represent sources of surface and body waves. To take into account the reflection of body waves at thesurface, the sytem (50) is supplemented with the following boundary condition for the energy density e b : e b ( t, r , , ˆk ) = e b ( t, r , , ˆk r ) (51)where ˆk r is the mirror image of the incident direction ˆk ( ˆk · ˆ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 , the terms s s and s b are given respectively by: s s ( t, r , z, ˆn ) = 2 αR ( z ) e − αz δ ( r )2 π (1 + R ( z )) s b ( t, r , z, ˆk ) = δ ( z − z ) δ ( r )4 π (1 + R ( z )) (52)where R ( z ) is the energy partitioning ratio defined in Eq. (20). Note that s s follows the same exponential decay as e s withdepth z (see Eq. 48), which takes into account the vertical dependence of the surface wave eigenfunction. In Eq. (52), thecomplex dependence of the body wave radiation with depth has been simplified by averaging the exact result (17) over therandom phase of the reflected wave. As a consequence, the energy radiation at the source covers uniformly the whole sphereof 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 surfacewas absent in Eq. (50). Again, this approximation is admissible if the scattered upgoing and downgoing energy fluxes arestatistically independent, a condition which is guaranteed by the randomization of the phase of the waves upon reflection atthe surface in our model. As discussed in the previous section, in the case of point scatterers the phase averaging proceduremakes all scattering processes isotropic which allows us to simplify the system of Eqs (50) by evaluating the scattering integralson the right-hand side of Eq. (50):( ∂ t + v g ˆn · ∇ ) e s ( t, r , z, ˆn ) = − e s ( t, r , z, ˆn ) τ s + E s ( t, r , z )2 πτ s → s + E b ( t, r , z )2 πτ b → s ( z ) + s s ( t, r , z, ˆn ) (cid:16) ∂ t + c ˆk · ∇ (cid:17) e b ( t, r , z, ˆn ) = − e b ( t, r , z, ˆk ) τ b ( z ) + E b ( t, r , z )4 πτ b → b + E s ( t, r , z )4 πτ s → b + s b ( t, r , z, ˆk ) (53)where we have used the normalization condition of the phase functions. The decreasing efficacy of scattering conversions withdepth is guaranteed by the exponential decay of τ b → s ( z ) − and E s ( t, 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 τ b ( z ) and τ b → s ( z ) is caused by the coupling with surface waves and not by a stratification of heterogeneity. A self-consistent formulation of coupled transport equations should verify two elementary principles: energy conservation andequipartition. 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 ˆn (in 2-D) or ˆk (in 3-D)yields: ∂ t E s ( t, r , z ) + ∇ · J s ( t, r , z ) = − E s ( t, r , z ) τ s → b + E b ( t, r , z ) τ b → s ( z ) ∂ t E b ( t, r , z ) + ∇ · J b ( t, r , z ) = − E b ( t, r , z ) τ b → s ( z ) + E s ( t, r , z ) τ s → b ,, (54) adiative transfer of body and surface waves where we have used the definition of the mean free times on the RHS of (54). Note that we have dropped the source termssince they are not essential to our argumentation. In Eq. (54), we have introduced the energy flux density vector of surfaceand body waves: J s = (cid:90) π e s ( t, r , z, ˆn ) v g ˆn d ˆ n J b = (cid:90) π e b ( t, r , z, ˆk ) c ˆk d ˆ k (55)Note that J s is contained in a horizontal plane. Upon integration over the whole half-space, the terms which contain thecurrent density vectors can be converted into surface integrals that vanish. Denoting by a double bar an integration over r and z , and summing the two Eqs of the system (54) leaves us with: ∂ t (cid:16) ¯¯ E s ( t ) + ¯¯ E b ( t ) (cid:17) = 0 (56)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 andfrom 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 αh (cid:29)
1, the mediummay be still be considered as a half-space in the treatment of surface waves. This leads us to: ∂ t ¯¯ E s ( t ) = − ¯¯ E s ( t ) τ s → b + (cid:90) h ¯ E b ( t, z ) τ b → s ( z ) dz∂ t (cid:90) h ¯ E b ( t, z ) dz = − (cid:90) h ¯ E b ( t, z ) τ b → s ( z ) dz + ¯¯ E s ( t ) τ s → b , (57)where the single bar denotes an integration over r only. Note that ¯¯ E s is the total energy of surface waves (up to an exponentiallysmall correction term), in contrast with ¯ E b which represents the body wave energy per unit depth. Our goal is not to solvethe system of Eq. (57) in its full generality but rather to exhibit an equipartition solution. In this regime, we expect thedistribution of body wave energy to become independent of depth. Hence we look for solutions of the system (57) of the form: (cid:32) ¯¯ E s ( t )¯ E b ( t, z ) (cid:33) = (cid:32) ¯¯ E s ¯ E b (cid:33) e − λt (58)Reporting the ansatz (58) into the integro-differential system (57), one arrives at a linear and homogeneous system of algebraicequations. A non-zero solution is obtained only if the determinant of the sytem vanishes which implies: λ (cid:18) λ − τ s → b − h (cid:90) h dzτ b → s ( z ) (cid:19) = 0 (59)One of the solutions is given by λ = 0 which corresponds to the asymptotic equipartition state such as:¯¯ E s ¯ E b = (cid:90) h τ s → b τ b → s ( z ) dz ≈ (cid:90) + ∞ R ( z ) dz = R (60)where Eq. (25) and (45) have been used, and the depth integral is again extended to + ∞ thanks to the assumption αh (cid:29) R at long lapse-time. Note that the energy density is not perfectly homogenized spatially, even at equipartition, asa consequence of the decay of the surface wave eigenfunction with depth. This may be related to the depth-dependence ofthe density of states near the boundary (see e.g. Hennino et al. 2001). The result (60) could have been predicted using theusual concept of equipartition which states that when filtered around a narrow frequency band, all the propagating modes ofa diffuse field should be excited to equal energy (Weaver 1982). In the case where the medium is unbounded at depth, therewill be a flux of energy across the lower boundary z = h which vanishes as the lapse-time increases. As demonstrated throughnumerical simulations later in this paper, equipartition also sets in in this configuration, though probably more slowly thanin the slab geometry. 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 equationbecause it exhibits a decay with depth which is inherited from the modal shape and therefore independent of the scatteringproperties. To circumvent the difficulty, we derive a closed diffusion equation for E b from which we subsequently deduce theenergy density of surface waves. To simplify the calculations, we assume that the scattering is isotropic. Proceeding in the L. Margerin et al. usual fashion, we expand the specific energy density into its first two angular moments (Akkermans & Montambaux 2007): e s ( t, r , z, ˆn ) = E s ( t, r , z )2 π + J s ( t, r , z ) · ˆn πe b ( t, r , z, ˆk ) = E b ( t, r , z )4 π + 3 J b ( t, r , z ) · ˆk π . (61)Multiplying, respectively, the first and second line of Eq. (53) by the unit vectors ˆn and ˆk , integrating over all possibledirections and employing the moment expansion (61), we obtain the following set of Equations: ∂ t J s ( t, r , z ) v g + v g ∇ (cid:107) E s ( t, r , z ) = − J s ( t, r , z ) v g τ s ∂ t J b ( t, r , z ) c + c ∇ E b ( t, r , z ) = − J b ( t, r , z ) cτ b , (62)where ∇ (cid:107) denotes the gradient operator in the horizontal plane. Note that the expansions (61) are used only to evaluate theintegral of the gradient on the left-hand side of the RT Equation. All other terms follow directly either from the definitionof the energy flux density vector or the assumption of isotropic scattering. Our final approximation consists in neglecting thederivative of the current vector with respect to time which yields the equivalent of Ohm’s law for the mutiply-scattered waves: J s ( t, r , z ) = − D s ∇ (cid:107) E s ( t, r , z ) J b ( t, r , z ) = − D b ( z ) ∇ E b ( t, r , z ) (63)where the following notations have been introduced: D s = v g τ s D b ( z ) = c τ b ( z )3 (64)Note that the total reflection condition (51) imposes that there is no net flux of body waves across the surface z = 0. To makeprogress, we now invoke the equipartition principle to fix the ratio between the energy densities of surface and body waves: E s ( t, r , z ) E b ( t, r , z ) = R ( z ) (65)in agreement with Eq. (60). This allows us to express the total energy flux as: J = J b + J s = − ( R ( z ) D s + D b ( z )) ∇ (cid:107) E b ( t, r , z ) + D b ( z ) ˆz ∂ z E b ( t, r , z ) (66)which may be interpreted as a generalization of Ohm’s law for diffuse waves. Eq (66) demonstrates that the coupling be-tween surface and body waves renders the energy transport both depth-dependent and anisotropic. It is worth emphasizingthat 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 bodyand 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 thesum of the set of Eqs (54): ∂ t E ( t, r , z ) + ∇ · J ( t, r , z ) = δ ( t ) δ ( z − z ) δ ( r ) , (67)where the right-hand side now contains the source term with z the source depth. Making use of Eqs (65) and (63), we obtainthe following diffusion-like equation verified by the body wave energy density: ∂ t E b ( t, r , z ) − ∇ (cid:107) · (cid:18) D b ( z ) + R ( z ) D s R ( z ) ∇ (cid:107) E b ( t, r , z ) (cid:19) −
11 + R ( z ) ∂ z ( D b ( z ) ∂ z E b ( t, r , z )) = δ ( t ) δ ( r ) δ ( z − z )1 + R ( z ) (68)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 )) − factor in front of the derivative operators. Actually, this difference is purely formal as may be shown by the change of variable z → z (cid:48) where z (cid:48) is defined as: z (cid:48) = (cid:90) z (1 + R ( x )) dx, = z + πk (cid:0) e − αz − (cid:1) (69)where Eq. (20) has been used. In the new variables, Eq. (68) may be rewritten as: ∂ t E (cid:48) b ( t, r , z (cid:48) ) − ∇ (cid:107) · (cid:0) D (cid:107) ( z (cid:48) ) ∇ (cid:107) E (cid:48) b ( t, r , z (cid:48) ) (cid:1) − ∂ z (cid:48) (cid:0) D ⊥ ( z (cid:48) ) ∂ z (cid:48) E (cid:48) b ( t, r , z (cid:48) ) (cid:1) = δ ( t ) δ ( r ) δ ( z (cid:48) − z (cid:48) ) (70)In Eq. (70), we have introduced the notations E (cid:48) b ( t, r , z (cid:48) ) = E b ( t, r , z ), z (cid:48) = (cid:82) z (1+ R ( x )) dx , as well as the following definitionsof the horizontal and vertical diffusivities: D (cid:107) ( z (cid:48) ) = D b ( z ) + R ( z ) D s R ( z ) D ⊥ ( z (cid:48) ) = D b ( z )(1 + R ( z )) (71) adiative transfer of body and surface waves 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 )) − 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 (cid:48) → ∞ . According to Eq. (69), this implies z (cid:48) ≈ z .Since the partitioning ratio R ( z ) goes to zero at large depth in the medium, Eq. (71) indicates that the vertical and horizontaldiffusivities become equal and the diffusion tensor isotropic. Furthermore, because the coupling between surface and bodywaves is negligible at large depth, its magnitude tends to the constant value D bulkb = c τ b → b /
3, as expected on physicalgrounds. In other words, the diffusion process at depth is governed by a simple 3-D diffusion equation for body waves withdiffusion constant D bulkb . This in turn suggests that at long lapse-time, the coda should decay as t − / in a non-absorbingmedium. This point will be further substantiated by numerical simulations.In the vicinity of the surface z = O ( α − ), Eq. (71) shows that the diffusivity of coupled body and surface waves is bothdepth dependent and non-isotropic. The origin of the z -dependence is clear since the efficacy of the coupling between surfaceand body waves decays exponentially with depth. In the vicinity of the surface, the anisotropy stems from the transport of afraction 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 coefficientsdictated by the equipartition principle. The vertical diffusivity is -up to the (1 + R ( z )) pre-factor inherited from the change ofscale in the vertical direction- equal to the diffusivity of body waves. In the next section, we illustrate the transport processof coupled body and surface waves by numerically simulating the system of Eq. (50). In this section, we explore some of the key features of our model with the aid of numerical simulations. The approach toequipartition as well as the role of mode coupling in the coda excitation are illustrated.
As outlined in introduction, Monte-Carlo simulations have been used for more than thirty years in seismology to simulate thetransport 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 modificationswhich 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 modeof propagation is randomly selected, following the source energy partitioning ratio (20), and the initial propagation directionis a uniformly distributed random vector in 2-D (resp. 3-D) for surface (resp. body) waves. Note that when the particle isof surface type, the particle propagates in a horizontal plane and its exact depth is immaterial. In fact, we may say thata particle of surface type is present at all depth with a probability distribution given by p s ( z ) = 2 α exp( − αz ) inheritedfrom the modal shape. The lapse-time to the first scattering event is randomly determined and obeys a simple exponentialdistribution when the particle represents surface waves. In the case of body waves, the selection process is more complicatedbecause their scattering mean free time depends on the depth in the medium. To address this difficulty, we employ the methodof 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 → b ) / (1 /l s ) may be interpreted as the transition probability froma surface to a body wave mode. Note that when such an event occurs, the particle is reinjected at a random depth in themedium following the probability distribution p s ( z ). To obtain energy envelopes, the position and mode of the particle ismonitored on a cylindrical grid at regular time intervals. The local energy density is estimated by averaging the number ofparticles per cell over a sufficiently large number of random walks. For accuracy, it is important that the cells be relativelysmall compared to the shortest mean free path in the medium. Figure 3 illustrates the striking difference between the global and local partitioning of the seismic energy into surface and bodywaves. The following parameters have been employed in the simulation: α = 1km − , c = 3km/s and τ s → s = 20s, τ s → b = 30s, τ b → b = 30s, τ b → s ( z ) = 10 exp(2 z )s. Note that in our model the group velocity of surface waves v g ≈ . z = 1km)and a deep one ( z = 5km), 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 ¯¯ E s (see the remarks L. Margerin et al. L. Margerin et al.
Time ¯¯ E s / ¯ E b ShallowDeepAsymptote
Time ¯¯ E s / ¯¯ E b ShallowDeep / t / Figure 2.
Local and global evolution of energy partitioning for a shallow (solid line) and a deep (dash-dotted line) source in a heterogeneous half-space filled with point scatterers. The horizontal time scaleis the mean free time for surface wave scattering. The penetration depth of surface waves is ↵ = 1km and the angular frequency is ! = 2 ⇡ Hz. The shallow and deep source are located, respectively,at depth ↵ and 5 ↵ . (Left): Temporal evolution of the ratio between the total energy of surfacewaves ( ¯¯ E s ) and the energy density of body waves integrated over a horizontal plane ( ¯ E b ). The bodywave energy is evaluated at the surface and averaged over a depth z = 1km. The asymptote (left) isthe prediction of Eq. (25). (Right) Temporal evolution of the ratio between the total energy of surface( ¯¯ E s ) and body waves ( ¯¯ E b ). The dotted line shows an algebraic t / decay. uniformly distributed random vector in 2-D (resp. 3-D) for surface (resp. body) waves. Notethat when the particle is of surface type, the particle propagates in a horizontal plane and itsexact depth is immaterial. In fact, we may say that a particle of surface type is present at alldepth with a probability distribution given by p s ( z ) = 2 ↵ exp( ↵z ) inherited from the modalshape. The lapse-time to the first scattering event is randomly determined and obeys a simple Depth . . . . . . . E n e r g y D e n s i t y t = 1 t = 2 t = 5 t = 10 t = 15 Figure 3.
Horizontally-averaged energy density of body waves ¯ E b as a function of depth in a heteroge-neous half-space filled point scatterers. The energy is averaged over a depth range z = ↵ / . ↵ . The solid and dashed line correspond toa shallow source ( z = 1km) and a deep source ( z = 5km), respectively. The lapse time in the codain mean free time unit is indicated on the right of the corresponding curves. Figure 3.
Local and global evolution of energy partitioning for a shallow (solid line) and a deep (dash-dotted line) source in a heteroge-neous half-space filled with point scatterers. The horizontal time scale is the mean free time for surface wave scattering. The penetrationdepth of surface waves is α − = 1 km and the angular frequency is ω = 2 π Hz. The shallow and deep source are located, respectively,at depth α − and 5 α − . (Left): Temporal evolution of the ratio between the total energy of surface waves ( ¯¯ E s ) and the energy densityof body waves integrated over a horizontal plane ( ¯ E b ). The body wave energy is evaluated at the surface and averaged over a depth∆ z = 1km. The asymptote (left) is the prediction of Eq. (25). (Right) Temporal evolution of the ratio between the total energy of surface( ¯¯ E s ) and body waves ( ¯¯ E b ). The dotted line shows an algebraic t − / decay. before Eq. 56 for a reminder of the notations), and the horizontally-integrated energy density of body waves ¯ E b at the surface z = 0, averaged over a depth range ∆ z = 1km. Hence, the ratio ¯¯ E s / ¯ E b has unit of inverse length. Independent of the sourcedepth, we find that the partitioning of the energy density at the surface -into surfacic energy of surface waves and volumetricenergy of body waves- converges toward the prediction of equipartition theory, at long lapse-time in the coda (see Eq. 25 and60). This numerical result confirms that the analysis of equipartition given in the previous section in slab geometry extends tothe half-space geometry. Furthermore, we find that the surface-to-body energy ratio overshoots the prediction of equipartitiontheory for the two sources at short lapse-time, by a factor which decreases with the source depth z .The stabilization of the local energy density ratio of surface and body waves at the surface of the half-space is to becontrasted 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 Figurealso suggests that the global transfer of energy from surface waves to body waves occurs at a rate proportional to t − / atlong lapse-time.Further insight into the equipartition process is offered in Figure 4, where we show the depth dependence of thehorizontally-averaged body wave energy density at different lapse-time in the coda. All the parameters of the simulationare the same as in Figure 3, except for the much finer spatial resolution ∆ z = α − / . Depth . . . . . . . E n e r g y D e n s i t y ←− t = 1 ←− t = 2 ←− t = 5 ←− t = 10 ←− t = 15 Figure 4.
Horizontally-averaged energy density of body waves ¯ E b as a function of depth in a heterogeneous half-space filled with pointscatterers. The energy is averaged over a depth range ∆ z = α − / . α − .The solid and dashed lines correspond to a shallow source ( z = 1km) and a deep source ( z = 5km), respectively. The lapse time in thecoda in mean free time unit is indicated on the right of the corresponding curves. adiative transfer of body and surface waves E n e r g y D e n s i t y E n e r g y D e n s i t y Figure 5.
Snapshots of the energy density of surface waves (cid:15) s (left) and of the volumetric energy density of body waves E b (right) atthe surface of a heterogeneous half-space filled with point scatterers in the case of a shallow source ( z = α − ). The horizontal axes arein units of the scattering mean free path of surface waves (left, top), body waves (right, top) and in kms (bottom). For body waves, wetake the mean free path value in the bulk of the medium ( z → ∞ ). The energy is averaged over cylindrical shells of width ∆ r = 5kmand deph ∆ z = 5km. roughly 10 mean free times, the depth distribution of body wave energy becomes homogeneous over a depth range at leastas large as 10 α − , independent of the source depth. This simulation therefore confirms the theoretical analysis performed inslab geometry. The homogenization of the energy of body waves is a dynamic process: the energy density of surface wavesincreases exponentially near the surface, thereby generating a larger amount of body-wave converted energy; this processis compensated by the exponential increase of the conversion rate from body to surface waves, which eventually yields anequilibrium. Note that the total energy density does not homogenize with depth, due to the exponential decay of the surfacewave eigenfunction with depth.In Figure 5 and 7, we illustrate in greater details the multiple-scattering process by showing snapshots of the surfacic andvolumetric energy densities ε s ( t, r ) and E b ( t, r , z ) | z =0 at regular time intervals ∆ t = 1 τ s starting at a lapse time t = 0 . τ s fora shallow ( z = 1km) and a deep ( z = 5km) source, respectively. The scattering parameters are the same as in Figure 3 andthe energy is averaged over a range of epicentral distance ∆ r = 5km and depth ∆ z = 5km. We use a double horizontal axison Figures 5-7 to show simultaneously the epicentral distance in kms and in units of mean free path. Note that in the case of E n e r g y D e n s i t y E n e r g y D e n s i t y Figure 6.
Same as Figure 5 but the coupling between surface and body waves has been deactivated. See text for further explanations. L. Margerin et al. E n e r g y D e n s i t y E n e r g y D e n s i t y Figure 7.
Snapshots of the energy density of surface waves (cid:15) s (left) and of the volumetric energy density of body waves E b (right) atthe surface of a heterogeneous half-space filled with point scatterers in the case of a deep source ( z = 5 α − ). The horizontal axes arein units of the scattering mean free path of surface waves (left, top), body waves (right, top) and in kms (bottom). For body waves, wetake the mean free path value in the bulk of the medium ( z → ∞ ). The energy is averaged over cylindrical shells of width ∆ r = 5kmand depth ∆ z = 5km. body waves, we take the value of the mean free path in the bulk of the medium τ b → 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 ofsurface waves, this amounts to computing the solution of a conventional 2-D multiple-scattering process with the mean freetime τ s . In the case of body waves, we consider a conventional 3-D multiple-scattering process in a half-space with a constantmean free time τ b → b , i.e., we remove the boundary layer where the coupling with surface waves occurs. To facilitate thecomparison between Figure 5 and 6, we have adjusted the strength of the source term in the conventional multiple-scattteringsimulations so that they match exactly the energy released at the source in the form of body and shear waves in the coupledcase.We first analyze the transport of surface waves in the case of a shallow source. As compared to the conventional 2-Dcase, mode coupling has at least two visible effects on the spatial distribution of the surface wave energy. First, it lowers theenergy level in the coda. As an illustration, we observe that after 10 mean free times the coda intensity is reduced by a factor S c a tt e r i n g O r d e r S c a tt e r i n g O r d e r Figure 8.
Contribution of the different orders of scattering to the energy envelopes of surface waves shown in Figure 5 (left, solidlines) and 7 (right, solid lines). Snapshots of the mean scattering order are represented as a function of the distance from the source. Forreference, the dashed line show the same distribution for a conventional transport model in 2-D. The horizontal axes are in units of thescattering mean free path of surface waves (top) and in kms (bottom). adiative transfer of body and surface waves Time10 E n e r g y D e n s i t y ShallowDeep t t t Time10 E n e r g y D e n s i t y ShallowDeep t t t Figure 9.
Energy density of body waves E b (left) and surface waves (cid:15) s (right) at the surface of a heterogeneous half-space filled withpoint scatterers in the case of a shallow source ( z = α − ). The enegy is averaged over a depth ∆ z = 5km and an epicentral distancerange ∆ r = 5km. The station is located at an epicentral distance of 50km. The horizontal axis is in units of the scattering mean freetime of surface waves in logarithmic scale. Typical algebraic decays are also shown. at least equal to 10 in Figure 5 compared to Figure 6. Second, mode coupling appears to enhance the visibility of ballisticsurface waves. In Figure 6, the ballistic term is completely masked by the diffuse contribution at roughly 6 mean free pathfrom the source, whereas a small ballistic peak is still visible at roughly 10 mean free paths from the source in Figure 5. Itis worth noting that the ballistic contribution is exactly identical in Figures 5 (left) and 6 (left). Again, this is the strongdecrease of the energy of scattered coda waves which explains the difference between the two Figures. Figure 8, which displaysthe spatial distribution of the mean order of scattering in the coda, reveals that the coda of coupled surface waves is depletedin high-order multiply-scattered waves compared to a conventional 2-D transport process. In other words, mode conversionsentail a strong conversion of multiply-scattered surface waves into body waves which decreases the energy level of surfacewave coda and, by comparison, enhances the ballistic contribution. Examination of Figures 5 (right) and 6 (right) revealsthat the effects of mode coupling on body waves are opposite to the ones just described for surface waves. Thus, we observethat the energy level in the coda is slightly increased by the transfer energy from surface wave to body waves. An additionalcontribution comes from the increase of the scattering strength of body waves near the surface which attenuates the ballisticwaves and transfers their energy into the coda. Examination of the decay of the ballistic peak of body waves with epicentraldistance in Figures 5 and 6 confirms the increased attenuation entailed by the coupling with surface waves. Other more exoticphenomena are also visible in Figure 5 such as some precursory body waves arrival due to the coupling from surface wavesto body waves. However this process is a very peculiar feature of our model, due to the higher wavespeed of surface wavescompared to the one of body waves.Further differences between our coupled model for surface and body waves and conventional transport theory is illustratedin Figure 7 where we show snapshots of the energy distribution of surface and body waves in the case of a deep source. Notethat in that case, surface waves can only be generated by mode conversions so that ballistic arrivals are absent in Figure7 (left). Interestingly, our numerical simulations indicate that surface waves are rapidly excited to a non-negligible level inthe coda. Examination of Figure 8 (right) further indicates that multiple-scattering is at the origin of the generation ofsurface waves in the coda when the source is located at large depth. These observations agree with our theoretical analysis ofequipartition, which implies that, independent of the source depth, the coda at the surface of a half-space always appears asa 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 = 1km)and a deep source ( z = 5km) at an epicentral distance of 50km. The scattering parameters are the same as in all previousFigures and the spatial resolution of the computation is 5km again. The impact of the depth of the source on the excitation ofballistic waves is obvious in Figure 9 and confirms the analysis of Figures 5 and 7. In particular, it is apparent that the directbody waves are less attenuated in the case of the deep source, as a consequence of the exponential decay of the scatteringconversions from body to surface waves with depth. To facilitate the identification of different propagation regimes in thecoda, we have superposed on the graphs some typical algebraic decays: t − for scattering in 2-D (either single or multiple, seee.g. Paasschens 1997), t − / for multiple scattering in 3-D, and t − for single-scattering in 3-D. In Figure 9, we observe thatfor both body and surface waves the coda obeys a t − / decay law at long lapse-time, independent of the source depth, whichis characteristic of a 3-D diffusion process. This supports the predictions of the diffusion model and confirms the dominanceof body waves in the transport process at large lapse-time. At short lapse-time, we observe a distinct behavior between the L. Margerin et al. two kinds of waves, particularly in the case of a shallow source. After the passage of the ballistic waves, body waves appear todecay slightly more slowly than the asymptotic t − / behavior. This may reflect the conversion of surface waves to body wavesas discussed in the analysis of Figure 5. Two propagation regimes show up clearly on the surface wave energy envelope, witha transition between the two around a lapse-time of 10 mean free times. At short time, the decay of surface waves appearsto be faster than the one of body waves, probably as a consequence of the transfer of surface wave energy into the volume asdiscussed in relation with Figure 5. Taken together, Figures 5-9 illustrate the much richer behavior of the coupled system oftransport Eq. (50), compared to the conventional transport process without coupling between surface and body waves. This work represents a first attempt at formulating a self-consistent theory of RT of seismic waves in a half-space geometryincluding the coupling between surface and body waves. The main approximation underlying our work is that, upon reflectionat the surface, the phase of body waves is randomized so that upgoing and downgoing fluxes may be considered as statisticallyindependent. Our approach distinguishes itself from the standard Eqs of RT for scalar waves found in the literature in oneimportant way: it keeps track - to some extent - of the wave behavior in the vicinity of the surface. This has a number ofconsequences: (1) surface and body waves are coupled by conversion scattering (2) even in a statistically homogeneous mediumit requires that the scattering properties of body waves depend on the depth in the medium. Furthermore, the reciprocityrelation between the surface-to-body vs body-to-surface mean free times plays a prominent role in the establishment of anequipartition 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 scatteringproperties and wave velocities of body and surface waves. We also show that our RT Eqs are operational, in the sense thatthey are readily amenable to numerical solutions by Monte-Carlo simulations. These simulations could be used in the futureto study in more details the dynamics of equipartition, in particular, how the equipartition time varies as a function of theratio 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, weplan to address the following issues: (1) Evaluate the impact of neglecting the interference between upgoing and downgoingbody waves on the scattering cross-section and, if possible, go beyond this approximation. (2) Extend the theory to morerealistic finite size scatterers and more general spatial distributions of scatterers. (3) Incorporate polarization effects for elasticwaves at a free surface. (4) Absorption of energy is also a very important mechanism of attenuation, which has been entirelyneglected in this work for simplicity. Because the sub-surface of the Earth is thought to be very strongly attenuating due tothe 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 ofcoupled S and P waves. Special efforts should be devoted to this important topic before our formalism can be applied to realseismic data. ACKNOWLEDGMENTS
The authors wish to thank the Associate Editor S. Ni and an anonymous referee for their suggestions to clarify the presentationof the results. The careful comments and constructive criticisms of H. Sato contributed to significant changes and improvementsin the content of the manuscript. The authors acknowledge the European Research Council under the European Union Horizon2020 research and innovation program (grant agreement no. 742335 - F-IMAGE).
REFERENCES
Abubakirov, I. & Gusev, A., 1990. Estimation of scattering properties of lithosphere of Kamchatka based on Monte-Carlo simulationof record envelope of a near earthquake,
Physics of the Earth and Planetary Interiors , (1), 52–67.Aki, K. & Richards, P. G., 2002. Quantitative Seismology, 2d Ed. , University Science Book.Akkermans, E. & Montambaux, G., 2007.
Mesoscopic Physics of Electrons and Photons , Cambridge University press, Cambridge.Carcol´e, E. & Sato, H., 2010. Spatial distribution of scattering loss and intrinsic absorption of short-period s waves in the lithosphereof japan on the basis of the multiple lapse time window analysis of hi-net data,
Geophysical Journal International , (1), 268–290.Eulenfeld, T. & Wegler, U., 2017. Crustal intrinsic and scattering attenuation of high-frequency shear waves in the contiguous unitedstates, Journal of Geophysical Research: Solid Earth , (6), 4676–4690.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, Geophysical Journal International , , 787–800.Gaebler, P. J., Sens-Sch¨onfelder, C., & Korn, M., 2015. The influence of crustal scattering on translational and rotational motions inregional and teleseismic coda waves, Geophysical Journal International , (1), 355–371.Gelfand, I. & Fomin, S., 1963. Calculus of Variations , Prentice Hall, New Jersey. adiative transfer of body and surface waves Gusev, A. & Abubakirov, I., 1996. Simulated envelopes of non-isotropically scattered body waves as compared to observed ones: anothermanifestation of fractal heterogeneity,
Geophysical Journal International , , 49–60.Hein Hoernig, R. O., 2010. Green’s functions and integral equations for the Laplace and Helmholtz operators in impedance half-spaces ,Theses, Ecole Polytechnique X.Hennino, R., Tr´egour`es, N., Shapiro, N. M., Margerin, L., Campillo, M., van Tiggelen, B. A., & Weaver, R. L., 2001. Observation ofequipartition of seismic waves,
Phys. Rev. Lett. , , 3447–3450.Hoshiba, M., 1993. Separation of scattering attenuation and intrinsic absorption in japan using the multiple lapse time window analysisof full seismogram envelope, Journal of geophysical research , (B9), 15809–15824.Hoshiba, M., 1995. Estimation of nonisotropic scattering in western Japan using coda wave envelopes: Application of a multiplenonisotropic scattering model, Journal of Geophysical Research , (B1), 645–657.Hoshiba, M., 1997. Seismic coda wave envelope in depth-dependent S wave velocity structure, Physics of the Earth and PlanetaryInteriors , , 15–22.Jing, Y., Zeng, Y., & Lin, G., 2014. High-frequency seismogram envelope inversion using a multiple nonisotropic scattering model:Application to aftershocks of the 2008 wells earthquake, Bulletin of the Seismological Society of America , (2), 823–839.Kinoshita, S., 1993. Evaluation of site factor and propagation characteristics by means of earthquake observation, Zisin: Journal of theSeismological Society of Japan , (2), 161–170.Kittel, C. et al., 1976. Introduction to solid state physics , vol. 8, Wiley New York.Lagendijk, A. & Van Tiggelen, B. A., 1996. Resonant multiple scattering of light,
Physics Reports , (3), 143–215.Lux, I. & Koblinger, L., 1991. Monte Carlo Particle Transport Problems: Neutron and Photon Calculations , CRC Press, Boca Raton,Florida.Maeda, T., Sato, H., & Nishimura, T., 2008. Synthesis of coda wave envelopes in randomly inhomogeneous elastic media in a half-space:single scattering model including rayleigh waves,
Geophysical Journal International , (1), 130–154.Mancinelli, N., Shearer, P., & Liu, Q., 2016a. Constraints on the heterogeneity spectrum of earth’s upper mantle, Journal of GeophysicalResearch: Solid Earth , (5), 3703–3721.Mancinelli, N., Shearer, P., & Thomas, C., 2016b. On the frequency dependence and spatial coherence of pkp precursor amplitudes, Journal of Geophysical Research: Solid Earth , (3), 1873–1889.Mancinelli, N. J. & Shearer, P. M., 2013. Reconciling discrepancies among estimates of small-scale mantle heterogeneity from pkpprecursors, Geophysical Journal International , (3), 1721–1729.Margerin, L., 2005. Introduction to radiative transfer of seismic waves , vol. 157 of
Geophysical monograph , chap. Seismic Earth:Analysis of broadband seismograms, pp. 229–252, American Geophysical Union, Washington.Margerin, L. & Nolet, G., 2003. Multiple scattering of high-frequency seismic waves in the deep earth: PKP precursor analysis andinversion for mantle granularity,
J. Geophys. Res. , , 2514.Margerin, L., Campillo, M., & van Tiggelen, B. A., 1998. Radiative tranfer and diffusion of waves in a layered medium : new insightinto coda Q, , 596–612.Margerin, L., Campillo, M., & van Tiggelen, B. A., 2000. Monte Carlo simulation of multiple scaterring of elastic waves, J. Geophys.Res. , , 7873–7892.Margerin, L., Van Tiggelen, B., & Campillo, M., 2001. Effect of absorption on energy partition of elastic waves in the seismic coda, Bulletin of the Seismological Society of America , (3), 624–627.Margerin, L., Plan`es, T., Mayor, J., & Calvet, M., 2016. Sensitivity kernels for coda-wave interferometry and scattering tomography:theory and numerical evaluation in two-dimensional anisotropically scattering media, Geophysical Journal International , , 650–666.Morse, P. M. & Ingard, K. U., 1986. Theoretical acoustics , Princeton university press.Obermann, A., Plan`es, T., Larose, E., Sens-Sch¨onfelder, C., & Campillo, M., 2013. Depth sensitivity of seismic coda waves to velocityperturbations in an elastic heterogeneous medium,
Geophysical Journal International , (1), 372–382.Obermann, A., Plan`es, T., Hadziioannou, C., & Campillo, M., 2016. Lapse-time-dependent coda-wave depth sensitivity to local velocityperturbations in 3-d heterogeneous elastic media, Geophysical Journal International , (1), 59–66.Paasschens, J., 1997. Solution of the time-dependent boltzmann equation, Phys. Rev. E , , 1135–1141.Pacheco, C. & Snieder, R., 2005. Time-lapse travel time change of multiply scattered acoustic waves, Acoustical Society of AmericaJournal , , 1300–1310.Pacheco, C. & Snieder, R., 2006. Time-lapse traveltime change of singly scattered acoustic waves, Geophysical Journal International , (2), 485–500.Poupinet, G., Ellsworth, W., & Frechet, J., 1984. Monitoring velocity variations in the crust using earthquake doublets: an applicationto the Calaveras fault, California, J. Geophys. Res. , (B7), 5719–5731.Poupinet, G., Got, J.-L., & Brenguier, F., 2008. Earth Heterogeneity and Scattering Effects on Seismic Waves , vol. 50 of
Advancesin Geophysics , chap. Monitoring temporal variations of physical properties in the crust by cross-correlating the waveforms of seismicdoublets, pp. 374–399, Academic Press, New York.Przybilla, J., Wegler, U., & Korn, M., 2009. Estimation of crustal scattering parameters with elastic radiative transfer theory,
GeophysicalJournal International , (2), 1105–1111.Ryzhik, L., Papanicolaou, G., & Keller, J. B., 1996. Transport equation for elastic and other waves in random media, Wave Motion , , 327–370.Sanborn, C. J. & Cormier, V. F., 2018. Modelling the blockage of Lg waves from three-dimensional variations in crustal structure, Geophysical Journal International , (2), 1426–1440.Sanborn, C. J., Cormier, V. F., & Fitzpatrick, M., 2017. Combined effects of deterministic and statistical structure on high-frequencyregional seismograms, Geophysical Journal International , (2), 1143–1159.Sato, H., 1994. Multiple isotropic scattering model including P-S conversion for the seismogram envelope formation, , 487–494.Sato, H. & Emoto, K., 2018. Synthesis of a scalar wavelet intensity propagating through von k´arm´an-type random media: Radiativetransfer theory using the born and phase-screen approximations, Geophysical Journal International , (2), 909–923.Sato, H., Fehler, M. C., & Maeda, T., 2012. Seismic Wave Propagation and Scattering in the Heterogeneous Earth , Springer Science& Business Media. L. Margerin et al.
Sens-Sch¨onfelder, C., Margerin, L., & Campillo, M., 2009. Laterally heterogeneous scattering explains Lg blockage in the Pyrenees,
JGR , , B07309.Shearer, P. & Earle, P., 2004. The global short-period wavefield modelled with a Monte Carlo seismic phonon method, GeophysicalJournal International , (3), 1103–1117.Sheng, P., 2006. Introduction to wave scattering, localization and mesoscopic phenomena , vol. 88, Springer Science & Business Media.Snieder, R., 1986. 3-D linearized scattering of surface waves and a formalism for surface wave holography,
Geophysical Journal Inter-national , (3), 581–605.Snieder, R., 2006. The theory of coda wave interferometry, Pure and Applied Geophysics , (2), 455–473.Takeuchi, N., Kawakatsu, H., Shiobara, H., Isse, T., Sugioka, H., Ito, A., & Utada, H., 2017. Determination of intrinsic attenuation inthe oceanic lithosphere-asthenosphere system, Science , (6370), 1593–1596.Tr´egour`es, N. P. & Van Tiggelen, B. A., 2002. Quasi-two-dimensional transfer of elastic waves, Physical Review E , (3), 036601.Turner, J. A. & Weaver, R. L., 1994. Radiative transfer of ultrasound, The Journal of the Acoustical Society of America , (6),3654–3674.Weaver, R. L., 1982. On diffuse waves in solid media, J. Acoust. Soc. America , , 1608–1609.Weaver, R. L., 1990. Diffusivity of ultrasound in polycrystals, Journal of the Mechanics and Physics of Solids , (1), 55–86.Wegler, U. & Sens-Schoenfelder, C., 2007. Fault zone monitoring with passive image interferometry, , 1029–1033.Wu, R. S., 1985. Multiple scattering and energy transfer of seismic waves—-separation of scattering effect from intrinsic attenuation I.Theoretical modelling, , 57–80.Yamamoto, M. & Sato, H., 2010. Multiple scattering and mode conversion revealed by an active seismic experiment at asama volcano,japan, Journal of Geophysical Research: Solid Earth , (B7).Yoshimoto, K., 2000. Monte Carlo simulation of seismogram envelopes in scattering media, Journal of Geophysical Research , (B3),6153–6161.Zeng, Y., 1993. Theory of scattered P-and S-wave energy in a random isotropic scattering medium, Bulletin of the Seismological Societyof America , (4), 1264–1276.Zeng, Y., 2006. Scattered surface wave energy in the seismic coda, Pure and Applied Geophysics , (2-3), 533–548. 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 theclassic book by Gelfand & Fomin (1963), after which our treatment is modeled. For simplicity we consider a vibrating stringof density ρ , tension T and length L . For the moment, we do not specify the boundary conditions. The total kinetic energystored in the string at time t is given by: T [ u ]( t ) = (cid:90) L ρ ( ∂ t u ( x , t )) dx, (A.1)where u denotes the displacement field. The instantaneous potential energy stored in the string may be expressed as V [ u ]( t ) = (cid:90) L T ( ∂ x u ( x , t )) dx (A.2)According to Hamilton’s principle, among all possible displacement fields, the one that satisfies the actual equations of motionshould make the following action integral: I [ u ] = (cid:90) t t ( T − V )[ u ]( t ) dt (A.3)stationary. Mathematically, this principle of stationary action may be expressed as: ∂ (cid:15) I [ u + (cid:15)ψ ] | (cid:15) =0 = 0 (A.4)where ψ is an arbitrary function. This is sometimes written as δI = 0, where δI is known as the first variation of the actionintegral. Using integration by parts, Gelfand & Fomin (1963) establish that: δI = (cid:15) (cid:18)(cid:90) t t (cid:90) L ( − ρ∂ tt u ( x , t ) + T ∂ xx u ( x , t )) ψ ( x, t ) dxdt + T (cid:90) t t ( ∂ x u (0 , t ) ψ (0 , t ) − ∂ x u ( l, t ) ψ ( l, t )) dt (cid:19) (A.5)The arbitrariness of the function ψ in Eq. (A.5) implies both the governing wave equation for the vibrating string: ρ∂ tt u ( x, t ) − T ∂ xx u ( x, t ) = 0 (A.6)as well as the so-called natural boundary conditions: ∂ x u ( x, t ) | x =0 = 0 and ∂ x u ( x, t ) | x = l = 0 , (A.7) adiative transfer of body and surface waves which correspond to a string with free ends. In order to obtain mixed boundary conditions, it suffices to add to the potentialenergy (A.2), a term of the form χu (0 , t ) where χ is a constant. Eq. (A.5) must be modified accordingly by adding the newcontribution − (cid:15)χ (cid:82) t t u (0 , t ) ψ (0 , t ) dt which, in turn, implies a natural boundary condition of the mixed type at x = 0:( χ∂ x u ( x, t ) − T u ( x, t )) | x =0 = 0 (A.8)The total potential energy may be rewritten in integral form as follows: V [ u ]( t ) = (cid:90) L (cid:2) χu (0 , t ) δ ( x ) + T ( ∂ x u ( x , t )) (cid:3) dx (A.9)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 INA 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 twoparts and begin with the surface wave contribution: G s ( r , z, z ) = 2 αe − α ( z + z ) (2 π ) (cid:90) R e i k (cid:107) · r k + α − k (cid:107) + i(cid:15) d k (cid:107) (B.1)Introducing cylindrical coordinates ( k (cid:107) , φ ) and integrating over angle yields: G s ( r , z, z ) = 2 αe − α ( z + z ) π (cid:90) + ∞ J ( k (cid:107) r ) k + α − k (cid:107) + i(cid:15) dk (cid:107) (B.2)where J 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 (cid:107) axis using the Hankel function of the first kind instead of the Besselfunction: G s ( r , z, z ) = αe − α ( z + z ) π (cid:90) + ∞−∞ H (1)0 ( k (cid:107) r ) k + α − k (cid:107) + i(cid:15) dk (cid:107) (B.3)In the last step, we employ the residue theorem by closing the contour in the upper half of the complex plane with a semi-circleof radius R → + ∞ and note the presence of pole at k (cid:107) = √ k + α + iη , ( η → + ). Thanks to the exponential decay of theintegrand, the integral on the semi-circle vanishes which yields: G s ( r , z, z ) = − iα π e − α ( z + z ) H (1)0 ( (cid:112) k + α r ) . (B.4)The result (B.4) is exact. The far-field approximation (16) follows by application of standard asymptotic expansions to theHankel function.The computation of the body wave contribution can also be split into two parts: G b ( r , z, z ) = 1(2 π ) (cid:90) + ∞ dq (cid:90) R e i k (cid:107) · r ( e − iqz + r ( q ) e iqz )( e − iqz + r ( q ) e iqz ) ∗ k − k (cid:107) − q + i(cid:15) d k (cid:107) = 1(2 π ) (cid:90) + ∞−∞ dq (cid:90) R e iq ( z − z ) e i k (cid:107) · r k − k (cid:107) − q + i(cid:15) d k (cid:107) + 1(2 π ) (cid:90) + ∞−∞ dq (cid:90) R r ( q ) e iq ( z + z ) e i k (cid:107) · r k − k (cid:107) − q + i(cid:15) d k (cid:107) = G ∞ ( r , z, z ) + G r ( r , z, z ) (B.5)where the unitarity of the reflection coefficient has been used and the q integral has been extended from −∞ to + ∞ thanksto the relation r ( q ) ∗ = r ( − q ). The first term in the second equality of (B.5) may be recognized as the full-space solution tothe Helmholtz Eq.: G ∞ ( r , z, z ) = 1(2 π ) (cid:90) + ∞−∞ dq (cid:90) R e iq ( z − z ) e i k (cid:107) · r k − k (cid:107) − q + i(cid:15) d k (cid:107) = − e ikR πR , (B.6)where R = (cid:112) r + ( z − z ) . The second term in the second equality of (B.5) represents the waves reflected at the surface: G r ( r , z, z ) = 1(2 π ) (cid:90) + ∞−∞ dq (cid:90) R r ( q ) e iq ( z + z ) e i k (cid:107) · r k − k (cid:107) − q + i(cid:15) d k (cid:107) (B.7) L. Margerin et al.
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: G r ( r , z, z ) = − i π (cid:90) + ∞−∞ r ( q ) e iq ( z + z ) H (1)0 ( (cid:112) k − q r ) dq, (B.8)To approximate this last integral in the far-field of the source, we first remark that for | q | > k the cylindrical waves areevanescent so that we may legitimately take − k and + k as integration limits. We next make use of the far-field expansion ofthe Hankel function to obtain the following oscillatory integral representation: G r ( r , z, z ) ≈ − i π (cid:114) πr (cid:90) + k − k r ( q ) e iq ( z + z )+ i √ k − q r − iπ/ (cid:112) k − q dq (B.9)Further noting that the derivative of the phase term: φ ( q ) = q ( z + z ) + (cid:112) k − q r (B.10)vanishes at : q = k ( z + z ) R (cid:48) (B.11)with R (cid:48) = (cid:112) r + ( z + z ) , we apply the stationary phase formula to obtain after some straightforward algebra: G r ( r , z, z ) ≈ − r ( q ) e ikR (cid:48) πR (cid:48) . (B.12)This term may be interpreted as the contribution of the image point of the source with a strength given by the reflectioncoefficient evaluated at an incidence angle corresponding to the specularly reflected ray connecting the source to the detectionpoint (see Eq. B.11). To complete the far-field approximation, we first note the following expansions: R = R − z z/R + o (1 /R ), R (cid:48) = R + z z/R + o (1 /R ) where R = √ r + z . Neglecting all terms smaller than 1 /R for the amplitude, all terms smallerthan z /R for the phase and further approximating q as kz/Rkz/R