Shahram Khazaie

Régis Cottereau
email: regis.cottereau@centralesupelec.fr

Didier Clouteau

Numerical observation of the equipartition regime in a 3D random elastic medium, and discussion of the limiting parameters

Keywords: Diffusion regime, Equipartition time, Randomly heterogeneous media, Localization regime

published or not. The documents may come     Numerical observation of the equipartition regime in a 3D random elastic medium, and discussion of the limiting parameters Shahram Khazaie, Régis Cottereau, Didier Clouteau

Introduction

The high-frequency seismograms recorded on the Earth's free surface typically start with a ballistic part (first arrivals of the direct P, S or Rayleigh waves), followed by the so-called coda [START_REF] Aki | Analysis of the seismic coda of local earthquakes as scattered waves[END_REF][START_REF] Herraiz | Coda waves: a review, pure and applied geophysics[END_REF], which has more incoherent and randomized features. Since the pioneering work of Aki [START_REF] Aki | Analysis of the seismic coda of local earthquakes as scattered waves[END_REF], the origin of the coda was identified in scattering on internal heterogeneities (see also [START_REF] Hoshiba | Simulation of multiple-scattered coda wave excitation based on the energy conservation law[END_REF][START_REF] Turner | Scattering and diffusion of seismic waves[END_REF][START_REF] Margerin | Monte carlo simulation of multiple scattering of elastic waves[END_REF][START_REF] Margerin | Introduction to radiative transfer of seismic waves[END_REF]). The analysis of the coda reveals the predominance of multiply scattered waves of shear (S) type [START_REF] Sato | Seismic Wave Propagation and Scattering in the Heterogeneous Earth: Second Edition[END_REF]. The rate at which the energy of the envelope of the coda decreases is shown to be largely independent of the initial source and the local site conditions [START_REF] Aki | Origin of coda waves: Source, attenuation, and scattering effects[END_REF]. Another property of the coda is equipartition [START_REF] Weaver | On diffuse waves in solid media[END_REF][START_REF] Weaver | Diffuse elastic waves at a free surface[END_REF], where the total energy of a diffuse wave field is equally distributed between the energy densities of the P and S waves. This energy stabilization state often arises before the wave field becomes completely isotropic [START_REF] Sánchez-Sesma | Retrieval of the greens function from cross correlation: the canonical elastic problem[END_REF][START_REF] Margerin | Diffusion approximation with polarization and resonance effects for the modelling of seismic waves in strongly scattering small-scale media[END_REF]. Equipartition is experimentally observed at the surface of the Earth (in Mexico for instance [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]), and the ratio of the S-to-P energies is seen to stabilize at a value (about 7) which is different from the expected one for full spaces (2K 3 = 2(v p /v s ) 3 ∼ 10.4 for typical parameter values). This discrepancy is due, on the one hand, to anelastic dissipation, which favors P waves, and on the other hand, to the presence of Rayleigh waves. The influence of anelastic dissipation on the stabilization ratio is discussed in Margerin et al. [START_REF] Margerin | Effect of absorption on energy partition of elastic waves in the seismic coda[END_REF]. The tectonic activity of the considered region is influential on the stabilization values [START_REF] Kubanza | Spatial variation of lithospheric heterogeneity on the globe as revealed from transverse amplitudes of short-period teleseismic p-waves[END_REF], which is out of the scope of this study.

Theoretical asymptotic analysis in the weakly scattering regime (wavelength λ of the source and the order of magnitude of the typical scale of the heterogeneities of the medium c , both much smaller than the propagation length (source-station distance) L, and weak fluctuations of the medium properties) describes in a precise manner how diffusion of the energy arises from the wave equation [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][START_REF] Zeng | Theory of scattered P-wave and S-wave energy in a random isotropic scattering medium[END_REF][START_REF] Fouque | Wave Propagation and Time Reversal in Randomly Layered Media[END_REF], with radiative transfer as an intermediate step [START_REF] Chandrasekhar | Radiative transfer[END_REF][START_REF] Ishimaru | Wave Propagation and Scattering in Random Media[END_REF][START_REF] Sheng | Introduction to Wave Scattering, Localization, and Mesoscopic Phenomena[END_REF]. Diffusion models reproduce the experimental observations of the coda well, including equipartition. The pressure (P) and S waves lose track of the source that created them and the S waves become depolarized irrespective of their initial state of polarization. Although these asymptotic models are mainly valid for full spaces, additional assumptions can provide interesting results in limited cases. For instance, the influence of Rayleigh waves on equipartition can be modeled (see Section 3.3, and also Sánchez-Sesma et al. [START_REF] Sánchez-Sesma | Energy partitions among elastic waves for dynamic surface loads in a semi-infinite solid[END_REF]), and that of topography can be approached through a spectral decomposition of the elastodynamic wave operator for layered half-spaces [START_REF] Margerin | Generalized eigenfunctions of layered elastic media and application to diffuse fields[END_REF][START_REF] Margerin | Energy partition of seismic coda waves in layered media: theory and application to pinyon flats observatory[END_REF]. Imperatori and Mai [START_REF] Imperatori | The role of topography and lateral velocity heterogeneities on near-source scattering and ground-motion variability[END_REF] investigated the influence of the topography and the random heterogeneities on the scattering of elastic waves and showed that both of them produce considerable amount of scattered wave energy with different characteristics.

With the increasing availability of computer tools, stochastic models of the medium properties are being widely used [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] Pitarka | Simulating forward and backward scattering in viscoelastic 3D media with random velocity variations and basin structure[END_REF][START_REF] Hartzell | Effects of 3D random correlated velocity perturbations on predicted ground motions[END_REF][START_REF] Ta | Modeling of random anisotropic elastic media and impact on wave propagation[END_REF][START_REF] Imperatori | Broad-band near-field ground motion simulations in 3-dimensional scattering media[END_REF][START_REF] Takemura | Scattering of high-frequency seismic waves caused by irregular surface topography and small-scale velocity inhomogeneity[END_REF] and diffusion of the energy can be observed directly from large-scale wave equation simulations. Margerin et al. [START_REF] Margerin | Monte carlo simulation of multiple scattering of elastic waves[END_REF] simulated the propagation of full elastic waves in a 3D full space with discrete random spherical heterogeneities. They estimated the time required to reach the equipartition t eq in both weakly and strongly heterogeneous regimes and discussed the influence of the source type, the source-station distance and the ratio a/λ, where a is the diameter of the inclusions. They were however limited to the simple scattering cross section of spheres and did not consider the existence of the free surface. In similar media [START_REF] Margerin | Effect of absorption on energy partition of elastic waves in the seismic coda[END_REF][START_REF] Trégourès | Generalized diffusion equation for multiple scattered elastic waves[END_REF][START_REF] Nakahara | Radiative transfer of elastic waves in two-dimensional isotropic scattering media: Semi-analytical approach for isotropic source radiation[END_REF], equipartition onset time were also numerically evaluated.

The objective of this work is to construct a numerical model to reproduce the onset of the equipartition regime in realistic media. We therefore study the influence of various parameters, both physical (ratio of background Pand S-velocities, propagation length, variance of the heterogeneities) or numerical (influence of Perfectly Matched Layers), on the onset of the equipartition regime. The construction of this numerical model is validated with respect to known theoretical results. In this paper, the propagation medium is considered to be continuously heterogeneous, non-dissipative and isotropic, and the full vectorial nature of the elastic waves is accounted for. Numerical results are not averaged over different realizations of the random medium, and only space averages over a single realization of the randomly heterogeneous medium are considered for the energies. The results are therefore coherent with what is actually observed on real recordings.

The outline of the paper is the following. In Section 2, the wave equation, the characteristics of the propagation medium and their random description along with the scattering cross sections are discussed. We also introduce some properties of a diffusion regime like the equipartition of energies on which we are focusing in this paper. This phenomenon is introduced as an equipartition regime in Section 3 wherein the limit ratio of the body wave energy is discussed. The influence of the presence of a free surface is then discussed. In Section 4, we introduce a numerical simulation in which the equipartition regime is observed. Other simulations in which this regime is not observed along with some explanations about the influencing parameters are presented in Section 5.

Elastic wave propagation and diffusion in a randomly heterogeneous medium

The general concepts related to the elastic waves propagating through three dimensional media of which mechanical parameters are modeled via stochastic approaches will be discussed in this section.

Wave equation and mechanical parameters

We will consider non-dissipative elastic media in which the wave energy losses due to friction, viscosity and radiation are negligible. The wave field u(x, t), being the displacement vector in space and time, is governed by the following wave equation:

ρ ∂ 2 ∂t 2 u(x, t) -∇ x • (C(x) : ∇ x ⊗ u(x, t)) = S (x , t), (x, x , t) ∈ Ω × Ω s × R + (1)
in which Ω s , Ω ⊆ R 3 are respectively the source support and the propagation medium characterized by a constant density ρ and a local fourth-order elastic tensor C(x) = C i jkl (x) that depends continuously on the position x. It should be pointed out that the spatial fluctuations of the density are considered as negligible compared to that of the elastic tensor for the sake of simplicity. Hence, we neglect it even though the former has been shown to be influential on the attenuation of wave energy and therefore on its propagation regime [START_REF] Frankel | Energy-flux model of seismic coda: Separation of scattering and intrinsic attenuation[END_REF][START_REF] Turner | Scattering of elastic waves in heterogeneous media with local isotropy[END_REF]. In this paper, we consider the case of locally isotropic material behavior so that the elastic tensor can be completely described by the local values of the bulk κ(x) and shear µ(x) moduli (or equivalently by the Lamé parameter fields λ(x) and µ(x)):

C i jk (x) = κ(x)δ i j δ k + µ(x) δ ik δ j + δ i δ jk - 2 3 δ i j δ k , (2)
wherein δ i j is the Kronecker's delta (δ i j = 1 if i = j, and δ i j = 0 if i j). The local propagation velocities of the body waves are related to the local values of elastic parameters as v p (x) = (κ(x) + 4µ(x)/3)/ρ and v s (x) = µ(x)/ρ. These phase velocities relate the local frequency and the local wavenumbers via the dispersion relations:

ω(x) = v p (x)|k p | = v s (x)|k s |.
Note that the local frequency remains always constant during the scattering processes.

In the sequel, we will see that one of the most relevant mechanical parameters in this study is the ratio between the spatially-averaged values of the phase velocities, defined as:

K = v p v s 4 3 + κ µ , (3)
in which denotes the spatial averaging. The measured values of K for the upper mantle vary between 1.65 and 1.8 [START_REF] Sato | Attenuation and envelope formation of three-component seismograms of small local earthquakes in randomly inhomogeneous lithosphere[END_REF]. The values between 1.7 and 2.3 are observed in laboratory tests under pertinent pressure and temperature conditions [START_REF] Sato | Attenuation and envelope formation of three-component seismograms of small local earthquakes in randomly inhomogeneous lithosphere[END_REF]. The Poisson materials with κ = 5 µ /3 (or equivalently λ = µ ) correspond to a characteristic value of K = √ 3, are frequently used in the geophysics literature [START_REF] Sato | Seismic Wave Propagation and Scattering in the Heterogeneous Earth: Second Edition[END_REF]. The positivity condition of the bulk and shear moduli imposes the minimum theoretically acceptable value of K min = √ 4/3 ≈ 1.16 which is also stated by the thermodynamical theory.

Random description of the propagation medium and scattering cross sections

A probabilistic approach is employed to model the mechanical parameters κ(x) and µ(x). Statistical description of these random fields is elaborated in what follows. In this study, the bulk and shear moduli (eigenvalues of the elastic tensor) are considered as two independent random fields following Gamma distributions (see Guilleminot and Soize [START_REF] Guilleminot | On the statistical dependence for the components of random elasticity tensors exhibiting material symmetry properties[END_REF] for more details). We model the mechanical parameters (κ(x), µ(x)) as second-order statistically homogeneous random fields with specified mean values κ and µ , autocorrelation functions

(ACF) R κ (x-x ) = E[(κ(x)-κ )(κ(x )- κ )] and R µ (x-x ) = E[(µ(x)-µ )(µ(x )-µ )], and cross-correlation function R κµ (x-x ) = E[(κ(x)-κ )(µ(x )-µ )].
We assume that the structure of the spatial correlation functions is isotropic which means that they depend only on the distance between the points r = |xx | and hence can be considered as 1D functions. Another characteristic of a random field is its correlation length which can be imagined as the typical size of the heterogeneities and is defined in this presentation as:

c = 2 ν 2 R + R(r) dr, (4)
in which ν 2 = R(0) is the variance of the corresponding random field. Finally, we introduce the normalized power spectral density function (PSDF) Φ(|k| c ) as the Fourier transform of R(r/ c )/ν 2 . Note that the wavenumber |k| is the Fourier relative of the space variable r.

In the case where the elastic waves propagate through the random medium described above, and in a weakly heterogeneous scattering regime ( c ∼ λ, L λ and ν 2 1) a transport regime occurs in which the wave energy densities can be analytically described. In this regime, the radiative transfer equations (RTE) of elastic waves describe the spatio-temporal evolution of the ensemble-averaged energy densities of the P and S waves. These equations are analogous to the energy conservation equations of the body wave energies. They are derived in the weakly heterogeneous scattering regime for random materials with local isotropy using a multiscale expansion of the Wigner transform of the wave field [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 these equations, the complex interactions between the heterogeneities and the waves are described by the differential and total scattering cross sections σ IJ (k, k ) and Σ IJ (k) (I, J ∈ {P, S }). The differential scattering cross sections σ IJ (k, k ) (with unit m 3 /s) specify the rate at which incident energy of mode type J with wave vector k is scattered to the wave energy of mode type I with wave vector k. The total scattering cross sections Σ IJ (k) (with unit 1/s)

Σ IJ (k) = R 3 σ IJ (k, k )dk , I, J ∈ {P, S}, (5)
represent the rate at which incident energy of mode type J with wave vector k is scattered to the wave energy of mode type I in any direction.

Since the main focus of this paper is on the equipartition regime which is narrowly related to the interchanges between the P and S wave energies during the scattering processes, here we will only introduce the P-S and S-P total scattering cross sections:

Σ PS (k)I 2 = 2K 3 Σ SP (k/K), (6)
in which I 2 is a 2 × 2 identity matrix. It is worth noticing that the total scattering cross sections corresponding to the incident S wave (Σ PS , Σ SS ) and to the incident P wave (Σ PP , Σ SP ) are of dimension 2 and 1, respectively. All these total scattering cross sections are functions of the wavenumber |k| and their definition is known explicitly based on second-order statistics of the wave parameters, [START_REF] Khazaie | Influence of the spatial correlation structure of an elastic random medium on its scattering properties[END_REF]. For instance:

Σ PS (|k|) = 4π 2 √ 1 + K 2 ν 2 µ v p 3 c |k| 4 +1 -1 χ 2 (1 -χ 2 )Φ 1 + K 2 -2Kχ |k| c dχ. ( 7
)
wherein χ = cos θ in which θ is the scattering angle (angle between the incident and scattered waves). We will use the values of Σ PS (|k|) to calculate the global equipartition time in a diffusion regime as we will see later in this paper.

In the next section we study the diffusion regime as the asymptotic limit of the transport equations at long lapse times and some fundamental characteristics of this regime will be mentioned.

Characteristics of the diffusion regime

Over long propagation paths and after many scattering events, the waves lose their memory about their source or initial conditions. The energy evolution of the multiply scattered waves can be described by a diffusion equation that can be derived through an asymptotic analysis of the RTE [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 a diffusion regime we have

∂φ ∂t = ∇ x • (D(|k|)∇ x φ) , (8)
in which φ is the ensemble-averaged energy density of the P or S waves. In this regime, due to the depolarization of S waves, their energy density (described by four Stokes parameters as in the case of electromagnetic waves) becomes proportional to the identity matrix and thus can be written as φI 2 . D is the elastic diffusivity which determines the degree of scattering, i.e. lower values of D imply highly scattering medium and vice versa. In the transport regime, the average distance (resp. time) before the propagation direction of each mode is significantly altered is defined as the scattering mean free path sc m , m ∈ {P, S } (resp. time t sc m , m ∈ {P, S }) of that mode. Several mean free path should be traveled by each wave mode in order for that mode to reach a diffusion state. The typical distances (resp. times) beyond which a diffusion regime sets in for each propagation mode is defined as the transport mean free path * m , m ∈ {P, S } (resp. time t * m , m ∈ {P, S }) of that mode. They are functions of the statistics of the underlying propagation medium and their analytical formulas can be found in [START_REF] Turner | Scattering and diffusion of seismic waves[END_REF][START_REF] Margerin | Attenuation, transport and diffusion of scalar waves in textured random media[END_REF][START_REF] Khazaie | Influence of the spatial correlation structure of an elastic random medium on its scattering properties[END_REF] for instance. The transport mean free times t * P and t * S are related to * P and * S via the corresponding mean phase velocities. In a diffuse wave field, the total energy is equipartitioned between body waves and their ratio follows an equipartition law which will be elaborated in the next section.

Equipartition of elastic waves

The energy stabilization state of a diffuse wave field, which will be hereinafter referred to as an equipartition regime, means that the total energy is equally partitioned between the energy densities of the body waves regardless of the initial conditions (location and magnitude of the source) and the details of the scattering processes. If we denote the spatial averages of the body waves in each point x and at any time t by E p (x, t) and E s (x, t) , in a full elastic medium the so-called equipartitioning law reads ([9, 17, 42, 18, 4])

R = E s (x, t) E p (x, t) = 2K 3 = 2 v p v s 3 2 4 3 + κ µ 3 2 . ( 9
)
This equation is true at all points of a boundaryless elastic medium for times larger than max(t * P , t * S ). As a result, on the one hand, each of the energies E p and E s are individually related to the statistical properties of the medium through the scalar diffusion equation and on the other hand, the ratio R = E s / E p is directly related to the first-order statistics of the underlying randomly heterogeneous medium. The equipartition law was first derived by Weaver [START_REF] Weaver | On diffuse waves in solid media[END_REF] by counting the number of normal wave modes in the context of room acoustics. Then, Weaver [START_REF] Weaver | Diffusivity of ultrasound in polycrystals[END_REF] and Papanicolaou et al. [START_REF] Papanicolaou | Stability of the P to S energy ratio in the diffusive regime[END_REF] derived it rigorously based on the transport theory of elastic waves. A comparison between equations ( 6) and [START_REF] Weaver | On diffuse waves in solid media[END_REF] shows that the stabilization ratio is exactly the value of the ratio between P-to-S and S-to-P total (and forward) scattering cross sections, i.e. R = Σ PS /Σ SP .

When an equipartition regime sets in, a stabilization ratio of about R 10.4 is expected for Poissonian crustal materials (K = √

3). However, Shapiro et al. [START_REF] Shapiro | The energy partitioning and the diffusive character of the seismic coda[END_REF] reported a stabilization value of 7 from local earthquakes in Mexico. This discrepancy from the expected value can be related on the one hand to the effects of the anelastic dissipation being in favor of the P wave energy and on the other hand to the role of the free surface which leads to the dominance of the surface Rayleigh waves. A discussion about the influence of the anelastic dissipation is out of the scope of this paper. However, the influence of the presence of a free surface on the stabilization values will be discussed in Section 3.3.

Global equipartition time

Rather than considering different realizations of a random medium, we will consider only a single realization of a large enough medium (see [START_REF] Bal | Accuracy of transport models for waves in random media[END_REF] for more discussion about the statistical stability of the wave energy density in weakly heterogeneous scattering regime). We then approximate all ensemble averages by the corresponding spatial ones. Hence, the equipartition of the energies in the whole propagation medium and/or over a surface is investigated in this presentation. In this regard, the required time to reach an equipartition regime, denoted by t eq , is called the global equipartition time. In discrete random media with spherically symmetric scatterers, Trégourès and van Tiggelen [START_REF] Trégourès | Generalized diffusion equation for multiple scattered elastic waves[END_REF] derived a sophisticated form of the diffusion equation as a coupled system of equations for the P and S wave energy densities. The solution of this equation in an unbounded medium shows an exponential decay of the body wave energy densities toward the steady-state solution with the same rate. This decay rate is called the global equipartition (stabilization) time and is defined as the time needed to reach an equipartition regime over a slice or the whole medium. It can be defined as

t eq (|k|) = 1 Σ PS (|k|) + Σ SP (|k|/K) = 2K 3 |k| -4 4π 2 √ 1 + K 2 (1 + 2K 3 )ν 2 µ v p 3 c +1 -1 χ 2 (1 -χ 2 )Φ 1 + K 2 -2Kχ |k| c dχ . ( 10
)
Figure 1 depicts the values of the normalized global equipartition time t eq (|k|) = t eq (|k|) v p / c for different correlation models (as elaborated in [START_REF] Khazaie | Influence of the spatial correlation structure of an elastic random medium on its scattering properties[END_REF]) and for two different values of K which will be subsequently used in our numerical simulations. Hence, for instance in high frequencies (i.e. |k| c 1) the influence of the correlation model on the values of t eq is substantial. Furthermore, one can observe that lower values of K result in faster global equipartitioning of the energies. For instance, in high-frequencies for the exponential ACF, increasing the value of K from 1.16 to √ 3 will increase the global equipartitioning time by a factor of 10. The same factor for the case where c ∼ λ (weakly heterogeneous scattering regime) is 5. The influence of the parameter K will be investigated later on via numerical simulations in this paper.

Equipartitioning without free surface

The calculation of the energy densities of wave modes requires the knowledge about the local values of the medium parameters. The objective of this section is therefore to express the equipartitioning law only in terms of the recorded wave field in the case of an open (unbounded) medium in which the surface waves do not appear.

The total energy density of elastic waves E(x, t) is the sum of the kinetic and strain energies. The latter can be decomposed into the P-wave energy E p (x, t), the S-wave energy E s (x, t) and a cross energy. Hence, E(x, t) can be expressed in terms of the wave field u(x, t) and the medium parameters (κ(x), µ(x)) (see [START_REF] Shapiro | The energy partitioning and the diffusive character of the seismic coda[END_REF] or [START_REF] Hennino | Observation of equipartition of seismic waves[END_REF] for instance): 3 (b) and ν κ = ν µ = 0.1: exponential (thick solid line), power-law (thin dashed-dotted line), Gaussian (thick dashed line), triangular (thin solid line) and low-pass white noise (thin dashed line). For more details on the definition of these spatial correlation functions see [START_REF] Khazaie | Influence of the spatial correlation structure of an elastic random medium on its scattering properties[END_REF].

E(x, t) = E k (x, t) + E p (x, t) + E s (x, t) + µ(x) H(x, t) = 1 2 ρ ∂u(x, t) ∂t 2 + 3κ(x) + 4µ(x) 6 (div u(x, t)) 2 + µ(x) 2 ||curlu(x, t)|| 2 + µ(x) H(x, t), (11)
where E k (x, t) is the kinetic energy density and the cross term H reads:

H(x, t) = 2
Let us now consider the coda part of a wave field u(x, t) which is assumed to be diffuse, meaning that its energy is governed by a diffusion equation. In this case, the wave field can be assumed to be a superposition of a priori uncorrelated plane P and S waves coming from different directions [START_REF] Shapiro | The energy partitioning and the diffusive character of the seismic coda[END_REF]. Consequently, the total energy densities of P and S waves could be written as the sum of their individual components. Since the wave field is diffuse, the displacement components will be uncorrelated. As a result, the spatial average of the cross term H will cancel so that the average of the total energy density will be:

E = E k + E p + E s = ρ 2 ∂u(x, t) ∂t 2 + 3κ(x) + 4µ(x) 6 (div u(x)) 2 + µ(x) 2 ||curl u(x)|| 2 . ( 13
)
The assessment of the stabilization ratio R can now be done using a grid of closely spaced sensors which enables us to calculate the derivatives of the displacement vector following all three spatial coordinates (divergence and curl of u in equation ( 13)):

R = µ(x) 2 ||curl u(x, t)|| 2 3κ(x)+4µ(x) 6 (div u(x, t)) 2 = 3 µ 3 κ + 4 µ • ||curl u(x, t)|| 2 (div u(x, t)) 2 = 1 K 2 • ||curl u(x, t)|| 2 (div u(x, t)) 2 , ( 14
)
in which t > max(t * P , t * S ) and we assumed that each of the pairs (µ(x), ||curl u(x, t)|| 2 ) and (3κ(x) + 4µ(x), (div u(x, t)) 2 ) are independent random variables so that the average of their product is the product of the respective averages. The authors verified this assumption through the numerical simulations presented later on and confirm its validity. Since in an open medium at long lapse times we have R = 2K 3 , the equipartitioning law in terms of the spatial derivatives of the recorded wave field for t > max(t * P , t * S ) reads:

||curl u(x, t)|| 2 (div u(x, t)) 2 = 2K 5 . ( 15
)
At the end of this section we point out that in an equipartition regime on an open Poissonian medium (full-space with K = √ 3), the energies are also equally distributed between E k and E p + E s so that the values of the [START_REF] Hennino | Observation of equipartition of seismic waves[END_REF]. In sections 4 and 5 this stabilization is investigated via numerical simulations.

Equipartitioning with a free surface

In practice, the equipartition regime can be observed mostly using seismic sensors on the Earth's surface. The difference with respect to the full space lies in the existence of surface waves and the complex interactions between them and the body waves. In this case, the calculation of the derivative of the wave field in z direction is rather difficult because it requires to install the seismic receivers in depth. Hence, the free surface boundary condition will be used as a tool to simplify the calculation of the stabilization ratio R over the surface of the Earth. The components of the stress tensor in z direction should vanish i.e. σ xz = σ yz = σ zz = 0. These equations result in three equations relating the vertical derivatives of the wave field to the horizontal ones. Using these relations and considering the fact that over the free surface the wave field is a Rayleigh plane wave with u y (z = 0) = 0, u x (z = 0) =sin(ωtkx) and u z (z = 0) = r cos(ωtkx) in which r is the vertical-to-horizontal axis ratio on the free surface, yields the following equipartition law (for t > max(t * P , t * S )) over the free surface (see [START_REF] Shapiro | The energy partitioning and the diffusive character of the seismic coda[END_REF] for more details):

R surface = E s surface E p surface = 1 4 K 2 4 ∂u z (x, t) ∂x 2 surface 4 ∂u x (x, t) ∂x 2 surface = K 2 r 2 = K 2 K 4 r 1 - K 2 r K 2 2 -K 2 r -2 1 -K 2 r 1 -K 2 r K 2 2 , (16)
in which K r = v r /v s is the ratio between Rayleigh and S wave speeds and is the solution of the characteristic equation of Rayleigh waves and surface specifies an averaging over the free surface. This equation for the case of Rayleigh waves propagating through piece-wise homogeneous medium is (see [START_REF] Auld | Acoustic fields and waves in solids[END_REF] for example) K 6 r -8K 4 r + (24 -16K -2 )K 2 r + (16K -2 -16) = 0. As an example, for Poissonian materials (K 2 = 3), the characteristic Rayleigh equation has three roots for K 2 r , among which the acceptable one is 2 -2/ √ 3. This value results in a propagation velocity of Rayleigh waves v r 0.92v s and inserting in equation ( 16) gives a stabilization ratio of R 6.5 that is consistent with 7 which is reported by [START_REF] Shapiro | The energy partitioning and the diffusive character of the seismic coda[END_REF]. Table 1 shows the analytical values for the stabilization of the S to P energy ratio in an unbounded medium as well as over the free surface (assumed to be piece-wise homogeneous) for different values of K. Finally, equation ( 15) over the free surface becomes:

||curl u(x)|| 2 surface (div u(x)) 2 surface = K 2 R surface = K 4 K 4 r 1 - K 2 r K 2 2 -K 2 r -2 1 -K 2 r 1 -K 2 r K 2 2 . ( 17
)
4. Observation of the equipartition regime in a 3D random elastic medium

The numerical simulations of elastic wave propagation in 3D randomly heterogeneous elastic media in this work are carried out via the spectral element-based code, SPEC3D. This parallel code is developed with initial works of [START_REF] Festa | The newmark scheme as velocity-stress time-staggering: an efficient pml implementation for spectral element simulations of elastodynamics[END_REF][START_REF] Delavaud | Simulation numérique de la propagation d'ondes en milieu géologique complexe: application à l'évaluation de la réponse sismique du bassin de Caracas (Venezuela)[END_REF][START_REF] Cupillard | Simulation par la méthode des éléments spectraux des formes d'onde obtenues par corrélation de bruit sismique[END_REF] by the seismology group of Institut de Physique du Globe de Paris1 [START_REF] Cupillard | RegSEM: a versatile code based on the spectral element method to compute seismic wave propagation at the regional scale[END_REF] and has been modified in order to take into account the random elasticity matrix [START_REF] Cottereau | Validation of software for 3D propagation of waves in heterogeneous and random media[END_REF][START_REF] Ta | Modeling of random anisotropic elastic media and impact on wave propagation[END_REF]. A global mesh is generated via a simple MATLAB script and the mesh partitioning is carried out via the software package METIS. We use the spectral representation method ( [START_REF] Shinozuka | Response variability of stochastic finite element systems[END_REF][START_REF] Shinozuka | Simulation of stochastic processes by spectral representation[END_REF]) to simulate Gaussian stochastic germs and then a so-called inverse transform sampling technique [START_REF] Devroye | Non-Uniform Random Variate Generation[END_REF] is employed to generate the stochastic germs following a Gamma probability density function to generate random values of the elastic moduli (κ, µ) (see [START_REF] Khazaie | Influence of the statistical parameters of a random heterogeneous medium on elastic wave scattering : theoretical and numerical approaches[END_REF] for a detailed discussion). We further assume that the coefficients of variation of κ and µ are equal, i.e. δ = ν κ / κ = ν µ / µ . This parameter δ is hereinafter called the dispersion level of the propagation medium and should be less than 70% in order that the random elasticity matrix along with its inverse be second-order, [START_REF] Soize | A nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF][START_REF] Soize | Maximum entropy approach for modeling random uncertainties in transient elastodynamics[END_REF].

The objective of this section is to investigate the equipartitioning phenomenon via the numerical simulations carried out using SPEC3D. The numerical example in this section consists of a propagation medium as a cube defined by Ω = {x ∈ R 3 | -1500 m ≤ x, y ≤ 1500 m; -3000 m ≤ z ≤ 0 m} with continuous random heterogeneities in elastic moduli and a constant density of 2000kg/m 3 . The average value of µ is assumed to be µ = 2 × 10 9 Pa, and the standard deviation of µ is parameterized as ν µ = µ δ = 2 × 10 9 δ. Hence, for a given value of K = v p / v s , the average value of κ and the average phase velocities will be:

κ = K 2 - 4 3 µ ; v s = µ /ρ = 1000 m/s ; v p = K v s . (18)
A half-space Ω with a dispersion level of δ = 0.40 is subjected during 5 s to an explosion source at its central point (x, y, z) = (0, 0, -1500)(m). The source is characterized by a Ricker pulse whose delay time t 0 and central frequency f 0 are 0.3s and 10Hz, respectively. Hence, the dominant wavelengths of the P and S waves will be λ p = 100K (m) and λ s = 100 (m), respectively. In the framework of the SEM, in order to accurately estimate the wave field, at least five control points per minimum wavelength should be used [START_REF] Faccioli | 2D and 3D elastic wave propagation by a pseudo-spectral domain decomposition method[END_REF][START_REF] Fichtner | Full Seismic Waveform Modelling and Inversion[END_REF][START_REF] Cristini | Some illustrative examples of the use of a spectral-element method in ocean acoustics[END_REF]. Hence, the propagation medium is divided by 60 × 60 × 60 elements (the element size in each direction being 50 m) each of which contains 8 × 8 × 8 Gauss-Lobatto-Legendre (GLL) points. As a result, the total number of points of the mesh is (7 × 60 + 1) 3 = 74618461. It should be noted that the shape functions are the Lagrange polynomials of order 7. The time integrations in the framework of the SEM are solved via an explicit Newmark scheme. In addition, the Courant number is 0.27 corresponding to a time discretization of dt = 0.0005 s. Using 96 processors, about 4 hours of run-time is required to compute the wave field during 5 seconds of wave propagation.

The spatial correlation function of the random elastic moduli is considered to be either low-pass white noise or exponential whose normalized PSDFs are respectively Φ(k) = 2 9π 4 H( 3π 2k) (H being the heaviside function) and

Φ(k) = 1 8π 2 (1 + k 2 4 ) -2 .
The value of K is considered to be its minimum theoretically acceptable, i.e. 1.16. The correlation length is c = 100m which is equal to the dominant shear wavelength to maximize the efficiency of the interactions between the waves and the heterogeneities in order to have a highly scattering medium.

Figure 2 shows the temporal variation of the space-averaged body wave energies E Ω p (t), E Ω s (t) and the kinetic energy E Ω k (t), calculated via the following equations: From now on, all figures showing the temporal evolution of the whole space-averaged energies include a shaded window indicating the time interval during which the decay of the energies E Ω p and E Ω s begins. The starting and the end points of this window specify respectively the lapse times corresponding to the first arrivals of the direct P and S waves (the waves traveling in the homogeneous background) to the boundaries of the medium. The bottom left plot of the Figure 2 shows that the energy ratio E Ω s /E Ω p converges to 2.9 and 3.1 respectively for the low-pass white noise and exponential correlation models. Since the integrals in Equation 19 are done over the whole domain including its boundaries, the ratio E Ω s /E Ω p accounts for the region near the free surface (where the stabilization value is R surface < R) along with the region far from the free surface (where the stabilization value is R). As a result, the ratio E Ω s /E Ω p should be stabilized at a value between R surface and R. In this example, both stabilization values lie between R surface = 1.5 (thin black line) and R = 3.12 (thick black line) (see Table 1 The objective of this paragraph is to investigate the equipartition regime more locally. To this end, an energy analysis is done over some horizontal slices at five different depths, i.e. z = 0, -250, -500, -1000m. The results are depicted in Figure 3 for both low-pass white noise (left plot) and exponential (right plot) correlations. This figure shows the temporal variation of the slice-averaged ratio ||curl u(x)|| 2 slice / (div u(x)) 2 slice . It can be observed that for the slice at the free surface (red curve), the energy ratio converges to the analytical value of K 2 R surface = 2 (thin black line). For other slices, stabilization occurs at a value between K 2 R surface = 2 and 2K 5 = 4.2 (thick black line) but closer to the full-space stabilization value as was also the case for the averaged energies over the whole medium.

               E Ω p (t) = Ω E p (x, t)dx = Ω 3κ(x)+4µ(x) 6 (div u(x)) 2 dx E Ω s (t) = Ω E s (x, t)dx = Ω µ(x) 2 ||curl u(x)|| 2 dx E Ω k (t) = Ω E k (x, t)dx = Ω ρ 2 ∂u(x,t) ∂t 2 dx (19)
In conclusion, for low values of K and a dispersion level of 40%, an equipartition regime is established at global and local scales for different correlation models. Note that according to the relation Σ PS /Σ S P = 2K 3 , lower values of K lead to lower contrasts between Σ PS and Σ S P (irrespective of the source type) which means that during the scattering process, the normal modes exchanges occur with closer probabilities compared to the case of higher values of K. The values of the ratios Σ PP /Σ PS and Σ S S /Σ S P are not significantly far from unity. Furthermore, the values of the transport mean free paths of P and S waves are respectively * P = 435m and * S = 450m so that: • The P and S waves reach a diffusion regime almost simultaneously since * P * S .

• The body waves can reach an equipartition regime before leaving the medium (since * P , * S < 1500m) and at enough distances (at least a transport mean free path) from the boundaries which can a priori pollute the wave energies by the unwanted reflections toward the physical medium.

These are some reasons that favor the onset of an equipartition regime in this case. In Section 5, we will study less favorable cases where this regime does not set in. slice / (div u(x)) 2 slice calculated over the whole slices including the intersections with the PMLs. Red, green, cyan and purple colors correspond respectively to the slices at z = 0m, z = -250m, z = -500m and z = -1000m. The thick and thin black lines correspond respectively to the analytical values on full space and over the free surface. Left and right plots correspond respectively to the low-pass white noise and exponential correlation models.

Limitations to reaching the equipartition regime

In this section we introduce numerical simulations in which an equipartition regime will not be observed. Table 2 lists all configurations that are considered in this paper along with the ratios between the total scattering cross sections for two different correlation kernels in terms of c , δ and K. Note that the "Reference" and "Reference-exp" configurations correspond to the case introduced in Section 4.

Setting up of a localization regime

In this case we increase the value of the mean phase velocity contrast to K = √ 3 (the "Localization" configuration in Table 2). To attain a more favorable case to reach an equipartition regime, one should increase the value of the dispersion level δ. We will consider δ = 0.6. The differences between the values of * P and * S in this case and the corresponding values for the cases with the same K and lower δ (0.4 or 0.15) is remarkable. In the following paragraph, the case of a propagation medium with K = √ 3 and δ = 0.6 will be treated. The half-space Ω is excited with two different source types: an explosion and a unidirectional force, both of which having a Ricker pulse time function. The simulation time is t = 5s and a low-pass white noise model is used as the correlation kernel of the heterogeneities.

In Figure 4, the top left plot depicts the temporal variations of E Ω p and E Ω s using the same color as in the previous section. The solid and dashed curves correspond respectively to the explosion and unidirectional sources. The unidirectional source which creates initially the S wave energies more than the P wave ones, results at first (before arriving to the gray shaded area) in higher values of E Ω s and lower values of E Ω p compared to the case of explosion source. However, at longer lapse times, the body wave energy densities become source-independent. We observe from the left plots of the Figure 4 that regardless of the source type, even before the complete evacuation of the direct wave energies, both of the phase energies are almost stabilized at a lapse time of t = 2s and their decay rate tend to zero at longer lapse times. This is in contradiction with a transport regime wherein E The bottom left plot of the Figure 4 shows the stabilization of the space-averaged energy ratio for both source types. We guess that this stabilization is related to the fact that a major part of the wave energy is trapped within the propagation medium so that it does not reach the boundaries. To verify this hypothesis, we do the same simulation (with an explosion source) but with Neumann boundary conditions all around the medium and the corresponding results are shown with dashed-dotted curves. From the top left plot, we observe a slight increase in E Ω p and E Ω s , resulting in a small change in the stabilization value (dashed-dotted blue curve in lower left plot). The latter is observed to be more than twice the analytical value in full-space. As it has been already discussed, this ratio must be even less than the analytical value 2K 3 which is not the case here. The right plots also show the evolution of E This regime is called the strong (or Anderson) localization regime [START_REF] Weaver | Anderson localization in the time domain: numerical studies of waves in two-dimensional disordered media[END_REF][START_REF] Sheng | Introduction to Wave Scattering, Localization, and Mesoscopic Phenomena[END_REF][START_REF] Tiggelen | Localization of waves[END_REF][START_REF] Friedrich | Localization of seismic coda at merapi volcano (indonesia)[END_REF] in which the random heterogeneities of the medium are too strong (δ = 0.6) so that the body wave energies are locally trapped before they reach the boundaries. In this regime, the interference effects are dominant such that they block the transport of energy in the medium. Hence, in this regime the transport equations can no longer describe properly the multiple scattering phenomenon, [START_REF] Friedrich | Localization of seismic coda at merapi volcano (indonesia)[END_REF]. It should be pointed out that in a localization regime there is no theoretical formula which relates the stabilization value to the statistical parameters of the medium. However, Weaver [START_REF] Weaver | Anderson localization in the time domain: numerical studies of waves in two-dimensional disordered media[END_REF] phenomenologically found an equation to describe the variation of the total energy density in space and time domain in 2D. Sepehrinia et al. [START_REF] Sepehrinia | Numerical simulation of the localization of elastic waves in two-and three-dimensional heterogeneous media[END_REF] also studied numerically the onset of a strong localization regime in 2D and 3D random elastic media. A transition to a localization regime happens when the Ioffe-Regel criterion, i.e. |k| i sc i < 1, i ∈ {P, S } (where sc i , i ∈ {P, S } is the scattering mean free path of mode i) is satisfied. In this simulation, the scattering mean free paths are sc P = 22m and sc S = 13m. Hence the Ioffe-Rege criterion is met for both wave modes and a localization regime is established. The visualization of the energy snapshots will help to observe the localization phenomenon. The snapshots of the adimensional ratios (div u(x)) 2 / (div u(x)) 2 slice and ||curl u(x)|| 2 / ||curl u(x)|| 2 slice calculated over the corresponding slice at the lapse times of t = 2s and t = 4s are shown in Figures 5 and6 respectively. The pattern of these fields is similar to a collection of points distributed over the slices containing the essential part of the slice's average energy. According to the Table 2, the transport mean free paths of the P and S waves are respectively * P = 25m and * S = 130m. The main point in these figures is that the body wave energies are localized within almost the same zone over each slice between the lapse times t = 2s and t = 4s. The results of the slice-averaged energy ratio gives similar information as in the lower left plot of the

Influence of mean free path

The propagation medium Ω is again subjected to an explosion source and the simulation time is considered to be t = 5s. It has a low-pass white noise ACF, the values of K and δ are respectively √ 3 and 0.15, and the medium is a half-space (see the "Large * P/S " configuration in Table 2). Left plots in Figure 7 show the temporal variations of E ) is globally stabilized over the entire medium in some time intervals, a global equipartition regime does not occur (as already seen in the bottom left plot of Figure 7).

From Table 2, the values of the transport mean free paths of P and S waves are respectively * P = 4000m and * S = 6500m. These distances cannot be achieved by the waves before reaching the boundaries in this simulation. A medium whose sides are at least 2 * S = 13000m which is more than 4 times bigger (in each side) than the propagation medium in this case, is required in order to observe the equipartitioning of the energies.

In conclusion, for Poissonian materials (K = √ 3) with low degrees of fluctuation (δ = 0.15 in this case), the equipartition regime was not observed because of the small propagation lengths compared to the transport mean free paths. In the next section, we increase the value of δ from 0.15 to 0.4 with the aim of decreasing the transport mean free paths and make them less than the propagation length.

Influence of inefficient mixture of body waves

In this section, we will be concerned with the influence of inefficient mixture between the body wave energies during the scattering events in the numerical observation of equipartition. For this purpose, the dispersion level of the elasticity matrix is hold to be the same as in the Reference case (δ = 0.4) and the same value of K = √ 3 will be used (the "Inefficient mixture" configuration in Table 2).

The propagation medium Ω is subjected to three different sources: explosion, unidirectional and their combination all of which having a Ricker pulse time function. The medium is a half-space with a low-pass white noise ACF. The simulation is performed over t = 5s. The top left plot in Figure 8 shows the temporal variation of E Ω p (red curves) and E Ω s (blue curves) where the solid, dashed and dashed-dotted curves represent respectively the results relative to the explosion, unidirectional and combined sources. This plot shows that the waves are not localized since both of the body wave energies are decaying with time. The bottom left plot depicts the corresponding S-to-P energy ratios which states no global equipartitioning for all source types. From Table 2, the values of the transport mean free paths of the P and S waves are respectively * P = 60m and * S = 580m. Hence, both wave modes can travel at least a transport mean free path before the arrival of the direct waves to the boundaries. Following Table 2 we have Σ PP /Σ PS = 32 and Σ S S /Σ S P = 33 implying that the P wave energies have much more tendency to scatter into the P waves and likewise the initial S wave energies are more likely to scatter into the S waves. For the sake of comparison, we remind the values of the ratios Σ PP /Σ PS and Σ S S /Σ S P being respectively 0.5 and 3.2 for the Reference case. These values are closer to 1 and imply a faster establishment of an equipartition regime [START_REF] Margerin | Monte carlo simulation of multiple scattering of elastic waves[END_REF]. Moreover, the ratio Σ PS /Σ S P is farther from 1 compared to the Reference case where it is 0.32. As a consequence, the asymmetry between the P-to-S and S-to-P mode exchanges is higher. The fact that the values of these three ratios are not as close to 1 as in the Reference case, can be imagined as another reason of discrepancy between the observed and the analytical stabilization values.

As a result, for K = √ 3 and δ = 0.4 even if the transport mean free paths are small compared to the propagation length, large values of the ratios Σ PP /Σ PS and Σ S S /Σ S P and the fact that the ratio Σ PS /Σ S P is not close enough to 1 prevent the onset of an equipartition regime. In previous sections, we put our focus on the physical parameters related to the propagation medium and no attention has been paid to the influence of the numerical parameters such as the boundary conditions to the stabilization values. In the following section, we discuss this influence via numerical simulations.

Influence of the Perfectly Matched Layers (PMLs)

The propagation domain being subjected to an explosion has now three different boundary conditions of types half-space (PMLs all around the medium except over the free surface), full-space (PMLs surrounding the propagation medium) and Neumann (reflecting or stress-free boundaries). The medium has the same statistical properties as that introduced in Section 5.3. Figure 9 summarizes all the results regarding the spatially-averaged energy densities. The results corresponding to the half-space, full-space and Neumann cases are differentiated by solid, dashed and dasheddoted curves. The lower left plot shows the corresponding stabilization values. A direct result that can be drawn from this plot is that the existence of the PMLs implies an increase in the stabilization values. This is assumed to be an effect of PML boundary conditions which seems to behave in favor of the S wave energy. We assume that the PML domains absorb the P waves more than the S waves and the small reflections towards the medium by them are almost the S waves. However, this observation could not necessarily be generalized to all numerical case studies. In these simulations, the stabilization values are about 14.5, 14 and 12 respectively for half-space, full-space and Neumann boundary condition cases. These values are not theoretically acceptable since they do not belong to the interval [6.44, 10.4]. This example shows that the use of PMLs, and probably of any other kind of absorbing layer or boundary, clearly has an influence of the equipartition regime. This should obviously be kept in mind when designing numerical experiments to observe the transition to the equipartition regime.

Conclusion

In this paper the influence of some physical and numerical parameters on the onset of an equipartition regime is investigated. In these analyses and through the numerical simulations, we identified the stabilization of the energy ratios along with the time at which this stabilization occurs. We showed theoretically that the correlation model of the underlying randomly heterogeneous medium has a noticeable influence, especially at high frequencies. Then, four different case studies were considered in which we change the physical parameters. The main results that have been extracted are: (1) For low values of K, the dispersion level of the random medium can be up to 40% in order to be capable of observing an equipartition regime as contrary to the theory which imposes low variances. (2) For higher values of K (close to the reference value √

3) and high values of δ (close to the maximum theoretically acceptable value, i.e. 70%) a localization regime is dominant because of high strengths of fluctuations. Hence, transition to an equipartition regime is impossible. (3) For higher values of K and medium dispersion levels (around 40% for instance), due to the inefficient mixture between the body waves, it is impossible to observe an equipartition regime. (4) For higher values of K and low dispersion levels (less than 15%), an equipartition regime could be observed but the large values of the transport mean free paths imply large sizes for the medium and hence high computational costs. At the end, we observed that the influence of the PMLs in our simulations is in favor of the S wave such that they increase the stabilization ratios.

Although Radiative Transfer and Diffusion equations provide very interesting insights into the behavior of waves in complex heterogeneous media and coda waves, their study is theoretically limited to full spaces, and to some extent, to half spaces. For more complex problems, numerical simulation is an interesting alternative. However, comparison with asymptotic analyses is difficult, because the latter are valid when given ratios of parameters (wavelength to correlation length for instance) are very small, whereas the exact meaning of "very small" is rarely discussed. This paper tried to provide quantitative guidelines on physical and numerical parameters for which equipartition can be observed, beyond the "very small" or "very large" given by the literature. By helping in understanding these questions, the paper participates in building a tool that can be used in the future to understand and model seismic coda for more real media. While the paper considered the controlled setting of full and half-spaces (if it were not for the PMLs), future work will obviously introduce more complicated geometries (topography and discontinuous material properties).

3 Figure 1 :

 31 Normalized global equipartition time t eq (|k|) = t eq (|k|) v p / c in terms of |k| c for K = 1.16 (a) and K = √

  in which Ω is the whole propagation medium excluding the PMLs. The individual variation of E Ω p (red curves), E Ω s (blue curves) are depicted in the top left plots. The solid and dashed curves in this Figure indicate the results corresponding to the medium with low-pass white noise and exponential correlation model, respectively.

  ) and are closer to the analytical value for an open unbounded medium which implies that the body waves are dominant over the entire medium. The top right plot of the Figure2

Figure 2 :

 2 Figure 2: Top left: temporal variation of E Ω p (red curves) and E Ω s (blue curves). Bottom left: temporal variation of E Ω s /E Ω p . Top right: temporal variation of KΩ (blue curves) and E Ω p + E Ω s (red curves). Bottom right: temporal variation of KΩ /(E Ω p + E Ω s ). Solid and dashed curves correspond respectively to the low-pass white noise and exponential correlation models. The starting and end points of the shaded window indicate respectively the lapse times in which the direct P and S waves propagating in homogeneous background arrive to the boundaries. shows that the space-averaged kinetic energy densities E Ω k (t) (blue curves) reach their maximum value at the source's delay time t = t 0 = 0.3s. From the bottom right plot, stabilization values of 0.8 and 1.1 are observed respectively for the low-pass white noise and exponential correlation models. The corresponding stabilization times are observed to be highly dependent on the correlation model. In this example, both of the ratios E Ω s /E Ω p and E Ω k /(E Ω p + E Ω s ) are stabilized for both correlation types and hence an equipartition regime is established.The objective of this paragraph is to investigate the equipartition regime more locally. To this end, an energy analysis is done over some horizontal slices at five different depths, i.e. z = 0, -250, -500, -1000m. The results are depicted in Figure3for both low-pass white noise (left plot) and exponential (right plot) correlations. This figure shows the temporal variation of the slice-averaged ratio ||curl u(x)|| 2 slice / (div u(x)) 2 slice . It can be observed that for the slice at the free surface (red curve), the energy ratio converges to the analytical value of K 2 R surface = 2 (thin black line). For other slices, stabilization occurs at a value between K 2 R surface = 2 and 2K 5 = 4.2 (thick black line) but closer to the full-space stabilization value as was also the case for the averaged energies over the whole medium.In conclusion, for low values of K and a dispersion level of 40%, an equipartition regime is established at global and local scales for different correlation models. Note that according to the relation Σ PS /Σ S P = 2K 3 , lower values of K lead to lower contrasts between Σ PS and Σ S P (irrespective of the source type) which means that during the scattering process, the normal modes exchanges occur with closer probabilities compared to the case of higher values of K. The values of the ratios Σ PP /Σ PS and Σ S S /Σ S P are not significantly far from unity. Furthermore, the values of the transport mean free paths of P and S waves are respectively * P = 435m and * S = 450m so that: • The P and S waves reach a diffusion regime almost simultaneously since *

u|| 2 /

 2 (Div u)2

Figure 3 :

 3 Figure 3: Temporal variation of ||curl u(x)|| 2slice / (div u(x)) 2 slice calculated over the whole slices including the intersections with the PMLs. Red, green, cyan and purple colors correspond respectively to the slices at z = 0m, z = -250m, z = -500m and z = -1000m. The thick and thin black lines correspond respectively to the analytical values on full space and over the free surface. Left and right plots correspond respectively to the low-pass white noise and exponential correlation models.

Ω

  p and E Ω s should decay in time because of the existence of PML domains around the medium.

s

  and their ratio for all above-mentioned cases. The top right figure shows that after a certain time, the values of E Ω s and E Ω p + E Ω s become time-independent for all three cases and their ratio tends to 0.35. This convergence occurs as soon as the direct P waves begin to leave the medium.

Figure 4 :

 4 Figure 4: Top left: temporal variation of E Ω p (red curves) and E Ω s (blue curves). Bottom left: temporal variation of E Ω s /E Ω p . Top right: temporal variation of KΩ (blue curves) and E Ω p + E Ω s (red curves). Bottom right: temporal variation of KΩ /(E Ω p + E Ω s ). Solid, dashed and dashed-dotted curves correspond respectively to half-space with explosion source, half-space with unidirectional source and Neumann boundary conditions with explosion source. The starting and end points of the shaded window indicate respectively the lapse times in which the direct P and S waves propagating in homogeneous background arrive to the boundaries.

Figure 4 .

 4 The objective of the next section is to deal with a simulation in which we keep the value of K = √ 3 but we decrease δ to 0.15 with the aim of getting rid of the localization regime which was dominant in this section. (a) t = 2(s), z = 0(m) (b) t = 4(s), z = 0(m) (c) t = 2(s), z = -1000(m) (d) t = 4(s), z = -1000(m) Figure 5: Snapshots of the adimensional ratio (div u(x)) 2 / (div u(x)) 2 slice in percent for two different lapse times t = 2(s) (plots (a,c)) and t = 4(s) (plots(b,d)) and two different depths z = 0(m) (plots (a,b)) and z = -1000(m) (plots (c,d)).

  (a) t = 2(s), z = 0(m) (b) t = 4(s), z = 0(m) (c) t = 2(s), z = -1000(m) (d) t = 4(s), z = -1000(m) Figure 6: Snapshots of the adimensional ratio ||curl u(x)|| 2 / ||curl u(x)|| 2 slice in percent for two different lapse times t = 2(s) (plots (a,c)) and t = 4(s) (plots(b,d)) and two different depths z = 0(m) (plots (a,b)) and z = -1000(m) (plots (c,d)).

Ω

  p and E Ω s and their ratio. From the bottom left plot, it can be observed that until about t = 2s (the end of the gray shaded area), the ratio E Ω s /E Ω p tends to increase. At this time, both of the direct P and S waves have already arrived to the boundaries and the role of the boundary conditions begins. The top right plot of the Figure7shows the temporal

Figure 7 :

 7 Figure 7: Top left: temporal variation of E Ω p (red curves) and E Ω s (blue curves). Bottom left: temporal variation of E Ω s /E Ω p . Top right: temporal variation of KΩ (blue curves) and E Ω p + E Ω s (red curves). Bottom right: temporal variation of KΩ /(E Ω p + E Ω s ). The starting and end points of the shaded window indicate respectively the lapse times in which the direct P and S waves propagating in homogeneous background arrive to the boundaries.

Figure 8 :

 8 Figure 8: Top left: temporal variation of E Ω p (red curves) and E Ω s (blue curves). Bottom left: temporal variation of E Ω s /E Ω p . Top right: temporal variation of KΩ (blue curves) and E Ω p + E Ω s (red curves). Bottom right: temporal variation of KΩ /(E Ω p + E Ω s ). Solid, dashed and dashed-dotted curves correspond respectively to explosion, unidirectional and combined source. The starting and end points of the shaded window indicate respectively the lapse times in which the direct P and S waves propagating in homogeneous background arrive to the boundaries.

Figure 9 :

 9 Figure 9: Top left: temporal variation of E Ω p (red curves) and E Ω s (blue curves). Bottom left: temporal variation of E Ω s /E Ω p . Top right: temporal variation of KΩ (blue curves) and E Ω p + E Ω s (red curves). Bottom right: temporal variation of KΩ /(E Ω p + E Ω s ). Solid, dashed and dashed-dotted curves correspond respectively to half-space, full-space and Neumann boundary conditions. The starting and end points of the shaded window indicate respectively the lapse times in which the direct P and S waves propagating in homogeneous background arrive to the boundaries.

Table 1 :

 1 Analytical values of the equipartitioning ratio for unbounded media and over the medium's traction-free surface for different values of K /( E p + E s ) tend to 1. The corresponding stabilization values for this ratio with different hypotheses for Poissonian materials are summarized in

		K	K 2 r	vr vs	r	R	R surface
	4 3	1.16 0.49	0.7	1.06 3.12	1.5
	√	1.5	0.80 0.89 1.34 6.75	4.02
	3 1.73 0.85 0.92 1.47 10.4	6.44
		2	0.87 0.93 1.57	16	9.80
	ratio E k

Table 2 :

 2 Scattering parameters for a random medium with low-pass white noise and exponential correlation kernels Config.

	P [m]

c [m] δ K Σ PS /Σ S P Σ PP /Σ PS Σ S S /Σ S P * S [m] *

http://www.ipgp.fr/ ~paulcup/RegSEM.html

Acknowledgments

The authors would like to thank Paul Martin Mai (KAUST) and the anonymous reviewer for their constructive comments that greatly contributed to the improvement of the final version of the paper. This work done within the SINAPS @ project benefited from French state funding managed by the National Research Agency under program RNSR Future Investments bearing reference No. ANR-11-RSNR-0022-04.