Grégoire Heller
email: gregoire.heller@irap.omp.eu

Ludovic Margerin
email: ludovic.margerin@irap.omp.eu

Olivier Sèbe
email: olivier.sebe@cea.fr

Jessie Mayor
email: jessie.mayor@edf.fr

Marie Calvet
email: marie.calvet@irap.omp.eu

Paul Sabatier

Revisiting multiple-scattering principles in a crustal waveguide: equipartition, depolarization and coda normalization

Keywords: coda waves, polarized waves, crustal waveguide, radiative transfer, scattering, energy partition, seismic source

Introduction

Coda waves are scattered waves by random small-scale heterogeneities of the medium [START_REF] Aki | Origin of coda waves: source, attenuation, and scattering effects[END_REF] and contain valuable information on the fluctuations of elastic and anelastic parameters in the Earth (see [START_REF] Sato | Seismic wave propagation and scattering in the heterogeneous earth[END_REF]Sato, 2019, for a review). Numerous studies extracted fundamental propagation properties from coda waves, such as scattering and intrinsic attenuation [START_REF] Wu | Multiple scattering and energy transfer of seismic waves-separation of scattering effect from intrinsic attenuation-i. theoretical modelling[END_REF][START_REF] Fehler | Numerical basis of the separation of scattering and intrinsic absorption from full seismogram envelope. a monte-carlo simulation of multiple isotropic scattering[END_REF][START_REF] Lacombe | Separation of intrinsic absorption and scattering attenuation from lg coda decay in central france using acoustic radiative transfer theory[END_REF][START_REF] Sens-Schönfelder | Radiative transfer theory for estimation of the seismic moment[END_REF][START_REF] Przybilla | Estimation of crustal scattering parameters with elastic radiative transfer theory[END_REF][START_REF] Eulenfeld | Measurement of intrinsic and scattering attenuation of shear waves in two sedimentary basins and comparison to crystalline sites in germany[END_REF][START_REF] Gaebler | Seismic scattering and absorption parameters in the w-bohemia/vogtland region from elastic and acoustic radiative transfer theory[END_REF][START_REF] Mayor | Tomography of crustal seismic attenuation in metropolitan france: Implications for seismicity analysis[END_REF], while others focused on the extraction of source properties from these multi-scattered waves [START_REF] Mayeda | Moment, energy, stress drop, and source spectra of western united states earthquakes from regional coda envelopes[END_REF][START_REF] Mayeda | Stable and transportable regional magnitudes based on coda-derived moment-rate spectra[END_REF][START_REF] Sens-Schönfelder | Radiative transfer theory for estimation of the seismic moment[END_REF][START_REF] Denieul | M w estimation from crustal coda waves recorded on analog seismograms[END_REF][START_REF] Sèbe | Spectral factorization of the source time function of an earthquake from coda waves, application to the 2003 rambervillers, france, earthquake[END_REF]. However, the separation of source, propagation and site effects can be difficult, as a consequence of our incomplete knowledge of the small-scale physical properties of the Earth's interior in addition to oversimplified models or assumptions.

There are different approaches to model energy propagation in multiple scattering media containing small-scale heterogeneities. In seismology, full wavefield simulations have been fruitfully employed [START_REF] Frankel | Finite difference simulations of seismic scattering: Implications for the propagation of short-period seismic waves in the crust and models of crustal heterogeneity[END_REF][START_REF] Obermann | Lapsetime-dependent coda-wave depth sensitivity to local velocity perturbations in 3-d heterogeneous elastic media[END_REF][START_REF] Leng | 3-d scattering of elastic waves by small-scale heterogeneities in the earth's mantle[END_REF]. Here, we adopt a radiative transfer approach which has been widely used in the past and is a common tool across disciplines.

Radiative transfer theory has been initially developed for electromagnetic waves [START_REF] Chandrasekhar | Radiative Transfer[END_REF], and then, transposed to seismology and acoustics [START_REF] Wu | Multiple scattering and energy transfer of seismic waves-separation of scattering effect from intrinsic attenuation-i. theoretical modelling[END_REF][START_REF] Sato | Multiple isotropic scattering model including ps conversions for the seismogram envelope formation[END_REF][START_REF] Weaver | Diffusivity of ultrasound in polycrystals[END_REF][START_REF] Ryzhik | Transport equations for elastic and other waves in random media[END_REF]. The Monte-Carlo method is often used to solve numerically the radiative transfer equations. Starting from a simplified unbounded heterogeneous medium in a scalar wave approximation [START_REF] Gusev | Monte-carlo simulation of record envelope of a near earthquake[END_REF][START_REF] Hoshiba | Simulation of multiple-scattered coda wave excitation based on the energy conservation law[END_REF], several improvements of the numerical models have been made in order to consider more realistic Earth models. Gradually more complex deterministic structures were considered, including depth-dependent velocity profiles consisting of several layers [START_REF] Hoshiba | Seismic coda wave envelope in depth-dependent s wave velocity structure[END_REF][START_REF] Margerin | Radiative transfer and diffusion of waves in a layered medium: new insight into coda q[END_REF][START_REF] Lacombe | Separation of intrinsic absorption and scattering attenuation from lg coda decay in central france using acoustic radiative transfer theory[END_REF], velocity gradients (Yoshimoto, 2000) or a 3D velocity model with a variable crustal thickness [START_REF] Sanborn | Modelling the blockage of lg waves from three-dimensional variations in crustal structure[END_REF]. These studies showed that deterministic structure have great impact on the coda amplitude and shape, which is controlled by partial guiding and leakage. The neglect of these effects could lead to a biased estimation of source or attenuation parameters, as shown in [START_REF] Margerin | Radiative transfer and diffusion of waves in a layered medium: new insight into coda q[END_REF]; [START_REF] Bianco | Separation of depth-dependent intrinsic and scattering seismic attenuation in the northeastern sector of the italian peninsula[END_REF].

Depth-dependent scattering and intrinsic attenuation were first implemented by [START_REF] Hoshiba | Simulation of coda wave envelope in depth dependent scattering and absorption structure[END_REF], who also pointed out some possible limits of the coda normalization method. More recently, lateral variations of scattering or intrinsic attenuation parameters were also taken into account [START_REF] Sens-Schönfelder | Laterally heterogeneous scattering explains lg blockage in the pyrenees[END_REF][START_REF] Sanborn | Modelling the blockage of lg waves from three-dimensional variations in crustal structure[END_REF], showing that a zone of strong attenuation affects considerably the coda shape (especially the Lg coda which can be greatly attenuated).

An attractive feature of radiative transfer theory is that it allows one to keep track of mode conversions and shear wave polarization in heterogeneous elastic me-dia. The theory was established by [START_REF] Weaver | Diffusivity of ultrasound in polycrystals[END_REF] and the first seismological applications were developed by Zeng (1993); [START_REF] Sato | Multiple isotropic scattering model including ps conversions for the seismogram envelope formation[END_REF], including the P -to-S and S-to-P conversions by scattering. In later numerical works, a complete description of polarization was included through a 5-dimensional Stokes vector [START_REF] Turner | Radiative transfer and multiple scattering of diffuse ultrasound in polycrystalline media[END_REF][START_REF] Margerin | Monte carlo simulation of multiple scattering of elastic waves[END_REF]. Such a description of shear waves allows one to separate the contributions of SV and SH polarizations which behave differently at both internal interface discontinuities and at the free surface. [START_REF] Turner | Ultrasonic radiative transfer in polycrystalline media: Effects of a fluid-solid interface[END_REF] developed a full treatment of wave polarization in the case of incident plan waves in a scattering medium limited by a solid-liquid interface. [START_REF] Sens-Schönfelder | Laterally heterogeneous scattering explains lg blockage in the pyrenees[END_REF] and [START_REF] Sanborn | Combined effects of deterministic and statistical structure on high-frequency regional seismograms[END_REF] considered the contributions of SV and SH polarizations in addition to deterministic reflections in the case of a crustal waveguide. To our understanding, these two studies neglected the ellipticity aspect of wave polarization which develops when S-waves are incident on a boundary at post-critical incidence angles. Since polarization affects the scattering pattern of S-waves, a complex feedback can be expected between volume scattering and deterministic reflections. The initial polarization state is controlled by the source radiation patterns, which have been included into Monte-Carlo simulations in order to simulate realistic shear dislocations such as small earthquakes [START_REF] Rachman | Separation of intrinsic and scattering attenuation using single event source in south korea[END_REF][START_REF] Sanborn | Combined effects of deterministic and statistical structure on high-frequency regional seismograms[END_REF]. The source radiation influences the shape of the seismic signal, especially around the direct arrivals, where the wavefield is not yet randomized. However, it is well accepted that there is no signature of the source mechanism on the late coda excitation.

The literature referenced above amply illustrates that both deterministic structures and random heterogeneities can greatly influence the coda shape and its polarization state. Yet, the complex interactions between bulk scattering and reflection/transmission at interfaces have not yet been systematically explored. The pur-pose of our study is to fill this gap in our understanding of the coda, with a particular view on some well established wavefield properties in infinite scattering media: equipartition and coda normalization. The use of these concepts in seismology is summarized hereafter.

The principle of equipartition is well understood in unbounded random media, where all directions of propagation become equiprobable after several scattering events and the partition of body wave energy among different polarizations tends to a fixed value which is independent of the source and scattering properties [START_REF] Weaver | On diffuse waves in solid media[END_REF][START_REF] Ryzhik | Transport equations for elastic and other waves in random media[END_REF]. This concept is fundamental for some domains of seismology, such as coda wave interferometry [START_REF] Snieder | Coda wave interferometry and the equilibration of energy in elastic media[END_REF] used for imaging or monitoring (Sens-Schönfelder and Wegler, 2006a) or Green's function retrieval from coda [START_REF] Campillo | Long-range correlations in the diffuse seismic coda[END_REF] or ambient noise correlations [START_REF] Shapiro | Emergence of broadband rayleigh waves from correlations of the ambient seismic noise[END_REF].

In seismology, the equipartition principle has been supported by several numerical or experimental studies. [START_REF] Margerin | Monte carlo simulation of multiple scattering of elastic waves[END_REF] show numerically that energy equilibration is reached after a few scattering events in Monte-Carlo simulations. [START_REF] Souriau | Multimethod characterization of the french-pyrenean valley of bagnères-de-bigorre for seismic-hazard evaluation: observations and models[END_REF] shows that a stabilization of the H/V ratio is rapidly reached for coda waves in the Pyrenees. This indicates that the diffuse wavefield has reached a form of equilibrium akin to the predictions of equipartition theory. Using receiver arrays, [START_REF] Campillo | Seismic wave diffusion in the earth lithosphere[END_REF]; [START_REF] Shapiro | The energy partitioning and the diffusive character of the seismic coda[END_REF]; [START_REF] Hennino | Observation of equipartition of seismic waves[END_REF]; [START_REF] Margerin | Energy partition of seismic coda waves in layered media: theory and application to pinyon flats observatory[END_REF] observed an equilibration of the ratio between compressional and shear deformation energies in the coda. Furthermore, these studies revealed that the P -to-S energy ratio can be different from the theoretical full-space value. In the case of [START_REF] Hennino | Observation of equipartition of seismic waves[END_REF], the discrepancy was convincingly reconciled with equipartition theory by incorporating Rayleigh waves at the surface of a half-space. But other studies confirmed analytically or numerically that the P -to-S energy ratio is shifted in the presence of intrinsic attenuation [START_REF] Margerin | Effect of absorption on energy partition of elastic waves in the seismic coda[END_REF]Trégourès and van Tiggelen, 2002a), interfaces (Trégourès and van Tiggelen, 2002a) or sub-surface layering [START_REF] Margerin | Energy partition of seismic coda waves in layered media: theory and application to pinyon flats observatory[END_REF]. Based on a radiative transfer model for scalar waves in a wavguide geometry, [START_REF] Margerin | Breakdown of equipartition in diffuse fields caused by energy leakage[END_REF] shows that the energy flux at the free surface is never isotropic, even at long lapse-time. Hence, while a key concept, the conditions of application of equipartition to seismological data needs to be carefully examined.

The principle of coda normalization which originates from the work of [START_REF] Aki | Origin of coda waves: source, attenuation, and scattering effects[END_REF] is widely used to separate source, propagation and site properties [START_REF] Rautian | The use of the coda for determination of the earthquake source spectrum[END_REF][START_REF] Tsujiura | Spectral analysis of the coda waves from local earthquakes[END_REF][START_REF] Aki | Attenuation of shear-waves in the lithosphere for frequencies from 0.05 to 25 hz[END_REF]. This principle assumes that, at a fixed lapse time, seismic coda amplitude becomes independent of the source-receiver distance. Therefore, the energy in the coda should be proportional to the radiated energy from the seismic source. However, this method shows some limitations, particularly for media with variations of attenuation [START_REF] Hoshiba | Simulation of coda wave envelope in depth dependent scattering and absorption structure[END_REF] or elastic parameters [START_REF] Hoshiba | Seismic coda wave envelope in depth-dependent s wave velocity structure[END_REF]Yoshimoto, 2000). These studies show that the coda excitation does not only depend on the source size, but also on the source depth and medium parameters. As a result, the estimated radiated energy or magnitude may be biased, which can lead to uncertainties in seismic hazard estimation.

In this paper, we will clarify the role of wave polarization in the complex interplay between scattering by random heterogeneities and deterministic reflections in a simple waveguide geometry. The theoretical aspect of the transport of polarized elastic waves in a crustal waveguide and the Monte-Carlo implementation will be outlined in the next section. Next, we discuss the shape of the modeled energy envelopes, the crustal energy decay and the partitioning energy ratios between the different polarizations. The relative impacts of interfaces, scattering and intrinsic attenuation on the crustal energy decay, the angular distribution of the flux and the energy partition is clarified with the aid of two parametric studies. Finally, the effects of the source radiation on the coda excitation is examined through 2 additional parametric studies, where the coda excitation by explosion or shear dislocations are compared.

2 Radiative transfer in a multiple-scattering crustal waveguide

Model presentation and important scale lengths

In this paper, we introduce a stratified medium consisting of a heterogeneous crust overlying a homogeneous mantle, as a simple but nevertheless realistic model of the lithosphere. The geometry is depicted in Fig. 1a. Crust and mantle are separated by a horizontal velocity contrast (Moho) at a constant depth H where the energy can be either reflected or transmitted. In the numerical examples we set H = 30 km.

The crust is bounded at the top by a free surface where only reflection is possible.

In Fig. 1a, the subscript 1 and 2 designate the elastic parameters of the crust and mantle, respectively. The shear and longitudinal velocities (β, α) and the density (ρ) of the crust and mantle are given by: β 1 =3.5 km/s, α 1 = √ 3β 1 ≈ 6.06 km/s, ρ 1 =2.8, β 2 =4.7 km/s, α 2 = √ 3β 2 ≈ 8.14 km/s and ρ 2 =3.3. In our model, the crust consists of a random medium with fluctuations of velocity and density δα, δβ and δρ, superposed on a homogeneous background. Following Birch's law, the continuous fluctuations of velocity are assumed to be proportional to the density perturbations:

δα α = δβ β = 1 ν δρ ρ . ( 1
)
In this study, the parameter ν of Birch's law is set to 0.8 as suggested for lithospheric rocks [START_REF] Sato | Attenuation and envelope formation of three-component seismograms of small local earthquakes in randomly inhomogeneous lithosphere[END_REF]. To describe the spatial correlation of the fluctuations, we choose an exponential function with scale length l c . This simple hypothesis is supported by borehole observations (e.g. Wu et al., 1994). In a model of radiative transport, the strength of heterogeneities is quantified by the mean free path which represents the average distance between two scattering events. Equivalently, the mean free path l can be defined as the characteristic scale length of attenuation of the coherent wave as a consequence of scattering by heterogeneities. This quantity is a sensitive function of the total variance of the fluctuations 2 , the correlation length l c and the (circular) frequency of the probing waves ω. It may be evaluated approximately with the aid of the Born approximation [START_REF] Sato | Attenuation and envelope formation of three-component seismograms of small local earthquakes in randomly inhomogeneous lithosphere[END_REF] or more rigorously from the Dyson equation [START_REF] Weaver | Diffusivity of ultrasound in polycrystals[END_REF]. In this work, we employed the former approach to calculate all the parameters of the transport model. It is valid below the high-frequency geometrical limit, i.e., the a-dimensional parameter ωl c 2 /c must be smaller than 1 (c stands for α or β) (e.g. [START_REF] Weaver | Diffusivity of ultrasound in polycrystals[END_REF]. Note that in the case of elastic waves, the mean free path also depends on the propagation mode (P or S) and will be denoted by l P and l S .

It is important to point out that the scattering pattern is never isotropic in elastic media. This property has profound effects on the overall transport of energy in the heterogeneous medium and calls for the introduction of another scale length known as the transport mean free path l * . This quantity can be understood as the average distance over which a coherent wave loses "memory" of its initial direction of propagation as a consequence of multiple non-isotropic events. The ratio between l * and l depends critically on the scattering mechanism with l * /l < 1 for preferentially backward scattering and l * /l > 1 for preferentially forward scattering. In the case of elastic waves, the derivation of the transport mean free paths of elastic waves l * P and l * S has been summarized by [START_REF] Turner | Scattering and diffusion of seismic waves[END_REF].

Earth crust is not perfectly elastic and this gives rise to irreversible conversion of seismic energy into other forms. In our model, these complex absorption processes are described at a phenomenological level. We assume that the energy carried by a wave decreases exponentially as it propagates through the medium. The characteristic scale length of attenuation depends again on the propagation mode and will be denoted by l a P and l a S for P -and S-waves, respectively.

As illustrated in Fig. 1a, we suppose for simplicity that there is no scattering in the mantle (hence the infinite value of the mean free paths l P 2,S2 ). Actually, scattering in the mantle is required to explain P -coda observations at teleseismic distances [START_REF] Mancinelli | Constraints on the heterogeneity spectrum of earth's upper mantle[END_REF]. Our model will nevertheless be applicable when the ratio between the transport mean free paths in the crust and the mantle is sufficiently small. More precisely, we require that the return time of waves transmitted to the mantle be large compared to our observation time window. This will be the case if the fluctuations in the mantle are either very weak (to decrease back-scattering) or large scale (to promote forward scattering), or both. The results by [START_REF] Mancinelli | Constraints on the heterogeneity spectrum of earth's upper mantle[END_REF] suggest that the later condition is likely fulfilled. It is worth noting that, due to the geometry of our model, the transmitted energy in the mantle cannot come back to the crust. Therefore, even without absorption the crustal energy can only decrease over time, a phenomenon called energy leakage [START_REF] Korn | A modified energy flux model for lithospheric scattering of teleseismic body waves[END_REF][START_REF] Margerin | Radiative transfer and diffusion of waves in a layered medium: new insight into coda q[END_REF]. As shown in [START_REF] Margerin | Radiative transfer and diffusion of waves in a layered medium: new insight into coda q[END_REF] and [START_REF] Margerin | Breakdown of equipartition in diffuse fields caused by energy leakage[END_REF], one can distinguish two regimes for leakage. In the case of strong scattering in the crust (l * /H 1), the typical leakage time scales like H 2 /(cl * ) [START_REF] Margerin | Residence time of diffuse waves in the crust as a physical interpretation of coda q: application to seismograms recorded in mexico[END_REF][START_REF] Wegler | Diffusion of seismic waves in a thick layer: Theory and application to vesuvius volcano[END_REF]. This formula can be interpreted as the typical time required for diffuse energy to propagate through the crust. In the opposite weak scattering regime (l * /H 1), [START_REF] Margerin | Breakdown of equipartition in diffuse fields caused by energy leakage[END_REF] shows numerically that the leakage time increases like l * /c. Physically, this means that the randomization of propagation directions by multiple scattering controls the leakage out of the waveguide in the case where reflections from the waveguide boundaries play a prominent role. In the following, we will employ the adimensional parameter l * S /H to compare the relative importance of scattering in the bulk of the crust to deterministic reflections at the Moho and the free surface.

Before presenting in details our physical model, we note that an alternative approach to the treatment of scattering in a waveguide has been proposed by Trégourès and van Tiggelen (2002b); [START_REF] Borcea | Onset of energy equipartition among surface and body waves[END_REF]. These authors developed a transport theory based on a surface wave representation of the wavefield, rather than a body wave representation in our case. In the quasi-2D theory of Trégourès and van Tiggelen (2002b), the surface eigenmodes of the waveguide exchange energy as they propagate and interact with the inhomogeneities. A distinctive advantage of this approach is that it incorporates naturally Rayleigh and Love waves. In our model, Rayleigh waves are still lacking although recent developments let us foresee that this could be remedied (Xu et al., 2022). Love waves, which consist of a coherent superposition of multiply-reflected horizontally polarized S-waves (SH), are represented as an incoherent superposition of SH energy in our case. These physical limitations of our model should be noted. Nevertheless, our approach can cover a broad range of cases that are not yet accessible with the quasi-2D theory. Indeed, in the strong scattering regime, i.e. (l < H), the coherent waves are strongly attenuated in between two waveguide boundary reflections, which implies that Love waves cannot develop in this case (Trégourès and van Tiggelen, 2002b). Another important limitation of the quasi 2-D theory is that the lower boundary acts as a perfect reflector, so that energy leakage cannot be modeled yet in this approach. In the next section, we recall the definition of Stokes parameters and briefly discuss their physical interpretation.

Stokes vector and polarization

The central quantity in radiative transfer theory is the Stokes vector which represents both the intensity and the polarization of the waves propagating in a given direction of space (at a given position and time). It should be noted that polarization is entirely determined by the correlations among the components of the wavefield (e.g. [START_REF] Wolf | Introduction to the Theory of Coherence and Polarization of Light[END_REF]. Potentially, there are 9 independent components to fully represent the polarization of a 3-component wavefield. But as shown by [START_REF] Weaver | Diffusivity of ultrasound in polycrystals[END_REF]; [START_REF] Ryzhik | Transport equations for elastic and other waves in random media[END_REF], 4 terms of the full correlation tensor in fact represent cross terms between P -and S-waves. The contribution of these terms may be shown to be negligible thanks to the large wavespeed difference between the two types of waves. This leaves us with 5 independent parameters. In this work, we adopt the standard Stokes description of polarization (4 parameters for S-waves) augmented with an additional parameter for P -waves.

We will first recall the definition of the elastic Stokes parameters for a single plane wave following the approach of [START_REF] Turner | Radiative transfer and multiple scattering of diffuse ultrasound in polycrystalline media[END_REF]. We introduce a Cartesian frame (x, y, z) and, for simplicity, consider waves propagating in the z direction. The field of a plane wave u(r, t) can be expressed as follows:

u(r, t) =        u Sx u Sy u P        =        A Sx (t)e -iω(t-z/β) A Sy (t)e -iω(t-z/β) A P (t)e -iω(t-z/α)        (2)
where A P , A Sx and A Sy stand for the complex displacement amplitudes of Pwaves and the x and y components of S-waves, respectively. The introduction of complex numbers for the amplitudes allows one to incorporate possible phase shifts between the x and y components of shear waves. As we shall see later, such phase shifts may occur upon reflection at medium boundaries. In Eq. ( 2), we have allowed the amplitudes A P , A Sx and A Sy to depend on time. In transport theory, such a time dependence should be slow compared to the central frequency of the waves ω [START_REF] Weaver | Diffusivity of ultrasound in polycrystals[END_REF][START_REF] Ryzhik | Transport equations for elastic and other waves in random media[END_REF]. In other words, the amplitudes can be considered as approximately constant over a large number of periods 2π/ω. In order to allow for a complete description of S-wave polarization at the level of transport theory, a 5-dimensional Stokes vector S of specific intensities is introduced:

S = (I P , I S , Q, U, V ) T (3)
We recall that the term "specific" refers to an angularly resolved flux of energy and, as such, depends on the space direction. The components of the Stokes vector depend on the complex amplitudes introduced in Eq. (2) in the following way:

I P = ραω 2 2 |A P | 2 (t) I S = ρβω 2 2 |A Sx | 2 (t) + |A Sy | 2 (t) Q = ρβω 2 2 |A Sx | 2 (t) -|A Sy | 2 (t) U = ρβω 2 e(A Sx (t)A * Sy (t)) V = ρβω 2 m(A Sx (t)A * Sy (t)) , (4)
where the brackets denote averaging over several periods. In Eq. ( 4) I P and I S are recognized as the total intensities of P -waves and S-waves and Q as the difference of S intensities measured along the x and y axes. U and V are related to the cross-correlation between the two components of the S-wave motion. It can be shown that U is similar to Q in the sense that it measures the difference of S intensity along two axes x and y that can be deduced from x and y by a 45 degrees rotation. Finally

V can be defined as the difference of circular left and circular right S intensities.

This last parameter is non-zero when the phases of the two amplitudes A Sx and A Sy are different. This situation can occur in practice and gives rise to elliptically polarized shear waves. The last four components of S are analogous to the Stokes parameters of electromagnetic waves and we refer to [START_REF] Chandrasekhar | Radiative Transfer[END_REF] for a thorough treatment.

Radiative transfer theory

To model energy propagation in the multiple-scattering crustal waveguide, we use the theory of radiative transfer for polarized elastic waves. As shown by previous authors, [START_REF] Weaver | Diffusivity of ultrasound in polycrystals[END_REF][START_REF] Turner | Radiative transfer and multiple scattering of diffuse ultrasound in polycrystalline media[END_REF][START_REF] Ryzhik | Transport equations for elastic and other waves in random media[END_REF][START_REF] Margerin | Introduction to radiative transfer of seismic waves[END_REF], the full Stokes vector S is necessary to formulate the radiative transfer equation, that we recall hereafter:

C -1 ∂ ∂t + p.∇ S(t, r, p) = -l -1 + l a-1 S(t, r, p) + l -1 4π P (p, p )S(t, r, p )d 2 p + S 0 (t, r, p).
(5)

The left-hand side represents the variation of the Stokes vector followed along a small bundle of rays in a scattering and absorbing medium. The first term on the right-hand side corresponds to the decay of the incident specific intensity due to both scattering and intrinsic absorption. The second term is the gain of specific intensity from a beam originally propagating in direction p and scattered into direction p. Eq.

(5) can be interpreted as local balance of energy for polarized elastic waves [START_REF] Turner | Radiative transfer and multiple scattering of diffuse ultrasound in polycrystalline media[END_REF]. The last term S 0 (t, r, p) represents the source and depends on the focal mechanism. In this work we consider both earthquakes and explosions but limit ourselves to point-like sources (see section 2.6 for more detail). The quantity (α, β, β, β, β) is the diagonal 5 × 5 velocity matrix which contains the mean velocities α and β of the P -and S-waves. The matrices l = diag(l P , l S , l S , l S , l S ) and l a = diag(l a P , l a S , l a S , l a S , l a S ) are the mean free path and absorption length 5 × 5 matrices respectively. We note that l -1 = g, the latter being the diagonal matrix for the integrated scattering coefficients. In the gain term, P (p, p ) is the 5 × 5 Mueller matrix that transforms the Stokes vector of an incident wave into the scattered Stokes vector. This matrix describes all the mode conversions and polarization effects in single scattering.

C = diag
Here, we will recall the most important aspects of the theory, in particular the treatment of polarization at scattering events in the bulk or upon reflection at medium boundaries. To motivate the discussion, examples of snapshots of P -and S-wave energies for an explosion source in the crust are displayed in Fig. 1b and Fig. 1c, at two different lapse time t = 0.75τ S and t = 1.5τ S , where τ S is the S-wave mean free time (set to 20 s here). In this Figure, we can recognize ballistic wavefronts, followed by diffusive waves forming the coda. For seismological application, the P and S energy densities are by far the most important quantities. But transport equations give us access to far more detailed information such as the polarization state or the angular spectrum of the wavefield, which are very helpful to understand the physics of multiple-scattering. Isotropy and depolarization of the wavefield are also central to important applications such as the Green's function reconstruction from coda waves [START_REF] Campillo | Long-range correlations in the diffuse seismic coda[END_REF][START_REF] Emoto | Asymmetry of coda cross-correlation function: Dependence of the epicentre location[END_REF]. We will therefore devote some space to the treatment of energy exchange between different polarizations in the following. The waveguide is bounded by a free surface at the top and the Moho at the bottom (represented by the dashed horizontal lines on the snapshot). The energy density is averaged over cylindrical shells with radius and thickness equal to 5 km.

Scattering and Mueller matrix for volume scattering

With the aid of definition (4), we can now write down the matrix relating incident and outgoing Stokes vectors in the single-scattering process. For clarity, we begin with the so-called scattering matrix. This requires the introduction of a few more notations which are depicted in Fig. 2a. The incident wave propagates along direction p and is scattered into direction p . We introduce the scattering plane spanned by the vectors p and p . Two right-handed reference frames denoted by (r, l, p) and (r , l , p ) are employed to express the incident and scattered Stokes vectors, respectively. The directions r (r ) and l (l ) are perpendicular and parallel to the scattering plane, respectively. Θ is the angle between the directions p and p, also called scattering angle. Φ is the rotation angle between the reference frames (x, y, z) and (r, l, p)

around the z-axis. For simplicity, we imagine a scattering volume V sc at the center of the coordinate system, whose size is much larger than the correlation length l c , yet much smaller than the scattering mean free path. We isolate V sc from the rest of the heterogeneous medium and calculate its response to an incident plane wave. For simplicity, we also assume that the scatterers obey spherical symmetry. Following a method similar to [START_REF] Margerin | Introduction to radiative transfer of seismic waves[END_REF], the scattered Stokes vector S sc can be related to the incident vector S i by the scattering matrix F as follows:

S sc = 1 R 2 F S i , ( 6
)
where R is the distance from the scatterer to the observer. F can be further explicited like so:

F =                   |f P P | 2 α/(2β)|f SP | 2 -α/(2β)|f SP | 2 0 0 β/α|f P S | 2 (|f ⊥⊥ | 2 + |f | 2 )/2 (|f ⊥⊥ | 2 -|f | 2 )/2 0 0 -β/α|f P S | 2 (|f ⊥⊥ | 2 -|f | 2 )/2 (|f ⊥⊥ | 2 + |f | 2 )/2 0 0 0 0 0 e(f ⊥⊥ f * ) -m(f ⊥⊥ f * ) 0 0 0 m(f ⊥⊥ f * ) e(f ⊥⊥ f * )                   , ( 7
)
where f * * designate the scattering amplitudes for the different mode conversions.

More specifically, f ⊥⊥ and f designate the S-to-S scattering amplitudes for polarization respectively perpendicular and parallel to the scattering plane. Because the specific intensities are expressed in frames that are attached to the scattering plane, the scattering amplitudes depend only on the scattering angle when the scatterer is rotationally invariant. Assuming a continuous medium with small fluctuations, the amplitudes f * * , adapted from Sato et al. (2012, p. 142) under the Born approximation, can be found in Appendix A.1. Under the small fluctuation assumption, all the scattering amplitudes are real m(f ⊥⊥ f * ) = 0 . Therefore, the scattering process does not modify the phase of the two components of the S-wave.

The directional dependence of the scattered intensity is quantified by the scattering coefficient g * * (Θ, Φ), which is defined as the ratio between the energy scattered into direction (Θ, Φ) per unit solid angle per unit time to the incident energy per unit area per unit time, normalized by the scattering volume. These scattering coefficients depend on the incident polarization and the scattering amplitudes as follows:

g P P (Θ, Φ) = 4π V sc |f P P | 2 g P S (Θ, Φ) = 4π V sc β α |f P S | 2 g SP (Θ, Φ) = 4π V sc α 2β |f SP | 2 1 - Q i cos 2Φ -U i sin 2Φ I i S g SS (Θ, Φ) = 4π V sc |f ⊥⊥ | 2 + |f | 2 2 + |f ⊥⊥ | 2 -|f | 2 2 Q i cos 2Φ + U i sin 2Φ I i S ( 8
)
where V sc is the elementary representative scattering volume that we have singled out. The scattering coefficients have unit inverse length. Upon integration over all scattering angles, the scattering coefficients provide the attenuations caused by the scattering process. As an example, we write: g P = 1/l P = (4π) -1 4π (g P P (Θ, Φ) + g P S (Θ, Φ)) sin ΘdΘdΦ with l P the P -wave scattering mean free path.

All the quantities we have derived so far suffice to completely describe the singlescattering process. To formulate the multiple scattering process, it is important to fix a global reference frame (x, y, z) (see the geometry in Fig. 2b) and adopt a global convention to decompose the Stokes vector. Following [START_REF] Turner | Radiative transfer and multiple scattering of diffuse ultrasound in polycrystalline media[END_REF],

we introduce a spherical coordinate system (θ, φ, ρ) to address this last point. If we think of z as the normal to the interfaces of the medium (Moho and free surface), the spherical coordinate system coincides in fact with the usual (SV , SH, P ) basis for shear vertical, shear horizontal and longitudinal wave motions. In spherical coordinates, the incident and scattered directions of propagation are denoted by (θ 0 , φ 0 ) and (θ 1 , φ 1 ), respectively. Following the conventions outlined above for the reference frames, the scattered Stokes vector S sc is obtained by taking the product of the Mueller matrix P (θ 1 , φ 1 , θ 0 , φ 0 ) with the incident Stokes vector S i . The Mueller matrix is deduced from the scattering matrix as follows:

P (θ 1 , φ 1 , θ 0 , φ 0 ) = L(-ψ 2 )F (Θ)L(ψ 1 ) (9)
where L is a rotation matrix for the Stokes vector (see appendix A.2). The rotations allow us to first transform the coordinates of the incident Stokes vector to the (r, l, p) system attached to the scattering plane, and then take the scattered Stokes vector in the (r , l , p ) system back to the spherical basis. In Appendix A.3, we express the rotation angles in terms of the incident and scattered propagation directions (θ 0 , φ 0 ) and (θ 1 , φ 1 ) using analytical geometry. [START_REF] Turner | Radiative transfer and multiple scattering of diffuse ultrasound in polycrystalline media[END_REF] employ arguments

from spherical trigonometry to show that

ψ 1 = i 1 + π/2 and ψ 2 = -i 2 -π/2, with
i 1 and i 2 the dihedral angles shown in Figure 2b. i 1 is the angle between the vertical plane containing the incident direction of propagation and the scattering plane. i 2 is the angle between the vertical plane containing the scattered direction of propagation and the scattering plane. The two approaches (analytical vs spherical geometry)

obviously yield the same results (see Appendix A.3). It is important to note that in the Born approximation an elliptically polarized plane wave remains elliptically polarized upon scattering. This property may be readily verified from the analytical form of the Mueller matrix. We give an explicit example of such a verification process in Appendix B where we compute the Mueller matrix for reflection at medium boundaries, which is the topic of the next section.

Mueller matrix for interface reflection

When seismic energy reaches an interface, reflection and transmission result in mode conversions and changes in the polarization. The purpose of this section is to provide a concise derivation of the matrices relating the incident and scattered Stokes vector at a flat interface. First, the amplitude and phase of the reflected or transmitted waves are determined from the reflection and transmission coefficients given by [START_REF] Aki | Quantitative seismology[END_REF], adapted to our global reference frame (x, y, z).

For horizontal interfaces, only P -and SV -waves are coupled. Following Margerin (1998, p.125), the reflected Stokes vector S r at the free surface can be written as:

S r = R Surf S i (10)
where R Surf is the Stokes reflection matrix. This matrix can be written as:

R Surf =                R I Ṕ - P R I ŚV - P /2 R I ŚV - P /2 0 0 R I Ṕ -SV (R I ŚV -SV + 1)/2 (R I ŚV -SV -1)/2 0 0 R I Ṕ -SV (R I ŚV -SV -1)/2 (R I ŚV -SV + 1)/2 0 0 0 0 0 R Ú Ù R V Ù 0 0 0 R Ú V R V V                (11)
where R I * * are the reflection coefficients for the specific intensity, which are defined as products between the energy reflection coefficients R E * * [START_REF] Aki | Quantitative seismology[END_REF] and the ratio between the squared velocities of the incident and the reflected waves. This last factor takes into account the change of geometrical spreading associated with the expansion or contraction of a ray bundle upon reflection. This follows immediately from the application of Snell's law at the interface. Note that in absence of mode conversions, the coefficient R I * * and R E * * are equal. The expressions for the coupling coefficients between the parameters U and V of the incident and scattered waves are particularly simple because the reflection coefficient of SH-waves is unity. We find:

R Ú Ù = R V V = e( Ś SV ) and R Ú V = -R V Ù = m( Ś SV )
, where Ś SV is the amplitude reflection coefficient of the SV -to-SV reflection at the surface.

Likewise, we can develop Stokes reflection and transmission matrices R M oho and T M oho for the reflection or transmission of specific intensity at the Moho. Following a method similar to [START_REF] Turner | Ultrasonic radiative transfer in polycrystalline media: Effects of a fluid-solid interface[END_REF], the reflected and transmitted Stokes vectors S r and S t are written as:

S r = R M oho S i (12) S t = T M oho S i (13)
The reflection and transmission matrices are given by:

R M oho =                R I P - Ṕ R I SV - Ṕ /2 R I SV - Ṕ /2 0 0 R I P -ŚV (R I SV -ŚV + R I SH -ŚH )/2 (R I SV -ŚV -R I SH -ŚH
)/2 0 0

R I P -ŚV (R I SV -ŚV -R I SH -ŚH )/2 (R I SV -ŚV + R I SH -ŚH )/2 0 0 0 0 0 R Ù Ú R V Ú 0 0 0 R Ù V R V V                (14)
and

T M oho =                    T I P - P T I SV - P /2 T I SV - P /2 0 0 T I P -SV (T I SV -SV + T I SH -SH )/2 (T I SV -SV -T I SH -SH
)/2 0 0

T I P -SV (T I SV -SV -T I SH -SH )/2 (T I SV -SV + T I SH -SH )/2 0 0 0 0 0 T Ù Ù T V Ù 0 0 0 T Ù V T V V                    , ( 15
)
respectively. The coupling coefficients between U and The derivation of the main elements of the reflection matrix at the Moho is given in Appendix B. We notice that at a post-critical angles of incidence, a phase shift occurs between the SV and SH-waves both at the Moho and the free surface. This implies in particular that an initially linearly polarized S-waves will acquire an elliptical polarization upon reflection. In turn, the polarization influences the scattering pattern in the bulk of the medium as shown in Eq. ( 8). Our analysis therefore demonstrates that polarization introduces a complex interplay between scattering from interfaces and from volumetric heterogeneities. To our knowledge, this phenomenon has been largely ignored in the seismological literature. In particular, assuming that S-waves are linearly polarized, i.e. V = 0, is only valid in unbounded media. Consequences of polarization will be further explored through numerical simulations. The next section is devoted to the treatment of energy radiation at the source for polarized elastic waves.

V are R Ù Ú = R V V = e( S ŚV S Ś * H ), R Ù V = -R V Ú = m( S ŚV S Ś * H ), T Ù Ù = T V V =

Stokes vector for double couple sources

In this work, we will consider two basic types of sources: explosions (isotropic P -wave sources) and dislocations (double-couple sources). In a Stokes parameters framework, the treatment of explosions is particularly simple. The angular dependence of the Stokes vector at the source S 0 , is given by:

S 0 (θ, φ) = 1 4π                1 0 0 0 0                (16)
The Stokes vector has been normalized in a way such that the total intensity at the source equals 1.

For shear dislocations, we assume that they can be adequately represented by point-like double-couples. The dislocation parameters, which can be set arbitrarily, are φ s (strike), δ (dip) and λ (rake). Following [START_REF] Aki | Quantitative seismology[END_REF], the wave amplitude in the far-field of the source is proportional to the radiation patterns of P -, SV -and SH-waves F P , F SV and F SH which may be written as:

F P (θ, φ) = cos λ sin δ sin 2 θ sin 2(φ -φ s ) -cos λ cos δ sin 2θ cos(φ -φ s )
+ sin λ sin 2δ cos 2 θ -sin 2 θ sin 2 (φ -φ s ) + sin λ cos 2δ sin 2θ sin(φ -φ s )

F SV (θ, φ) = sin λ cos 2δ cos 2θ sin(φ -φ s ) -cos λ cos δ cos 2θ cos(φ -φ s ) + 1 2 cos λ sin δ sin 2θ sin 2(φ -φ s ) -1 2 sin λ sin 2δ sin 2θ 1 + sin 2 (φ -φ s ) F SH (θ, φ) = cos λ cos δ cos θ sin(φ -φ s ) + cos λ sin δ sin θ cos 2(φ -φ s ) + sin λ cos 2δ cos θ cos(φ -φ s ) -1 2 sin λ sin 2δ sin θ sin 2(φ -φ s ) (17)
Since there are no phase shifts between the two components of the S-wave, they are linearly polarized so that the last Stokes parameter V is zero. The Stokes vector in the global frame (x, y, z) can be written as:

S 0 (θ, φ) = 5 8π                   3 2 W P F P (θ, φ) 2 W S F SV (θ, φ) 2 + F SH (θ, φ) 2 W S F SV (θ, φ) 2 -F SH (θ, φ) 2 2W S F SV (θ, φ)F SH (θ, φ) 0                   . ( 18
)
Like for the explosion, the Stokes vector has been normalized in such a way that the energy at the source is 1. W P and W S are the fraction of energy release at the source, which are equal to:

W P = 2 2 + 3γ 5 W S = 3γ 5 2 + 3γ 5 , ( 19
)
where γ is the velocity ratio α/β.

Monte-Carlo implementation

The transport of polarized elastic waves in the crustal waveguide has been simulated with the Monte-Carlo method. In this approach, seismic energy is carried by discrete seismic phonons [START_REF] Shearer | The global short-period wavefield modelled with a monte carlo seismic phonon method[END_REF] who obey the laws of classical ray theory in between two scattering events, where their mode (P , S), polarization state (Stokes vector) and propagation direction change randomly. A major advantage of this approach is its flexibility, allowing the treatment of various types of sources, mode conversions and boundaries. In the numerical examples shown in this work, millions or billions of phonon trajectories have been simulated to obtain reasonably accurate results. The Monte-Carlo scheme, illustrated in Fig. 1a in the case of a multiple-scattering waveguide, is divided into 3 steps:

(1) The initial direction of propagation and mode of a unit energy phonon is randomly selected by interpreting Eqs. ( 16)-( 18) in a probabilistic way. In the case of the double-couple source, the initial direction of propagation is simulated by the rejection method [START_REF] Lux | Monte carlo particle transport methods: neutron and photon calculations[END_REF]. The double-couple diagrams simulated with 500 million phonons are shown in Fig. 3a for all the modes. We draw the attention of the reader to the difference of scale between the P -wave and S-wave radiation patterns, where the radiated S-wave energy is around 23.4 times greater than the radiated P -wave energy. Despite small statistical fluctuations, the simulated double-couple diagrams match well the theoretical formulas (Fig. 3b).

(2) Each phonon performs a "random walk" during which its direction of propagation and polarization are modified at each scattering process that occurs in the bulk of the medium or at an interface. The free path length between two scattering events is determined by an exponential probability law with scale length l P or l S for P and S modes, respectively. If the phonon reaches an interface in between two scattering events, its mode, polarization and direction of propagation are modified according to the laws of reflection/transmission outlined above. The new mode of propagation of the phonon is determined by interpreting probabilistically the energy conservation law at the interface. Subsequently, the propagation direction of the phonon is modified deterministically by application of Snell's law. The Stokes vector is then updated by application of Eq. ( 10) and Eq. ( 12)-( 13) at the free surface and Moho, respectively. Note that upon reflection/transmission, the free path length is adjusted to take into account the difference between the P and S scattering mean free paths in the case of mode conversions or of transmission to a medium with different scattering properties. When scattering occurs in the bulk of the medium, Eqs (8) are interpreted probabilistically to determine the new mode and propagation directions of the scattered phonon. This may be achieved by a decomposition into conditional probabilities as explained in [START_REF] Margerin | Monte carlo simulation of multiple scattering of elastic waves[END_REF]. Once the new state of the particle is determined, the scattered Stokes vector is updated by application of formula Eq. ( 9). Finally, we remark that as the phonon propagates through the medium, its weight decreases exponentially as a consequence of intrinsic attenuation. Again, the scale length of exponential decay (l a P , l a S ) depends on the propagation mode (P , S).

(3) As it propagates through the medium, the complete state of the phonon is monitored at discrete time instants i × dt, where dt is the temporal resolution. At each time step, the construction of the energy envelopes for each polarization state is achieved through a weighted average of the Stokes vectors of the phonons present in the sub-volume of the medium containing the receiver of interest. The resulting envelopes are normalized by the total number of phonons to simulate a unit energy source.

Examples of numerical simulation results

The purpose of this section is to explain, through a typical example, the general approach that we have adopted to analyze our results. A more detailed account of the numerical exploration is given in the following sections. Results of Monte-Carlo simulations in a heterogeneous crustal waveguide for an explosion source are presented in Fig. 4. l * S /H and l * S /l S are set to 3.25 and 2.78, respectively and intrinsic attenuation is neglected. These values correspond to a random medium with l c = 1 km, τ S = 10 s, H = 30 km and a central frequency of 1 Hz. We recall that the nondimensional parameter l * S /H measures the relative importance of reflection from interfaces versus scattering in the bulk of the medium. When this parameter is large, the former process has a very large impact on the propagation as illustrated later. The second parameter l * S /l S measures the level of scattering anisotropy. As this parameter increases, the randomization of propagation directions requires an increasingly large number of scattering events. This parameter controls the degree of angular anisotropy of the energy fluxes in the waveguide at long lapse-time.

Typical energy envelopes are shown in Fig. 4a, for epicentral distances ranging from 50 to 300 km. Energy envelopes in the waveguide exhibit a strong direct P -wave arrival, followed by weak multiples caused by boundary reflections. We note that the development of Lg-like phases is made possible by the guiding effect. At sufficiently large post-critical distances, it is apparent that the guided phases dominate the signal. Because the initial source is an explosion in the middle of the crust, the envelope shape of the trapped S-wave has a diffuse character with no visible ballistic arrival. In Fig. 4b, we observe that total energy in the crustal waveguide decays exponentially after a few mean free times. This exponential decay (readily apparent in linear-log scale in Fig. 4b) is a consequence of the leakage effect. A leakage time τ l can therefore be introduced as follows:

E tot (t) = E 0 exp - t τ l , ( 20
)
where E tot (t) is the total energy in the crustal waveguide and E 0 is a constant. Eq. ( 20) is typically valid after a few mean free times. In Fig. 4b, we also observe that the P -and S-waves energies trapped in the crustal waveguide decay exponentially at the exact same rate after a few mean free times. To understand this fact, we may introduce the intuitive idea that energy leakage plays a role equivalent to absorption: some energy is removed from the system by transmission through the Moho at a rate which is different for P -and S-waves. We may then refer to previous works by [START_REF] Margerin | Effect of absorption on energy partition of elastic waves in the seismic coda[END_REF] who discuss the impact of absorption on the equipartition process. These authors show that even if the dissipation rate of one the modes (P or S) is larger than the other, one still reaches a regime where the total energy of both modes decays at the same rate, in agreement with the result shown in Figure 4b.

Figures 4c and4d represent, respectively, the time evolution of the ratios between the P -and the S-wave energies (E P /E S ) and between the SV -and SH-wave energies (E SV /E SH ) in the whole crust. Interestingly, we observe that these ratios stabilize after 10 to 15 S-wave mean free time. However, the limits differ from the theoretical predictions of equipartition in infinite space, as derived by [START_REF] Weaver | On diffuse waves in solid media[END_REF][START_REF] Weaver | Diffusivity of ultrasound in polycrystals[END_REF]; [START_REF] Ryzhik | Transport equations for elastic and other waves in random media[END_REF]:

lim t→∞ E P E S = β 3 1 2α 3 1 lim t→∞ E SV E SH = 1 . ( 21
)
We stress that these results, valid in infinite space, do not depend on the scattering properties, on the source type or on the amount of energy initially released at the source. In this sense, equipartition ratios can be deemed "universal". Eq. ( 21) is the mathematical translation of the fact that S-waves become depolarized as a result of the multiple scattering process. However, in the example shown in Fig. 4, the E P /E S and E SV /E SH ratios stabilize around 35.3 percent below and 12.4 percent above the full space values (Eq. ( 21)), respectively. These results are in sharp contrast with the equipartition principle. The deviation of the E SV /E SH ratio from 1 means that a residual polarization subsists at long lapse-time in the waveguide, in sharp contrast with the full-space predictions. This result shows that a full description of wave polarization is required to model correctly the energy transport in the crust.

The reduction of the ratio E P /E S can be simply understood as a consequence of the leaky nature of P -waves. In this sense, leakage behaves like absorption at a phenomenological level. The equilibration ratio is simply shifted in favor of the mode which is least absorbed [START_REF] Margerin | Effect of absorption on energy partition of elastic waves in the seismic coda[END_REF]. The increase of the E SV /E SH ratio requires a slightly more sophisticated argument. To understand the origin of the observed shift, the guiding of SV and SH polarizations will be examined in Sections 5 and 6 through a parametric study.

Leakage

A parametric study has been carried out to investigate the effects of the scattering properties of the crust on energy leakage. We examined a wide range of mean free paths, with l S /H varying from 0.073 to 37.3 and l * S /l S ranging from 0.80 (predominantly backward scattering) to 41.95 (strong forward scattering). We examined the validity condition of the Born approximation (ωl c 2 /c < 1) for the set of parameters used in this parametric study and found that only the couple of param- 20)). For every sub-figure, note that the time is normalized by the S-wave mean free time

• S .

Leakage

A parametric study has been carried out to investigate the e ects of the scattering properties of the crust on energy leakage. We examined a wide range of mean free paths, with l S /H varying from 0.072 to 37.3 and l ú S /l S ranging from 0.80 (backward scattering) to 41.95 (strong forward scattering). Leakage times • l are estimated from the fit of the crustal energy decay with Fig. 4b in a time window where the stabilization of the energy ratios is reached. The results are shown in Fig. 5, where the inverse of • l is represented as a function of the parameters l ú S /H and l ú S /l S . We observe that leakage reaches maximum e ciency around l ú S = H in the case of forward scattering, as previously noted by [START_REF] Margerin | Residence time of diffuse waves in the crust as a physical interpretation of coda q: application to seismograms recorded in mexico[END_REF] and (Lacombe, 2001, p. 144). However, in the backward scattering case, the maximum e ciency of leakage seems to be reached when l ú S is slightly lower than H. This maximum marks the limit between the di usive 18 21)). For every sub-figure, note that the time is normalized by the S-wave mean free time τ S .

eters (l *

S /H = 0.073, l * S /l S = 0.80) are outside of the range of applicability. These parameters correspond to a medium with a very strong scattering and a dominantly backward mechanism. Leakage times τ l are estimated from the fit of the crustal en-ergy decay with Fig. 4b in a time window where the stabilization of the energy ratios is reached. The results are shown in Fig. 5, where the inverse of τ l is represented as a function of the parameters l * S /H and l * S /l S . We observe that leakage reaches maximum efficiency around l * S = H in the case of forward scattering, as previously noted by [START_REF] Margerin | Residence time of diffuse waves in the crust as a physical interpretation of coda q: application to seismograms recorded in mexico[END_REF]; [START_REF] Wegler | Diffusion of seismic waves in a thick layer: Theory and application to vesuvius volcano[END_REF] and (Lacombe, 2001, p. 144). However, in the backward scattering case, the maximum efficiency of leakage seems to be reached when l * S is slightly lower than H. This maximum marks the limit between the diffusive regime (l * S < H) where leakage is limited by the transit time of diffuse energy through the crust, and the guided regime (l * S > H) where the seismic energy is multiply-reflected at the medium boundaries and leakage occurs mostly through the scattering of post-critically reflected crustal waves. In the diffusive regime, leakage is independent of the details of the scattering mechanism, in the sense that very different levels of scattering anisotropy yield the same leakage rate provided only that the wave diffusivities are equal. The situation is quite different in the guided regime, where the leakage rate depends clearly on the scattering mechanism. Indeed, we observe in Fig. 5 that the leakage rate depends on both l * S /H and l * S /l S in the case l * S > H. This shows that in regions where scattering is mostly weak like in France (l * S ≈ 250 km at 3 Hz [START_REF] Lacombe | Separation of intrinsic absorption and scattering attenuation from lg coda decay in central france using acoustic radiative transfer theory[END_REF]), the scattering mechanism has a strong influence on the energy leakage and may influence the shape of the coda. In particular, leakage effects contribute to its overall decay rate which may lead to an overestimation of attenuation as determined for example from the coda quality factor Q c [START_REF] Margerin | Residence time of diffuse waves in the crust as a physical interpretation of coda q: application to seismograms recorded in mexico[END_REF][START_REF] Wegler | Diffusion of seismic waves in a thick layer: Theory and application to vesuvius volcano[END_REF]. Note that absorption generally dominates scattering in France so that the former mechanism is mostly responsible for the coda decay, except in tectonically active areas such as the Pyrenees [START_REF] Sens-Schönfelder | Laterally heterogeneous scattering explains lg blockage in the pyrenees[END_REF].

10 -1 10 0 10 1 10 2 l * S /H 10 -3 10 -2 τ -1 l (s -1 ) l * S /

Angular distribution of energy fluxes in the waveguide

For an unbounded medium, the multiple scattering process tends to homogenize the distribution of energy in phase space after a large number of scattering events [START_REF] Ryzhik | Transport equations for elastic and other waves in random media[END_REF]. This implies that all propagation directions become equally probable at long lapse-time. However, as shown by [START_REF] Margerin | Breakdown of equipartition in diffuse fields caused by energy leakage[END_REF] for scalar waves, the angular distribution of specific intensity (i.e. of energy fluxes) never reaches perfect isotropy in the long-lapse time limit in an open waveguide. In [START_REF] Margerin | Breakdown of equipartition in diffuse fields caused by energy leakage[END_REF], it is shown that such a breakdown of strict-sense equipartition is mainly due to the leakage effects. In the elastic case, we expect that the angular dependence of both mode conversions and reflectivity at the Moho and free surface eventually contribute to the anisotropy of the wavefield. Based on the same parametric study as described in Section 4, we now investigate the impact of medium boundaries on the anisotropy of the energy fluxes at long lapse time in the crust. Fig. 6 shows polar plots of the angular distribution of the P , S, SH and SV specific intensities averaged over the whole crustal thickness and all azimuths for a set of representative values of the key parameters l * S /H and l * S /l S . Note that the specific intensities have been normalized by the sum of integrated P and S intensities in the crust. Furthermore, the specific intensities have also been averaged temporally in a late time window, such that all energy ratios have reached equilibrium. As an example, in the case where l * S /l S = 2.78, the time window is between 20τ S and 40τ S . For each sub-figure, the asymptotic angular distributions are shown for a fixed scattering mechanism and different values of the parameter l * S /H. We try to keep the discussion as brief as possible and examine asymptotic regimes for the parameter l * S /H. Fig. 6 clearly shows that when l * S H (or in the diffusive regime), the wavefield is nearly isotropic for all the modes. Although leakage introduces a preferential vertical transport direction across the slab, it leaves little imprint on the specific intensities when the crust is thick compared to the mean free path. Furthermore it may be demonstrated that, while the reflection coefficients are angularly dependent at the free surface, they leave no imprint on an incident equipartitioned wavefield. In other words, if the incident fluxes are isotropic and equilibrated (i.e. E S /E P = 2(α 1 /β 1 ) 3 , E SV = E SH ), so will be the reflected wavefield. This fact is demonstrated in Appendix C. In the thick crust regime, the only imprint left by the presence of interfaces is the exponential decay of the coda, not to be confused with dissipation by anelastic processes.

In the thin crust regime l * S H (or in the guided regime), the situation is far more complex. We can clearly observe in Fig. 6 that the wavefield becomes increasingly anisotropic as the the parameter l * S /H increases. The effect is particularly spectacular for P -waves whose distribution of propagation directions becomes sharply peaked along the horizontal. This property may be traced back to their leaky nature. Indeed, the P -waves that survive in the waveguide at long lapse-time must minimize their interactions with the Moho, thereby favoring transport of energy along the free surface. This result directly contradicts the principle of equipartition.

In the next section, we will show that the development of a strongly anisotropic P wavefield is also associated with a very large shift of the E S /E P ratio in favor of

S-waves.

The SV and SH wavefields show similar behaviors with additional complexities for the former in relation to the strong coupling with P -waves at the free surface.

The most prominent feature visible for both SH-and SV -waves is in Fig. 6 the large jump in the energy fluxes between the pre-and post-critical angles of propagation (the critical S-to-S reflection is at θ ≈ 48.1 degrees in our model). We note that the post-critical SV and SH wavefields are nearly isotropic in the backward and non-preferential scattering (l * S /l S ≈ 1) cases. As the scattering becomes more and more peaked in the forward direction, the angular patterns tend to flatten out along the horizontal direction. It also apparent that the SV flux is slightly larger than the SH flux for large l * S /H. This finding is consistent with our previous examination of Fig. 4d where the SV -to-SH energy ratio was found to be greater than 1 at long lapse-time. We will further comment on this residual polarization effect when we analyze the partitioning ratio. Other noticeable differences between the energy flux distributions of SH and SV are the peaks visible around 30-35 degrees incidence for downgoing SV -waves. This extra intensity finds its origin in the P -SV conversions on the free surface for P -waves propagating at nearly grazing incidence. The sharpness of these peaks clearly attenuates as the scattering anisotropy increases. In the case of strong forward scattering, the vanishing of the peaks of downgoing SV intensity can be traced back to the very strong leakage of the P -wave energy out of the waveguide, leaving almost no P -wave energy to be converted to SV -waves at the free surface. In this case of extreme forward scattering, we observe that the angular distributions of SH and SV fluxes are almost identical. This result will resurface when we discuss the asymptotic energy partitioning ratios. As an additional difference of flux distribution between SV -and SH-waves, we also notice smaller amplitude SV intensity peaks in the upward direction (near θ = 150 degrees). They clearly correspond to the Moho reflection of the extra SV intensity propagating at θ = 30 -35 degrees. These peaks are mostly visible when scattering is preferentially backwards or weakly anisotropic, and disappear in the limit of strong forward scattering.

These results show that adding the Moho breaks the isotropy of the wavefield, and therefore, the equipartition principle can no longer be applied in the weak scattering regime. In Fig. 4c and Fig. 4d, we could observe that the breakdown of equipartition has also an impact on the long lapse time energy rations E P /E S and E SV /E SH . In the next section, we will see how we can link the anisotropy of the flux to the shift of the these ratios.

Energy partition

The results of the parametric study described in Sections 4 and 5 have also been exploited to study the influence of the scattering properties on the long lapse-time S /l S is fixed at 0.98 (non-preferential scattering). (c) l * S /l S is fixed at 2.78 (forward scattering). (d) l * S /l S is fixed at 41.95 (strong forward scattering). The angle θ = 0 corresponds to the downward direction (i.e. propagation from the free surface to the Moho). limits energy partitioning ratios averaged over the crust. The results for E P /E S and E SV /E SH are presented in Sections 6.1 and 6.2, respectively. An additional parametric study, which highlights the effects of intrinsic attenuation on energy partitioning in the coda has been carried out and results are presented in Section 6.3.

P -to-S energy ratios

In Fig. 7, the long lapse-time limit of the energy ratio E P /E S is represented as a function of l * S /H for a series of media with different levels of scattering anisotropy quantified by the ratio l * S /l S . The symbol colors correspond to the long lapse-time limit of the energy ratio E SV /E SH . In the thick crust regime (l * S H), the energy ratios E P /E S and E SV /E SH tend to the infinite space equipartition predictions (see Eq. ( 21)) as expected. The leakage rate is so small that it hardly affects the energy partitioning between the different polarizations. Furthermore, as already alluded to in the previous section, an equipartitioned wavefield is left invariant when it interacts with the free surface (see Appendix C). In the thin crust regime (l * S H), on the one hand, E P /E S strongly decreases as the scattering weakens and eventually tends to zero for very large l * S , due to the leaky nature of P -waves. It is only through the coupling with guided S-waves that some P -wave energy can be trapped at long lapse-time in the crust. As the scattering strength decreases, the coupling becomes weaker and weaker which explains the vanishing of the P -wave energy for thin crusts.

The transition between the thick and thin crust regimes occurs in a range of crustal thickness 1 < l * S /H < 20. Remarkably, the overall sigmoid shape of the curve representing E P /E S as a function of l * S /H depends little on the level of scattering anisotropy, be it forward or backward. As an application, the study of Lacombe S /H ≈ 8 in France. From Fig. 7, we deduce 0.25 E P /E S 0.4 which in turn implies that the average crustal E P /E S ratio is 2.5 to 4 times smaller than the value expected in full space. The far more complicated dependence of E SV /E SH on the scattering properties is examined in the next subsection.

SV -to-SH energy ratios

In Fig. 8, the long lapse-time limit of the energy ratio E SV /E SH is represented as a function of l * S /H for a series of media with different levels of scattering anisotropy, as quantified by the ratio l * S /l S . The results of the simulations for backward and forward scattering have been presented separately in Fig. 8a and Fig. 8b to enhance the readability.

Before going into details, it is worth noting that for all the random media that we have simulated the asymptotic E SV /E SH ratio always exceeds 1. We shall find some counterexamples to this result when we introduce absorption in the next section. In the thick crust regime l * S < H, we observe little dependence of the asymptotic SVto-SH ratio on the scattering mechanism. Deviations from the depolarized state (say by more than 10%) at long lapse-time only occur in the thin crust regime l * /H > 1, on which we will focus the discussion. In this regime, for a fixed l * S /H, the ratio E SV /E SH increases from the backward scattering case to the moderately forward scattering case (l * S /l S ≈ 1.5). From this point on, E SV /E SH starts to decrease. The increase of the E SV /E SH ratio (and the joint decrease of the E P /E S ratio) may be explained as follows. First we remark that in the thin crust regime, P -waves tend to propagate near the horizontal direction, as seen in Section 5. This corresponds to grazing incidence angles at the free surface which gives rise to large conversions to SV -waves (see e.g. [START_REF] Aki | Quantitative seismology[END_REF]. Furthermore, there are almost no pre-critical SV -waves in the waveguide (Section 5), which tends to block conversions from SV -to P -waves at the Moho or the free surface. We therefore suggest that part of the extra SV energy comes from the surface conversion of P -waves which tend to propagate mostly horizontally in the crust. Nevertheless, this explanation is not sufficient because if we assume that the whole P -wave energy was converted to SV energy, the SV -to-SH ratio would not exceed 1.19. This is in contradiction with results from our numerical simulations, which show that the SV -to-SH ratio can be as large as 1.28. By direct examination of the angular dependence of the SV and SH reflection coefficients at the Moho [START_REF] Aki | Quantitative seismology[END_REF], we found that pre-critical SV -waves are globally more guided than pre-critical SH-waves. This is far from evident as part of the SV energy couples to mantle P -waves but it is nevertheless true. Thus, the reflection properties of polarized S-waves at the Moho provide a second mechanism to explain the extra SV energy in the crust. In the case of strong forward scattering (l * S /l S > 5), the SV -to-SH ratio limit reaches a maximum and starts to decrease for increasingly large l * S . This behavior can again be understood in the light of reflection properties at the Moho. Indeed, in the thin crust regime, there is less and less S-wave energy propagating at pre-critical angles of incidence as l * S /H increases, hence no differential behavior between the two modes. This explanation is consistent with the similarity of the angular distribution of SV and SH fluxes in the waveguide reported in Fig. 5.

Role of intrinsic attenuation

An additional parametric study, which highlights the effects of intrinsic attenuation on energy ratios in the coda, has been carried out. In this work, absorption is introduced in a purely phenomenological manner and takes the form of a simple exponential decay in the time domain with characteristic absorption time τ a P,S for Pand S-waves, respectively. In Fig. 9, the long lapse time limits of the energy ratios E P /E S and E SV /E SH are represented as a function of the adimensional parameter τ a S /τ S . The simulation results are shown for 4 different couples of scattering parameters. We fixed the P -wave attenuation time as τ a P = 9/4τ a S , which corresponds to the usual assumption that all dissipation is caused by shearing in a Poisson solid (see e.g. Shearer, 1999, p. 114). The relative contribution of scattering to the total attenuation is often quantified with the aid of a seismic albedo defined as:

B = 1/τ S 1/τ S + 1/τ a S (22)
When scattering dominates intrinsic attenuation (i.e. B > 0.5 or τ S < τ a S ), the shift of both the E P /E S and E SV /E SH energy ratios with respect to the elastic case is rather small. As shown analytically in [START_REF] Margerin | Effect of absorption on energy partition of elastic waves in the seismic coda[END_REF] for an infinite medium and in Trégourès and van Tiggelen (2002a) for a multiple scattering waveguide in the diffusive regime, the P -to-S ratio should be shifted in favor of the least attenuated mode. This effect becomes visible in Fig. 9a for sufficiently small albedos (B < 0.5 or τ S > τ a S ). Interestingly, in this strong absorption regime the P -to-S ratio limit becomes almost independent of the scattering parameters and can become largely higher than the infinite space equipartition limit. There is in fact a competition between the leakage effect which tends to decrease the P -to-S ratio and intrinsic attenuation which acts oppositely. When the latter effect dominates it is in fact possible to reach extreme values of E P /E S as high as 100 times the full-space predictions.

Note however that in practice, the coda will likely return to the noise level before any equipartition is reached if absorption dominates scattering.

As illustrated in Fig. 9b, the SV -to-SH ratio is significantly shifted from the elastic medium value only in the strong absorption regime B < 0.5. Because there is no difference in intrinsic attenuation between the two components of S-wave polarization, the overall effect of absorption is less pronounced for the SV -to-SH than for the P -to-S partition ratio. Nevertheless, the dependence of the former on the albedo is more delicate to understand. In the simpler forward scattering case (l * S /l S = 2.78), we remark an increase of the SV -to-SH ratio as the albedo decreases. In the backward scattering case (l * S /l S = 0.80), however, the SV -to-SH ratio first passes through a minimum before increasing again as B → 0. The increase of the SV -to-SH ratio in the strong absorption limit is easy to understand on the grounds that the dominance of the P -wave energy in the waveguide (as seen in Fig. 9a) promotes the coupling with SV -waves through conversions at the free surface. The decrease of the SV -to-SH ratio in a small range of albedo is rather counterintuitive. To understand this behavior we represent in Fig. 10 the distribution of the fluxes of P -, S-, SV -and SH-waves for different ratios τ a S /τ S in the case of forward (Fig. 10a) and backward scattering (Fig. 10b). It is important to keep in mind that for small albedos, the energy is mostly transported by P -waves propagating horizontally in the waveguide. Hence, we expect the P -to-S coupling induced by scattering and mode conversion at interfaces to be at the origin of the observed complexity. To support this hypothesis, we remark that the azimuthally averaged angular patterns for P -to-SH and P -to-SV scattering generated by a horizontal P wavefield (Fig. 10c-d) can readily be identified in the SH and SV angular flux distributions (Fig. 10a-b). Interestingly, for horizontally propagating P -waves, the scattered SH energy tends to propagate predominantly near the horizontal direction, in contrast with SV -waves which can show large lobe of scattered intensity off the incident horizontal direction. In the case of predominantly backward scattering, the SV energy scattered at large angles leaks out of the waveguide, which induces a reduction of the SV energy flux with respect to the SH flux. In the case of forward scattering, the very large peaks of P -to-SV conversions at the free surface tend to mask this effect. Furthermore, the angular distribution of P -to-SV scattered energy is narrower in the forward scattering than in the backward scattering case, which in turn explains the different behavior of the SV -to-SH partition observed in the two scattering regimes. Again, the extreme absorption cases analyzed in this section are probably of academic relevance only, but they illustrate the complex polarization effects that come into play when deterministic reflections interact with volumetric scattering.

Effects of the source mechanism on the coda excitation

In this section, we will focus on the effects of the source mechanism on the coda excitation. We first recall the basic assumption underlying the coda normalization method, initially developed by [START_REF] Rautian | The use of the coda for determination of the earthquake source spectrum[END_REF] and [START_REF] Aki | Attenuation of shear-waves in the lithosphere for frequencies from 0.05 to 25 hz[END_REF], which has been widely used to separate source, propagation and site effects from coda waves.

For a given station, the method stipulates that the coda decay at long lapse-time is independent of epicentral distance and source radiation [START_REF] Aki | Origin of coda waves: source, attenuation, and scattering effects[END_REF].

Therefore, the coda energy level should be proportional to the source size (and site effect) only. We will examine to what extent this principle remains applicable in the moderate differences around the direct arrivals which are due to the changes in the source radiation patterns. Nevertheless, there are still significant differences in the coda excitation. With the scattering parameters adopted in Fig. 11a and Fig. 11b, the coda excitation by the strike slip fault is, respectively, 28.7% and 21.2% larger than the excitation by the horizontal and the reverse fault.

Fig. 11c and Fig. 11d show that the differences in the coda excitation depend on the scattering strength. Although the differences are weak in thick crust regime, the energy ratios in the coda can differ by as much as 60 percents in the thin crust regime for the strike-slip mechanism compared to the horizontal fault or the reverse fault mechanisms. If we suppose that the radiated energy is proportional to the seismic moment, these results show that the magnitude of earthquakes determined from coda waves, with different focal mechanisms, can be partly biased (by a maximum of 0.13 in M w unit). These differences are related to the dynamics of energy transport at early time, which depends on the source mechanism. For the strike slip fault, we find that the energy is preferentially radiated around the horizontal plane (or at postcritical angle), and is less susceptible to be transmitted into the mantle, compared to the horizontal or reverse faults.

In Fig. 12, we compare the energy envelopes generated by an isotropic source of P -waves and a double-couple (strike, dip and rake are all set to 0°), for two representative scattering strengths (l * S /H = 3.25 in Fig. 12a and l * S /H = 12.98 in Fig. 12b). The two sources release the same amount of energy.

As was to be expected, large differences of amplitude are noticeable for the direct arrivals, as a consequence of the differences in source radiation. These differences in the P and Lg phases can be used to discriminate an explosion source from a sheardislocation source (e.g. [START_REF] Sanborn | Combined effects of deterministic and statistical structure on high-frequency regional seismograms[END_REF]. the long-lapse time coda is much more strongly excited by a double-couple than by an explosion, but the decay rate is independent of the source mechanism. In Fig. 13, we represent the ratio between the total crustal energy generated by a double-couple and an explosion. This gives us access to a quantitative estimation of the differences in the coda excitation level. In the thick crust regime (l * S < H), the coda levels are similar for the two types of sources, implying that the coda normalization method is applicable in this case. The differences between the two types of sources become significant in the more common thin crust regime l * S > H, where the coda excitation can differ by more than one order of magnitude between a double couple and an explosion for l * S > 70H (in favor of the double-couple source). As an application, in France where l * S /H ≈ 8 [START_REF] Lacombe | Separation of intrinsic absorption and scattering attenuation from lg coda decay in central france using acoustic radiative transfer theory[END_REF], we expect the coda excitation to differ by a factor 2-3 between an earthquake and an explosion. More specifically, we found empirically that the ratio R between the energy of a double-couple and an explosion in the coda can be summarized by a power-law of the R ∝ (l * S /H) η with η ≈ 0.7. As a result, in the case of weak scattering, the coda normalization method has to be corrected for the difference in source mechanisms. If this cannot be accomplished, cross-comparison of magnitude estimates between different source types should be avoided. We investigated the mechanism at the origin of the large variation of coda excitation between the explosion and the double couple. We found that the differences observed in the coda are in fact inherited from the dynamics of energy transport at early time, when the crustal energy released by the explosion is mostly composed of P -waves. Since P -waves are leaky, a large fraction of the energy released at the source escapes the crust before scattering has time to convert the energy into partially guided S-waves.

This study suggests that for a given family of focal mechanisms, the details of the fault parameters can affect the excitation of the coda at long lapse-time in the thin crust regime l * S > H. However, the differences of energy level are modest and can cause only minor biases in magnitude estimation based on the coda normaliza- Since in practice the parameter l * S depends on frequency (e.g. Sens-Schönfelder and [START_REF] Sens-Schönfelder | Radiative transfer theory for estimation of the seismic moment[END_REF][START_REF] Przybilla | Estimation of crustal scattering parameters with elastic radiative transfer theory[END_REF][START_REF] Calvet | Lapse-time dependence of coda q: Anisotropic multiple-scattering models and application to the pyrenees[END_REF][START_REF] Eulenfeld | Crustal intrinsic and scattering attenuation of high-frequency shear waves in the contiguous united states[END_REF], we may expect a frequency-dependent bias in the coda normalization method. In the case of double couples with different mechanisms, the effect of l * on the amplitude of the bias remains small. Indeed, we find that the difference in coda excitation increases by at most 10 % when l * S increases by a factor of 2. By contrast, the difference in coda excitation for an explosion and a double-couple increases by a factor 1.66 when l * S increases by a factor 2, thereby introducing a potentially frequency dependent bias in the coda normalization method.

Conclusion

In this study, we evaluated the applicability of basic multiple scattering principles (equipartition, depolarization of shear waves and coda normalization) to seismic waves propagating in a heterogeneous crustal waveguide. To carry out this task, we performed Monte-Carlo simulations of transport of elastic waves taking into account deterministic reflection and transmission of waves at the medium boundaries, as well as a complete description of polarization effects based on a 5-dimensional Stokes vector representation. We analyzed the effects of the scattering properties on the energy decay, the energy partitioning and the angular distribution of energy flux in the crustal waveguide and showed that a full description of polarization is required to model correctly the energy transport. While the energy densities of the different wave modes stabilize at long lapse time in the coda, the equilibration ratios generally do not match the prediction of equipartition in full space, except when the crust is thick compared to the transport mean free path. Like in the scalar case [START_REF] Margerin | Residence time of diffuse waves in the crust as a physical interpretation of coda q: application to seismograms recorded in mexico[END_REF][START_REF] Wegler | Diffusion of seismic waves in a thick layer: Theory and application to vesuvius volcano[END_REF], the decay rate common to all the modes can be quantified by a leakage time τ l and leakage is optimum when the transport mean free path is of the order of the crustal thickness (l * s ≈ H).

The effect of crustal interfaces becomes particularly pronounced when the transport mean free path is greater than the crustal thickness, a common situation for the continental lithosphere [START_REF] Sato | Isotropic scattering coefficient of the solid earth[END_REF]. In this thin crust regime regime l * S > H, we find that the wavefield is never isotropic, even at long lapse-time in the coda. This confirms the breakdown of equipartition put forward by [START_REF] Margerin | Breakdown of equipartition in diffuse fields caused by energy leakage[END_REF] in the scalar case. The anisotropy of the fluxes becomes more and more pronounced as the scattering weakens and the number of reflection at interfaces between two scattering events increases. The energy partitioning is shifted to the detriment of P -waves and in favor of SV -waves (E P /E S < 1 2 (β 1 /α 1 ) 3 and E SV /E SH > 1). The magnitude of the shifts depends on the parameters l * S /H and l * S /l S , quantifying the strength of the scattering and the scattering anisotropy, respectively. In the limit of weak scattering, the P -to-S energy ratio in the coda can decrease to nearly zero, as a consequence of the leaky nature of P -waves. In the same limit, the SV -to-SH energy ratio can increase up to 1.28, meaning that the S-waves are not completely depolarized in the waveguide. This increase is simply explained by the fact that the energy reflection coefficient at the Moho is larger for SV than for SH polarizations, although the former mode can convert to P -waves. An additional mechanism favoring the dominance of SV -waves is the strong P -to-SV conversion at the free surface. These results are unaffected in the case of moderate levels of dissipation in the crust (seismic albedo greater than 0.5). If intrinsic attenuation becomes comparable or larger than crustal scattering, both the P -to-S and SV -to-SH energy ratios are affected, compared to the elastic case. We provide complete explanations for the observed ratios in terms of the scattering mechanism. While of academic interests, these results may not always be applicable to real data due to the presence of noise, which hampers the observation of the asymptotic long lapse-time limit.

We also studied the validity of the assumption underpinning the coda normalization method, namely that the source mechanism leaves no imprint on the coda excitation at long lapse time. Through a first parametric study, we find that there are small-to-moderate differences in coda excitation depending on the focal mechanism of the double couple. Through a second parametric study, we demonstrate that the overall energy level in the coda is higher for double couple sources than for explosions in the thin crust (or weak scattering) regime l * > H of the continental crust. The weaker coda excitation related to explosions is a consequence of the strong transmission of P -wave energy to the mantle at short lapse-time.

In weak scattering areas (l * S H), like in France where l * S ≈ 250 km at 3 Hz [START_REF] Lacombe | Separation of intrinsic absorption and scattering attenuation from lg coda decay in central france using acoustic radiative transfer theory[END_REF], earthquake data are expected to be globally far from the equipartition ratios of the infinite medium. A residual polarization in the coda is likely and the source mechanism could leave its imprint in the energy level of the coda, even at long lapse-time. Hence, while basic multiple-scattering principles like equipartition or coda normalization are important guides to our understanding of scattering amplitudes are expressed as:

f P P (Θ) = k 2 P 4π ν -1 + cos Θ + 2 γ 2 sin 2 Θ -2 + 4 γ 2 sin 2 Θ ξ 2k P 1 -cos Θ 2 f P S (Θ) = -k 2 S 4π sin Θ ν 1 - 2 γ cos Θ - 4 γ cos Θ ξ k P 1 -2γ cos Θ + γ 2 f SP (Θ) = k 2 P 4π sin Θ ν 1 - 2 γ cos Θ - 4 γ cos Θ ξ k P 1 -2γ cos Θ + γ 2 f ⊥⊥ (Θ) = k 2 S 4π [ν(cos Θ -cos 2Θ) -2 cos 2Θ] ξ 2k S 1 -cos Θ 2 f (Θ) = k 2 S 4π [ν(cos Θ -1) + 2 cos Θ] ξ 2k S 1 -cos Θ 2 . ( 23
)
ξ is the Fourier transform of the spatial fluctuations of the random medium. k P and k S are the P -wave and S-wave wavenumbers. γ is the velocity ratio α/β. We recall that ν is the Birch coefficient introduced in Eq. ( 1). We note the following reciprocity relation:

f SP = -β 2 α 2 f P S .

A.2 Scattering coefficients

To derive the scattering coefficients, we start from Eq. ( 6) which relates incident and scattered vectors in bases attached to the scattering plane. Suppose now that the scattering plane is oriented at an arbitrary angle Φ with respect to the fixed basis (x, y, z) with respect to which the incident Stokes vector is decomposed (see Figure 2a). Noting that the unit normal vector to the scattering plane may be written as r = r = p ∧ p/|p ∧ p|, we have:

tan Φ = r • y r • x . ( 24
)
The Stokes rotation matrix L(Φ), which expresses the Stokes vector in a new coordinate system (rotated by an angle Φ around the direction of propagation axis), can be written as [START_REF] Turner | Radiative transfer and multiple scattering of diffuse ultrasound in polycrystalline media[END_REF]:

L(Φ) =               
1 0 0 0 0 0 1 0 0 0 0 0 cos 2Φ sin 2Φ 0 0 0 -sin 2Φ cos 2Φ 0

0 0 0 0 1                . ( 25
)
The components I P , I S and V are invariant under such a rotation. Note that L(Φ) = L(Φ ± π) so that the rotation angle is defined modulo π. The scattered Stokes vector S sc in the (r , l , p ) basis may be expressed as (see Figure 2a):

S sc = F (Θ)L(Φ)S i ( 26
)
where Θ is the scattering angle and Φ is the rotation angle to bring the (x, y, z) axes in coincidence with the (r, l, p) axes. From Eq. ( 26) we can deduce the expressions of the scattering coefficients. We first compute the total scattered intensities:

I sc P = 1 R 2 |f P P | 2 I i P + α 2β |f SP | 2 (I i S -Q i cos 2Φ -U i sin 2Φ) I sc S = 1 R 2 β α |f P S | 2 I i P + |f ⊥⊥ | 2 + |f | 2 2 I i S + |f ⊥⊥ | 2 -|f | 2 2 (Q i cos 2Φ + U i sin 2Φ) (27)
In Eq. ( 27), the brackets indicate an average over the realization of the random medium. Note that for notational simplicity the dependence of the scattering amplitudes on Θ is implicit in Eq. ( 27) and following. Recalling that differential crosssections dσ * * /dΩ are defined as the ratio between the energy scattered in a given direction per unit of solid angle and time, and the incident energy flux, we write: where the surface element is given by dS = R 2 dΩ. In turn, the scattering coefficients g * * are related to the cross-sections by:

g * * = 4π V sc dσ * * dΩ ( 29
)
We then deduce the following expressions of the differential scattering coefficients:

g P P (Θ, Φ) = 4π V sc |f P P | 2 g P S (Θ, Φ) = 4π V sc β α |f P S | 2 g SP (Θ, Φ) = 4π V sc α 2β |f SP | 2 1 - Q i cos 2Φ -U i sin 2Φ I i S g SS (Θ, Φ) = 4π V sc |f ⊥⊥ | 2 + |f | 2 2 + |f ⊥⊥ | 2 -|f | 2 2 Q i cos 2Φ + U i sin 2Φ I i S ( 30
)
The case discussed in [START_REF] Sato | Seismic wave propagation and scattering in the heterogeneous earth[END_REF] corresponds to an incident Stokes vector of the form: (0, 1, 1, 0, 0) (S-waves polarized along the x-axis). It may be readily verified that formulas (30) agree with the one given in [START_REF] Sato | Seismic wave propagation and scattering in the heterogeneous earth[END_REF] in this case.

In the case of an incident circular polarization (for example), the Stokes vector is of the form: (0, 1, 0, 0, ±1). Formulas (30) then predict that there is no azimuthal dependence of the differential scattering coefficients. From a physical point of view, this makes sense as both the medium and the incident wave are invariant by rotation density of rays. In other words, the bundle contains N i S dΩ i S statistically independent rays. Up to the common pre-factor ρ 1 β 1 ω 2 N i S /2, we have

I i S =|A θ | 2 + |A φ | 2 Q i =|A θ | 2 -|A φ | 2 U i =2 e(A θ A * φ ) V i =2 m(A θ A * φ ) (38)
for the incident waves, where the coefficients A represent the typical amplitudes of the rays in the bundle. Thanks to the assumption of statistical independence, the Stokes parameters are proportional to the ray density. We denote by dΩ r S and N r S the solid angle subtended by the reflected ray bundle and the number density of reflected rays, respectively. Up to the common pre-factor ρ 1 β 1 ω 2 N r S /2, the Stokes parameters of the reflected rays may be written as: (39)

I r S =|A θ | 2 | S ŚV | 2 + |A φ | 2 | S ŚH | 2 Q r =|A θ | 2 | S ŚV | 2 -|A φ | 2 | S ŚH | 2
Note that the angular aperture of the ray bundle is conserved for S-to-S reflection, i.e. dΩ i S = dΩ r S . This implies in turn N i S = N r S since we keep track of a fixed set of rays. In Eq. ( 38), we have adopted a similar notation as [START_REF] Aki | Quantitative seismology[END_REF] for the displacement reflection coefficients. For instance S ŚV represents the reflection coefficient for a downgoing SV -wave incident at the Moho reflected to an upgoing SV -wave. Due to a different convention for the definition of the upgoing SV wavefield, our displacement reflection coefficients involving an upgoing SV -waves differ by a sign from those given in [START_REF] Aki | Quantitative seismology[END_REF]. From equations ( 38)-(39), S-waves ( Eq. 36). The reflected P -wavefield of a single plane wave may be written as:

u r P = r r rA θ S Ṕ e -iωt-iω(z cos i-x sin i)/α1 (42)

with sin i = sin j α 1 β 1 . Let us note that if any reflected (transmitted) ray is in the evanescent regime, the associated intensity reflection (transmission) coefficient is set to 0. To illustrate the idea, we note for instance that R I SV -P = 0 for 1/α 1 < p < 1/β 1 (p is the ray parameter) both at the Moho and the free surface. Note that in the same range of ray parameters the reflection coefficients R I P -SV is undefined (since P -waves are not propagating in this regime). With this in mind, we obtain for the intensity of the reflected set of P -wave rays:

I r P = ρ 1 α 1 ω 2 N r P | S Ṕ | 2 |A θ | 2 2 (43)
where S Ṕ is the displacement reflection coefficient at the Moho. Using Eq. ( 38), we obtain:

I r P = α 1 | S Ṕ | 2 (I i S -Q i )N r P 2β 1 N i S ( 44
)
Since we follow a given set of rays, the change in ray density is given by the ratio between the incident and reflected solid angles:

N r P N i S = dΩ i S dΩ r P ( 45
)
where dΩ r P denotes the solid angle subtended by the reflected P rays. Using Snell's law, we find:

dΩ i S dΩ r P = β 2 1 cos i α 2 1 cos j (46)
Reporting ( 46) into (44), we arrive at:

I r P = β 2 1 R E SV - Ṕ (I i S -Q i ) 2α 2 1 (47)
as follows:

I P (n o )|n o • z|d 2 n o = R E Ṕ - P (n i )I P (n i )|n i • z|d 2 n i + R E ŚV - P (n j )I SV (n j )|n j • z|d 2 n j (48)
where d 2 n i,o,j represent elementary solid angles (see Fig. 14) and the R E denote the energy reflection coefficients (see also the main text and Appendix B for further details). Using the notations of [START_REF] Aki | Quantitative seismology[END_REF], they can be expressed in terms of the traditional displacement reflection coefficients as follows:

R E Ṕ - P (n i ) = Ṕ P (n i ) 2 R E ŚV - P (n j ) =1 -Ś SV (n j ) 2 (49)
Note that the displacement reflection coefficients are always real for the propagating waves considered in our model. We may now put the ratio between the reflected and incident energy densities in the form:

I P (n o ) I P (n i ) = Ṕ P (n i ) 2 |n i • z|d 2 n i |n o • z|d 2 n o + I S (n j )(1 -Ś SV (n j ) 2 )|n j • z|d 2 n j I P (n i )|n o • z|d 2 n o (50)
We may now apply Snell's law to simplify this expression. For P -to-P reflection we clearly have: |n i • z|d 2 n i = |n o • z|d 2 n o . For S-to-P conversion we differentiate Snell's law to obtain:

cos idi α 1 = cos jdj β 1 ( 51
)
where j is the incidence angle of SV -waves and i the reflection angle of P -waves.

Combining (Eq. ( 51)) with usual Snell's law, we obtain:

|n j • z|d 2 n j |n o • z|d 2 n o = β 2 1 α 2 1 ( 52
)
Reporting the last result in Eq. ( 50) yields: For an equipartitioned upgoing wavefield, the specific intensities I S (n j ) and I P (n i ) are independent of the incidence angles and their ratio is equal to α 2 1 /β 2 1 . Noting the symmetry relation Ś SV (n j ) 2 = Ṕ P (n i ) 2 [START_REF] Aki | Quantitative seismology[END_REF], we arrive at the desired result:

I P (n o ) I P (n i ) = Ṕ P (n i ) 2 + β 2 1 I S (n j )(1 -Ś SV (n j ) 2 ) α 2 1 I P (n i ) (53)
I P (n o ) = I P (n i ) (54)
Following the same line of reasoning, the result can be extended to SV -waves. In conclusion, if the wavefield incident at the free surface is equipartitioned, so is the reflected wavefield.

Fig. 1 :

 1 Fig. 1: (a) Model geometry and parameters of the Monte-Carlo simulations. The seismic source is placed at a depth of 15 km for all simulations. There is no scattering in the mantle, which explains the infinite mean free paths. (b) and (c) Snapshot of P -and S-wave energy propagation in the lithosphere, at a time of 0.75τ S and 1.5τ S respectively, with τ S = 20 s. The energy is generated by an explosion (red cross).

Fig. 2 :

 2 Fig. 2: (a) Reference frames for incident and scattered waves. The vectors r and l are perpendicular and parallel to the scattering plane, respectively. (b) Angle conventions for the derivation of the Mueller matrix in a global reference frame. The scattering plane is represented by the shaded surface.

  e( S SV S S * H ) and T Ù V = -T V Ù = m( S SV S S * H ), where S ŚV and S SV are the amplitude reflection and transmission coefficients of the SV -to-SV reflection at the Moho and S ŚH and S SH are the amplitude reflection and transmission coefficients of the SH-to-SH reflection at the Moho.

Figure 3 :

 3 Figure 3: (a) Simulated radiation patterns of a double-couple (strike, dip and rake are all equal to 0°). From left to right are represented the energy radiation patterns of the modes P , S (sum of SH and SV modes), SV and SH. (b) Corresponding theoretical radiation patterns of a double-couple.

Fig. 3 :

 3 Fig. 3: (a) Simulated radiation patterns of a double-couple (strike, dip and rake are all equal to 0°). From left to right are represented the energy radiation patterns of the modes P , S (sum of SH and SV modes), SV and SH. (b) Corresponding theoretical radiation patterns of a double-couple.

Figure 4 :

 4 Figure4: Monte-Carlo simulation results for a P -wave explosion in the crustal waveguide described in Fig.1a. l ú S /H and l ú S /l S are set to 3.25 and 2.78, respectively. There is no intrinsic absorption. (a) Simulated energy envelopes for epicentral distances ranging from 50 km to 300 km. (b) Evolution of the total energy in the crustal waveguide. The P -waves and S-waves energy contributions are also shown. (c) Evolution of the P -to-S energy ratio in the crust. (d) Evolution of the SV -to-SH energy ratio in the crust. For (c) and (d), the theoretical limits of the ratios for an infinite multi-scattering space are shown by the horizontal dashed line (see Eq. (20)). For every sub-figure, note that the time is normalized by the S-wave mean free time

Fig. 4 :

 4 Fig. 4: Monte-Carlo simulation results for a P -wave explosion in the crustal waveguide described in Fig. 1a. l * S /H and l * S /l S are set to 3.25 and 2.78, respectively. There is no intrinsic absorption. (a) Simulated energy envelopes for epicentral distances ranging from 50 km to 300 km. (b) Evolution of the total energy in the crustal waveguide. The P -waves and S-waves energy contributions are also shown. (c) Evolution of the P -to-S energy ratio in the crust. (d) Evolution of the SV -to-SH energy ratio in the crust. For (c) and (d), the theoretical limits of the ratios for an infinite multi-scattering space are shown by the horizontal dashed line (see Eq. (21)). For

Fig. 5 :

 5 Fig. 5: Inverse of the leakage time τ -1 l as a function of l * S /H and l * S /l S .

Fig. 6 :

 6 Fig. 6: Polar plots of the angular distribution of the specific intensity in the crustal waveguide for P -, S-, SV -and SH-waves for different values of l * S /H. 4 different scattering regimes have been considered: (a) l * S /l S is fixed at 0.80 (backward scattering). (b) l * S /l S is fixed at 0.98 (non-preferential scattering). (c) l * S /l S is fixed at 2.78 (forward scattering). (d) l * S /l S is fixed at 41.95 (strong forward scattering). The angle θ = 0 corresponds to the downward direction (i.e. propagation from the free

Fig. 7 :

 7 Fig. 7: Long lapse-time limit of the P -to-S ratio as a function of l * S /H for different levels of scattering anisotropy as quantified by l * S /l S . The color scale provides the corresponding limit of the E SV /E SH ratio. The infinite space value of E P /E S is represented by the horizontal dashed line.

Fig. 8 :

 8 Fig. 8: Long lapse-time limit of the SV -to-SH ratio as a function of l * S /H for different levels of anisotropy as quantified by l * S /l S . (a) for different backward scattering regimes. (b) for different forward scattering regimes. The color scale provides the corresponding limit of the E P /E S ratio. The infinite space value of E SV /E SH is represented by the horizontal dashed line.

Fig. 9 :

 9 Fig. 9: (a) Long lapse-time limit of the P -to-S ratio as a function of the normalized absorption time of S-waves τ a S /τ S . (b) Long lapse-time limit of the SV -to-SH ratio as a function of τ a S /τ S . Each symbol corresponds to a fixed l * S /H and l * S /l S . The limit in the infinite multiple scattering and purely elastic space is represented by the horizontal dashed line. Vertical black lines indicate error bars (95 % confidence interval).

FigFig. 11 :

 11 Fig. 11: (a) Comparison between the energy envelopes generated by a horizontal fault (continuous blue lines) and by a strike-slip fault (dashed red line). (b) Comparison between the energy envelopes generated by a reverse fault (continuous blue lines) and by a strike-slip fault (dashed red line). The epicentral distances range from 50 to 300 km. For every simulation in (a) and (b), l * S /H and l * S /l S are fixed to 3.25 and 2.78 (forward scattering), respectively. (c) Coda excitation ratios between a strikeslip and a horizontal fault as a function of l * S /H. (d) Coda excitation ratio between a strike-slip and a reverse fault mechanism as a function of l * S /H. for every simulation in (c) and (d), l * S /l S is fixed to 2.78 (forward scattering).

Fig. 12 :

 12 Fig. 12: Comparisons between the energy envelopes generated by a double-couple source (continuous blue lines) and by an isotropic P -wave source (dashed red line). The epicentral distances range from 50 to 300 km. The scattering parameters are set to: (a) l * S /H = 3.25 and l * S /l S = 2.78. (b) l * S /H = 12.98 and l * S /l S = 2.78.

Fig. 13 :

 13 Fig. 13: Coda excitation ratios between an isotropic P -wave source and a doublecouple as a function of l * S /H. l * S /l S is fixed at 2.78 (forward scattering). The continuous line is the best-fitting algebraic function fitting the last 4 points. The horizontal dashed line indicates a coda excitation ratio equal to 1.

Fig. 14 :

 14 Fig. 14: Energy balance at the free surface. The downgoing P -waves propagating around direction n o have two origins: reflection of P -waves with incidence direction n i , or conversion of SV -waves with incidence direction n j . The incidence angles of P -and S-waves are denoted by i and j, respectively. The elementary solid angles of incident P , incident S and reflected P -waves are denoted by d 2 n i , d 2 n j and d 2 n o , respectively.

Acknowledgements

We would like to thank U. Wegler for his constructive criticism which helped to improve the manuscript. Additional comments by an anonymous reviewer were also appreciated. This study was part of the SIGMA-2 project and funded by EDF (Électricité de France) and CEA (Commissariat à l'énergie atomique et aux énergies alternatives). This work was granted access to the HPC resources of CALMIP supercomputing center under the allocation 2020-p20031.

the coda, they may not always be applicable to real data. Correct estimation of scattering (l * S /H, l * S /l S ) and intrinsic attenuation properties (l a S ) in a sufficiently realistic stratified model of Earth are therefore key to correctly interpret the information encoded in the coda. In the future, this work may be continued by improving on existing methods for the estimation of small to moderate earthquake magnitudes based on the physical modeling of the seismogram envelopes (Sens-Schönfelder and [START_REF] Sens-Schönfelder | Radiative transfer theory for estimation of the seismic moment[END_REF][START_REF] Eken | Moment magnitude estimates for central anatolian earthquakes using coda waves[END_REF].

Appendices

A Basic scattering formulas

A.1 Scattering amplitudes

We consider a continuous medium with small fluctuations of the elastic parameters (see the main text for further details). The scattering amplitudes f * * can be found in Sato et al. (2012, p. 142). Because the specific intensities are expressed in frames that are attached to the scattering plane, the scattering amplitudes depend only on the scattering angle Θ when the scatterer is rotationally invariant. These about the propagation direction. We conclude that the scattering coefficients g SP and g SS depend on the incident polarization.

A.3 Rotation angles for the Mueller matrix

The Mueller matrix relates the incident and scattered Stokes vectors using a global convention to decompose the wavefields. Following previous works [START_REF] Turner | Radiative transfer and multiple scattering of diffuse ultrasound in polycrystalline media[END_REF], we introduce a global reference frame (x, y, z) and a spherical basis (θ, φ, ρ) such that θ, φ φ φ and ρ are in the directions of increasing θ (co-latitude), φ (azimuth) and ρ (distance from origin), respectively. With this convention, the derivation of the Mueller matrix requires the following steps: (1) rotate the incident Stokes vector from the (θ, φ φ φ, ρ) frame to the (r, l, p) frame attached to the scattering plane (note that p = ρ); (2) Apply the scattering matrix F ; (3) rotate the scattered Stokes vector from the (r , l , p ) frame attached to the scattering plane to the (θ , φ φ φ , ρ ) frame (the primes refer to the scattered wave direction with ρ = p ).

The angle ψ 1 that realizes step (1) is defined by:

where the unit normal vector to the scattering plane is given by:

Using the definition of the scattering angle Θ we note:

From Eq. ( 31), we arrive after some algebra at:

The angle ψ 2 which brings (θ , φ φ φ , ρ ) in coincidence with (r , l , p ) is defined by

Note that rotation by an angle -ψ 2 needs to be applied to the scattered Stokes vector (see Eq. 9).

B Construction of the reflection matrix at the Moho

The elements of the reflection matrices (11, 14 and 15) can be derived from the classical treatment of [START_REF] Aki | Quantitative seismology[END_REF] by applying the superposition principle. Let us illustrate the method by deriving some typical Stokes parameters reflection coefficients for incident S-waves on the Moho. The vertical z-axis is oriented downwards. In the spherical basis (θ, φ, r), an arbitrarily polarized incident shear plane wave u i S may be written as:

u i S = θ θ θA θ e -iωt+iω(z cos j+x sin j)/β1 + φ φ φA φ e -iωt+iω(z cos j+x sin j)/β1 (36)

with j the angle between r and z. The θ θ θ, φ φ φ components correspond to the usual SV and SH polarizations, respectively. In general the coefficients A may be complex to take into account possible phase shifts between the two components of the shear wave motion. The reflected S-wavefield u r S may be written as:

with φ φ φ = φ r φ r φ r and θ r θ r θ r is the mirror image of θ θ θ by the interface. Let us now relate incident and reflected Stokes parameters. To do so, we consider an incident bundle of rays of typical angular aperture dΩ i S = sin jdjdφ and denote by N i S the number we deduce the following relations between incident and reflected Stokes parameters:

which agrees with the formulas given in the text. These results are valid for any incidence angle 0 ≤ j < π/2. From Eq. ( 40) we deduce:

This last Eq. shows that a fully polarized incident wave (i.e. such that

When all the transmitted and refracted rays are in the propagation regime, the reflection coefficients for Stokes parameters are real. Furthermore, the off-diagonal reflection coefficients for U and V are equal to zero (because they are related to the imaginary part of the displacement reflection coefficients). As a consequence, an initially linearly polarized S-wave (V i = 0) will remain linearly polarized upon reflection (V r = 0).

When some rays are in the evanescent regime (which occurs when sin j > β 2 /α 1 ) the displacement reflection coefficients of SV -waves are complex and the Stokes parameters of reflected S-waves will exhibit a non-zero parameter V characteristic of elliptic polarizations. This elliptical polarization stems from the fact that the reflected SV and SH components are phase shifted with respect to one another due to their different reflection coefficients at the Moho.

Following the same line of reasoning, mode conversions may be treated with a notable difference related to Snell's law. Consider again an incident bundle of plane

where R E SV -Ṕ = α 1 cos i| S Ṕ | 2 /β 1 cos j is the energy reflection coefficient as defined by [START_REF] Aki | Quantitative seismology[END_REF]. The formula is valid for any incident angle such that sin j < β 1 α 1 . The change of ray divergence at the Moho is at the origin of the factor β 2 1 /α 2 1 in Eq. ( 47). Note that this purely geometrical factor is automatically accounted for in the numerical simulations since the seismic phonons obey Snell's law. In the text, we have lumped together the geometrical factor with the energy

C Reflection of a diffuse field at the free surface

In this Appendix, we demonstrate that the free surface does not affect an equilibrated and isotropic wavefield. We recall that when this condition applies, we have equal amount of SV -and SH-waves incident at the free surface, which acts as a mirror for the latter. As a consequence, the isotropy of SH-waves is obviously preserved.

The critical aspect pertains to the mode coupling between P -and SV -waves which strongly depends on the incidence angles. Let us illustrate the reasoning in the case of a beam of dowgoing P -waves at the free surface (see Fig. 14). This beam originates from the reflection of a set of upgoing P -and SV -wave rays. To properly quantify the contribution of SV -waves, two important physical effects have to be taken into account: mode conversions and focusing/defocusing of the beam of rays. We denote by I P (n i ) and I SV (n j ) the specific intensities of P -and SV -waves incident at the free surface. For the reflected P -waves we use the notation I P (n o ). Note that the unit vectors n i , n j and n o denoting the propagating direction of the rays are related by Snell's law. The principle of energy conservation for P -waves may be expressed

Conflict of interest

The authors declare that they have no conflict of interest.