Etude des propriétés élastiques de la croûte

Measuring the scattering mean free path of Rayleigh waves on a volcano from spatial phase decoherence 

In heterogeneous media, after a sufficient amount of time of propagation, waves enter the multiple scattering regime. In this regime waves bounce on several heterogeneities before reaching the receivers. The characteristic length (resp. time) after which such a regime can be observed is the scattering mean free path ` (resp. time), defined as the distance (resp. time) between two successive scattering events. This distance also refers to the characteristic distance of attenuation of the coherent wavefront. The coherent wave, rigorously defined as the wave that resists ensemble averaging, roughly corresponds to the direct (or ballistic) wave. The scattering mean free path ` reflects the degree of heterogeneity of the medium: the longer the scattering mean free path, the weaker the scattering. ` depends on two features: the intensity of the fluctuations of the mechanical properties in the medium, and the spatial extension of the fluctuations. In nature, ` is found to vary over several orders of magnitude, depending on the frequency, and also on the nature of the material at test. The multiple scattering regime can be observed in many fields of wave physics. It was for instance demonstrated in optics and in acoustics with the observation of the coherent backscattering effect (or weak localization) (van Albada and Lagendijk, 1985; Wolf and Maret, 1985; Tourin et al., 1997). Seismic waves are also known to exhibit long lasting wave-trains that follow ballistic P- or S-waves: the so-called seismic coda. Since the pioneering works of Aki (1969), the coda is known to be reproducible and has been recognized to possibly originate from multiple scattering effects. Among other 4.1 Introduction applications, scattering parameters are excellent candidates for remotely characterizing geological media at depth, which is a key challenge in geosciences (Margerin and Nolet, 2003). Nevertheless, in most practical cases it is very hard to discriminate scattering effects from intrinsic absorption effects. For instance, the Spatial Auto-Correlation (SPAC) technique, widely used in near-surface geophysics (Aki, 1957), consists in fitting the seismic wave spatial-correlation by the product of a Bessel function and an exponential decay whose physical interpretation remains debated (Prieto et al., 2009; Tsai, 2011; Nakahara, 2012). At longer distance from the source, it is also possible to study the envelope of the diffusive coda, which in principle allows one to evaluate the scattering and intrinsic attenuation (Hoshiba, 1993; Carcolé and Sato, 2010). In continental areas however, an unbiased estimation of these two parameters may become difficult when the leakage of scattered waves is dominant (Margerin et al., 1999). Although a possible remedy has been recently proposed (Del Pezzo and Bianco, 2010), evaluating the scattering mean free path of a complex material without the bias of intrinsic absorption remains a challenging issue. Recently, Anache-Ménier et al. (2009) studied the phase fluctuations in the coda of earthquakes recorded during a temporary experiment at the Pinyon Flats Observatory in California. They proved that seismic coda waves obey Gaussian statistics in certain frequency bands, and suggested to use the correlation of the spatial phase derivative measured in the coda to estimate the scattering mean free path. The key point of their approach is that the phase fluctuations are caused by random phase shifts acquired at each scattering event, and are therefore independent of the absorption structure. Note that the quantity of physical interest is the spatial phase difference (or derivative) and not the phase itself, because the latter is dominated by a trivial ωt dependence, with ω the circular frequency. The determination of the scattering mean free path at Pinyon Flats was not completely conclusive because the aperture of the experimental network was less than one wavelength. In this article, we pursue the work of Anache-Ménier et al. (2009) and communicate the results of a temporary field experiment that was specially designed to measure the scattering mean free path from spatial phase decoherence. The field experiment was set up at the foot of a recent (late pleistocene & holocene) but inactive volcano in the French Auvergne in 2010. A volcanic area is particularly convenient for the study, as areas with tectonic and/or volcanic activity are highly heterogeneous and known to produce long lasting coda waves (Abubakirov and Gusev, 1990; Aki and Chouet, 1975; Aki and Ferrazzini, 2000; Goodman, 1985). Larose et al. (2004a) studied the weak localization of seismic waves at this site and found an estimate of 200 m for the mean-free path for seismic waves around 20 Hz. This gives us a typical order of magnitude of the scattering mean free path that we can expect in our experiment. In section 4.2 we first give an overview of the theoretical aspects that underlie our study. Then, in section 4.3 we describe the field experiment; prove that our signal obeys  circular Gaussian statistics, the requirement to determine the scattering mean free path ` from spatial phase decoherence; and show the good estimate of ` that we obtain with our experiment (section 4.3.5). 

Theoretical aspects 

Definition of the phase 

To define unambiguously the phase of the seismic wavefield u(t, r) recorded at position r and time t, we introduce the associated analytic signal as follows: Ψ(t, r) = u(t, r) + iHu(t, r), (4.1) where i is the imaginary unit and Hu(t) = P.V. R u(t 0 )dt0 t−t 0 denotes the Hilbert transform of the field u(t). Using the polar representation of complex numbers, the field Ψ can be expressed as: Ψ(t, r) = A(t, r)exp(iΦ(t, r)) (4.2) where Φ(t, r) is the wrapped phase of the field which takes values in (−π : π]. By correcting the phase Φ for the 2π jumps that occur at ±π, one obtains the unwrapped phase Φu, which is a continuous function taking values in R. The unwrapping operation can be performed either in the time or in the spatial domain. In this work, the field is analyzed on a linear array of seismic stations, which calls for the spatial unwrapping of the phase. In general, the unwrapping operation is not topologically invariant, i.e., it depends on the path from the initial to the final point. In our experiment, the path is dictated by the linear geometry of the array. Note that the phase of narrowly bandpassed signals is dominated by the term ωt -ω the central frequency-, which does not convey any interesting information on the medium. By evaluating the phase difference between two nearby stations, this trivial ωt dependence is removed and the interesting fluctuations of the phase caused by the presence of heterogeneities in the medium become accessible. Motivated by this observation, we introduce two possible definitions of the phase difference. The first one, denoted by ∆Φ ∈ (−2π, 2π], is obtained by subtracting the wrapped phases measured at two adjacent stations. The second one, denoted by ∆Φu, is obtained by unwrapping the phase spatially at each time step. Due to the finite separation between stations, it is impossible to distinguish between a large physical jump of the phase from an artefact caused by its mathematical definition. To circumvent this cycle skipping problem, we impose that the absolute value of a phase jump between two nearby stations cannot exceed π, i.e. ∆Φu takes values in (−π, π]. Using this convention, the following relations between the two definitions of the phase difference can be established: ∆Φu =∆Φ, ∆Φ ∈ (−π, π] ∆Φu =∆Φ − 2π, ∆Φ ∈ (π, 2π] ∆Φu =∆Φ + 2π, ∆Φ ∈ (2π, −π] (4.3) 46 4.2 Theoretical aspects 

 Circular Gaussian Statistics 

Although Gaussianity is a standard hypothesis in statistical wave propagation problems, we give some heuristic arguments in support of this assumption. In the multiple scattering regime, the field u measured at a point can be considered as a superposition of a large number of partial waves ψα that propagated along independent paths: u = ReX α ψα = ReX α aαexp(iφα), (4.4) where the subscript α labels the different trajectories. The hypothesis that the partial waves follow independent paths is valid when the average distance between two scattering events, i.e. the mean free path `, is much larger than the wavelength λ. Due to the phase shift that occurs at each scattering event, after a few scattering mean free paths, the phases of the partial waves become random and uniform in the interval (−π : π]. In these conditions, we may apply the central limit theorem, which stipulates that the arithmetic mean of a large number of identically distributed and independent random variables will be approximately normally distributed. Hence, we conclude that the seismic coda wavefield u(t) recorded at an arbitrary point of the medium can be modeled as a centered Gaussian random variable with variance σ 2 (t) = I(t), where I(t) is the intensity of the coda. Because the Hilbert transform is a linear operator, it can be demonstrated that the imaginary part of the associated analytic signal obeys the same Gaussian distribution and is independent from the real part(Goodman, 1985). A complex random variable which verifies the properties just enunciated is known as gaussian circular. Based on this observation, we now assume that the complex analytic signals recorded at an arbitrary number of points at time t in the coda are jointly gaussian circular: P(u1, v1, . . . , uN , vN ) = P(Ψ1, . . . , ΨN ) = 1 π N det(C) exp    −(Ψ∗ 1 . . . Ψ ∗ N )C −1   Ψ1 . . . ΨN       , (4.5) where P denotes the probability density, u1 = ReΨ1, v1 = ImΨ1, and C = D ΨiΨ∗ j E is the covariance matrix. (Ψi . . . ΨN ) are N different measurements of the field after application of a normalization procedure to be detailed below. In this work, we do not test the validity of Eq. (4.5) in its full extent. Instead, we will derive marginal distributions for the phase and its derivative based on assumption (4.5), and compare the experimental measurements with the theoretical predictions to obtain estimates of the wavenumber and mean free path under the array. We also assume that different times in the coda correspond to different realizations of the underlying random process and that temporal and statistical averaging are equivalent (ergodic hypothesis). To complete our definition of the statistical properties of the wavefield, the covariance matrix C must be specified. Each element of the matrix depends on the two-point correlation function of the recorded wavefield u, narrowly band-passed around circular frequency ω. In the multiple-scattering regime, the field-field correlation function in the   coda may be approximated as follows: hu(R + x/2, t)u(R − x/2, t)i = S(ω)e −R2/4Dt−t/ta (4πDt) d/2 ImGd(ω, |x|), (4.6) where S(ω) is proportional to the source spectrum, R is the position vector connecting the source to the midpoint of the two receivers, D is the diffusion constant of the waves, ta is their absorption time, and Gd is the average Green’s function of the multiple-scattering medium. The subscript d represents the space dimension and | · | the corresponding Euclidean distance. In our experiment u represents the vertical component of the ground displacement. Formula (4.6) has been demonstrated for scalar waves by (Barabanenkov and Ozrin, 1991) based on an eigenfunction expansion of the Bethe-Salpeter equation. Extensions to vectorial fields such as electromagnetic waves or coupled P and S waves have also been published (Barabanenkov and Orzin, 1995; Margerin, 2013). Physically, Eq. (4.6) expresses the diffusive transport of the energy from the source to the array, and displays explicitly the proportionality between field-field correlations and the imaginary part of the Green’s function. In seismology, the successful extraction of the surface wave part of Green’s function between two stations from coda waves has been previously reported by (Campillo and Paul, 2003; Paul et al., 2005). A very important point to be noted is that the Green’s function which appears in Eq. (4.6) depends solely on the mean free path, at least for sufficiently weak absorption. We must now determine the correct form of Gd to be inserted in Eq. (4.6). From our previous coherent back-scattering experiment in the same area (Larose et al., 2004a), we concluded on the dominance of Rayleigh waves in the coda. In particular, by measuring the lateral extension of the backscattering enhancement spot as a function of frequency, we were able to estimate the dispersion law of Rayleigh waves, in excellent agreement with independent estimates based on classical seismic techniques. After normalization of the two-point correlation function by the intensity received in the coda, we can write the elements of the covariance matrix as: hΨiΨ ∗ j i = C¯(|xj − xi |) = J0(k|xj − xi |)e −|xj−xi|/2l , (4.7) where k and ` denote the wavenumber and mean free path of the Rayleigh waves under the array, respectively. The function C¯ is recognized as the normalized two-point correlation function of the wavefield. The underlying 2-D character of the propagation is apparent in Eq. (4.7). Note that the details of the intensity normalization are unimportant because the distribution of energy in the coda is spatially homogeneous at the scale of the linear array. While relation (4.7) gives in principle access to the mean free path `, the oscillatory character of the Bessel function J0 makes it difficult to extract the rate of decay of the exponential term in practice. The key of the method proposed by Anache-Ménier et al. (2009) is to remove the oscillatory term by considering the correlation of the phase derivative. The main properties of the phase and of its derivative are therefore outlined in the next paragraph. 

 Statistical properties of the phase 

We first discuss briefly the one-point statistics (N = 1). We remind the reader that the real and imaginary part of the field are independent and normally distributed with zero mean and equal variance σ 2 = I(t), which depends on the time in the coda. Although our paper focuses on the properties of the phase, it is worth recalling the marginal probability densities of other quantities derived from the complex field Ψ, as they are readily measured experimentally: • The amplitude A follows a Rayleigh distribution: P(A) = 2A hIi exp − A2 hIi ! , (4.8) where hIi = hA2 i is the average intensity. Using a large dataset from northern Japan, Nakahara and Carcolé (2010) have demonstrated that the amplitude of coda waves is described by Rayleigh statistics with very good accuracy. • The distribution of intensity is described by a decaying exponential law: P(I) = 1 hIi exp − I hIi ! . (4.9) which follows straightforwardly from Eq. (4.8) after introduction of the new variable I = A2 . • Finally, as previously remarked, the phase Φ exhibits a uniform probability distribution: P(Φ) = 1 2π . (4.10) As any oscillating signal shows a uniform phase distribution, the last property does not constitute a conclusive test of Gaussianity. In addition, we wish to point out that a wavefield obeying one-point circular Gaussian statistics does not necessarily result from multiple scattering. On the one hand, in a single scattering medium with low-contrast inclusions, the phase of the waves is not significantly modified upon scattering, and remains essentially equal to ωt, which annihilates the independent and random character of the phase. On the other hand, in a dilute scattering medium containing high-contrast inclusions, each scattering event introduces a possibly large phase shift which depends on the scattering angle, the geometry, and the physical properties of the obstacle; since the phase of the singly- or doubly- scattered partial waves may also be considered as random, the resulting wavefield may obey Gaussian statistics and exhibit fluctuations similar to those observed in the high-order multiple scattering regime. We now explore in more details the two-point statistics (N = 2). Consider two wavefields (Ψ1 = A1e iΦ1 , Ψ2 = A2e iΦ2 ), recorded at geophones separated by a distance δ. Our objective is to derive the statistical properties of the phase difference ∆Φ = Φ2−Φ1, of the unwrapped phase difference ∆Φu, and of the phase derivative Φ 0 . The correlation matrix has the following simple form: C = 1 g g 1 ! , (4.11) where 0 < g = C¯(δ) < 1. The determinant and inverse of this matrix are readily calculated. Inserting these results into Eq. (4.5), one obtains the joined distribution of the two fields: P(Ψ1, Ψ2) = e −|Ψ1| 2−|Ψ2| 2+2gReΨ1Ψ∗ 2 π 2 (1 − g 2 ) . (4.12) Introducing the new variables (A1, A2, Φ1, Φ2) and integrating out the phase Φ1 yields: P(A1, A2, ∆Φ) = A1A2(2π − |∆Φ|) π 2 (1 − g 2 ) e −A2 1 − A2 2 + 2A1A2g cos ∆Φ 1 − g 2 , (4.13) where ∆Φ ∈ (−2π, 2π]. To eliminate the remaining variables A1 and A2, it is convenient to introduce polar coordinates A1 = r cos θ, A2 = r sin θ with θ ∈ [0, π/2]. The integration over (r, θ) is a straightforward computational exercise and the final result is: P(∆Φ) = N(1 − g 2 ) 1 − f 2 1 + f arccos −f √ 1 − f 2 ! , (4.14) where N = (2π − |∆Φ|)/4π 2 . Note that the distribution of the unwrapped phase difference ∆Φu is easily obtained using relations (4.3), and is formally identical to that of ∆Φ with N = 1/2π. In the limit δ → 0, we may expand the phase difference ∆Φ ≈ Φ 0 δ and the correlation coefficient g ≈ 1 − Qδ2/2, where Φ 0 represents the first spatial derivative of the phase and Q = −C¯00(0). Note that since the correlation function C¯ is even, the correlation coefficient g must be expanded to second-order in the small parameter δ. Upon reporting these approximations into Eq. (13), and applying a Taylor series expansion around δ = 0, we obtain the statistics of the first phase derivative: P(Φ0 ) = Q 2(Q + Φ02 ) 3/2 (4.15) Remarkably, the first spatial derivative of the phase depends on only one parameter Q. Using the expression (4.7) for C¯, we obtain Q = k 2/2(1 − 2/(kl) 2 ) ≈ k 2/2. The last approximation is valid because the mean free path is in general much larger than the wavelength. This means that within the experimental accuracy, the measurement of Q does not put constraints on the mean free path. For sufficiently large values of the phase derivative, the probability density (4.15) is independent of the parameter Q and follows a power law decay P(Φ0 ) ∼ Φ 0−3/2 , characteristic of circular gaussian statistics. We now outline the derivation of the correlation function of the spatial phase derivative along the array. The method is analogous to the one employed to derive Eq. (4.14) 50 4.2 Theoretical aspects but the calculations are far more complex. We provide the key ingredients and refer the interested reader to the work of van Tiggelen et al. (1999) for more details. The first step is the calculation of the correlation matrix of four fields Ψi (i = 1, 2, 3, 4) acquired at four positions: x1, x1 + δ, x3, x3 + δ. Using expression (4.7), we obtain: C(x, δ) =   1 C¯(δ) C¯(x) C¯(x + δ) −C¯(δ) 1 C¯(x − δ) C¯(x) −C¯(x) −C¯(x − δ) 1 C¯(δ) −C¯(x + δ) −C¯(x) −C¯(δ) 1   , (4.16) where x = x3 − x1. Using a perturbation approach, the eigenvectors and associated eigenvalues of the matrix (4.16) are calculated to first and second order in the small parameter δ, respectively. In this way, the joint probability distribution of the fields and of their spatial derivatives (Ψ(x1),Ψ(x3),Ψ0 (x1),Ψ0 (x3)) can be calculated. Introducing the polar representation Ψ(xi) = A(xi)e iΦ(xi) , the joint probability distribution P(A1, A0 1 , A3, A0 3 , Φ1, Φ 0 1 , Φ3, Φ 0 3 ) of the amplitude and phase at x1 and x3 together with their spatial derivatives is obtained. Note that this probability density depends on x = x3 − x1 only. After integration over the amplitudes, amplitude derivatives and phases, the correlation function of the phase derivative may be expressed as (van Tiggelen et al., 2006): CΦ0(x) = ZZ +∞ −∞ P(Φ0 1 , Φ 0 3 )Φ0 1Φ 0 3dΦ 0 1dΦ 0 3 = 1 2 log(C¯(x))00(log(1 − C¯(x) 2 )) (4.17) For sufficiently large x, the exact expression (4.17) can be approximated as follows: CΦ0(x) ≈ e −|x|/l 2|x| (4.18) As anticipated, the oscillatory term J0(k|x|) does not show up in equation (4.18). Using numerical simulations, Anache-Ménier et al. (2009) have demonstrated that the convenient formula (4.18) applies for x > λ/5. Hence, correlations of the spatial derivative of the phase provide direct access to the mean free path. Two very important remarks must be made at this point. To derive Eq. (4.18), we have implicitly assumed that translational invariance applies. If the statistical properties vary under the array, the result (4.18) may well be invalid. The second point pertains to the evaluation of the spatial derivative of the field. Throughout the derivation of Eq. (4.18), we have assumed that the derivatives are evaluated in the direction of the array. The correlation function of the directional derivative of the phase at an angle α with respect to the array direction gives a quite different result. In particular the undesirable oscillations of the correlation C¯ are not suppressed in this case (Ghysels, 2005). 

Experiment and results 

 Experiment setup and data acquisition 

We work at the foot of the neighboring volcanic craters “Grand Sarcoui” and “Puy des Goules” in the French Auvergne. Larose et al. (2004a) could show the evidence ofthe weak localization of seismic waves in this area (at 20Hz) and gave hence a proof of multiple scattering. To measure the scattering mean free path from phase decoherence, we set up an experiment that simultaneously fulfills two conditions: (1) The intergeophone spacing must be sufficiently small to ensure a good correlation of the signals at neighboring geophones; (2) The total aperture of the array should not be too small compared to the mean free path. From numerical simulations with a finite difference code (Derode et al., 2003b), we observed that the intergeophone spacing δ should be δ ≤ λ/10, where λ is the Rayleigh-wavelength, and that an increasing number of geophones significantly increases the precision of the determination of `. We see that it is essential to know the Rayleigh-wavelength in the area of interest. Therefore, in a preliminary study, we shot 20 refraction seismic profiles (hammer source) in the area of interest. We then computed the dispersion curves: once by determining the wavenumbers from the field correlation (Eq. 4.7) at different central frequencies, and once from a frequency-wavenumber analysis. From these analysis, we could determine an approximated wavelength of 24 m at 10 Hz, 11 m at 20 Hz and 7.5 m at 30 Hz. We hence decided to work at high-frequencies (30 Hz) and set up a linear array of 48 geophones along an unpaved road at the foot of the volcanoes with an inter-geophonespacing of δ =60 cm (≈ λ/10). The total length of the array is 28.8 m, which is approximately 4 times the dominant wavelength at 30 Hz. To decrease the ambient noise, we buried the geophones at 20 cm depth. We placed 18 explosives on unpaved roads around the receiver array to excite energy from all directions. The experimental setup is shown in Figure 4.1. The acquired data show long lasting coda waves with a high signal to noise ratio for the first 20 s in the coda (Fig. 4.2). For our analysis we assume a 2D medium. This approach is valid as long as the vertical component of our signals is largely dominated by Rayleigh-waves (Larose et al., 2004a). 

 Test of one-point Circular Gaussian Statistics 

To begin with, we test that the seismic coda waves that we acquired obey Circular Gaussian Statistics. Therefore, we filter the data with a narrow second-order bandpass filter of 5% around a central frequency of 30 Hz. We selected a short time window of 1.5 s starting after the surface wave arrivals. We compute the imaginary part (Hilbert transform) of the recorded signals u(t, r) that gives us access to the phase from the complex analytical signals (Eq. 4.4). We can define: 52 4.3 Experiment and results Figure 4.1: Area of interest in Auvergne. The blue line marks the array of 48 geophones with 0.6 m spacing. The red stars mark the positions of the 18 explosive sources placed around the geophones. 1. The wrapped phase Φ, which is the argument of the complex field u in the range (−π, π]. 2. The unwrapped phase Φu, that is obtained by correcting for the 2π jumps that occur when Φ goes through ±π. The result is a continuous function. We then test the criteria for Circular Gaussian Statistics, as mentioned in section 4.2.2. The results for the different tests are shown in Fig. 4.3 together with the theoretical predictions. We can see that: • The real and imaginary part of the field follow a Gaussian distribution (Fig. 4.3a,b) and are independent (Fig. 4.3c). • The amplitude presents a Rayleigh distribution (Fig. 4.3d). There are only small discrepancies for larger amplitudes that result from superimposed noise. • The intensity is calculated as the squared amplitude of the Hilbert transformed signal. Its distribution clearly follows a decaying exponential probability function (Fig. 4.3e). 50 5 10 15 20 −2 0 2 (a) time (s) Amplitude 0 5 10 15 20 10−4 10−2 100 (b) time (s) Energy Figure 4.2: (a) Raw data from an explosive source, showing a coda of 20 s. The square indicates the chosen time window in the coda. (b) Dimensionless energy of the signal. • The probability distribution of the phase is uniform with a value of 1/2π (Fig. 4.3f). We mentioned in section 4.2.2 that the phase distribution is an ambiguous test for the Gaussian character of a signal, as every oscillating signal shows a uniform phase distribution. We will hence also test the phase derivative distribution to verify that the seismic coda obeys Circular Gaussian Statistics.