Table Of Content

Phase statistics of seismic coda waves D. Anache-Ménier and B. A. van Tiggelen Laboratoire de Physique et de Modélisation des Milieux Condensés, Université Joseph Fourier/CNRS, BP 166, 38042 Grenoble, France L. Margerin Centre Européen de Recherche et d’Enseignement des Géosciences de l’Environnement, Université Aix Marseille, CNRS, Aix en Provence, France We report the analysis of the statistics of the phase fluctuations in the coda of earthquakes 9 recordedduringatemporaryexperimentdeployedatPinyonFlatsObservatory,California. Theob- 0 serveddistributionsofthefirst,secondandthirdderivativesofthephaseintheseismiccodaexhibit 0 universalpower-lawdecayswhoseexponentsagreeaccuratelywithcircularGaussianstatistics. The 2 correlation function of thespatial phasederivativeis measured and used toestimate themean free path of Rayleigh waves. n a J PACSnumbers: 46.65.+g,91.30.Ab,46.40.Cd 9 ] Intheshort-periodband(>1Hz),ballisticarrivalsof 50 h seismicwavesareoftenmaskedbyscatteredwavesdueto p small-scale heterogeneities in the lithosphere. The scat- 0 - tered elastic waves form the pronounced tail of seismo- o e gramsknownas the seismiccoda [1, 2]. Evenwhenscat- −50 g teringisprominent,itis stillpossibleto definethe phase cs. osifgtnhaelsψe(istm,ri)c=recAo(rtd,rb)yeiiφn(tt,rro)d,uwciitnhgAthtehceomampplelxituandaelyantidc y (m)−100 si φ the phase. In the past, many studies have focused on −150 y the modeling of the mean field intensity I(t) = A(t)2 h h i −200 [see 3, for review]. The goal of the present paper is to p study the statistics of the phase field in the coda. In the [ −250 coda, the measureddisplacements result fromthe super- −100 0 100 200 300 2 position of many partial waves which have propagated x (m) v alongdifferentpathsbetweenthesourceandthereceiver. 6 FIG. 1: Geometry of the seismic array Eachpathconsistsofasequenceofscatteringeventsthat 5 affect the phase of the corresponding partial wave in a 5 randomway. Fornarrow-bandsignals,thephasefieldcan 3 . therefore be written as φ(t,r)=ωt+δφ(t,r), where ω is configuredasagridandtwoorthogonalarmswithsensor 9 thecentralfrequency,andδφdenotestherandomfluctu- spacingsof7meterswithinthegridand21metersonthe 0 ations. Thetrivialcyclicphaseωtcancelswhenaspatial arms[5]. Weselected8earthquakesofmagnitudegreater 8 0 phase difference is considered between two neighbouring than 2 with good signal to noise ratio in the coda. Typ- : points. Spatially resolved measurements are facilitated ically, epicentral distances are less than 110 km and the v by dense arrays of seismometers that have been set up coda lastsmore than30 secondsafter the directarrivals. i X occasionally. We note that the phase of coda waves has To perform the statistical analysis, we filtered the sig- r not been given much attention so far. The advantage of nal in a narrow frequency band centered around 7 Hz a phase is that it is not affected by the earthquake magni- ( 5%)andselecteda15stimewindowstartingaround5 ± tude,andthatitcontainspureinformationonscattering, secondsafterthedirectarrivals. Inthistimewindow,the notblurredbyabsorptioneffects. Forthestatisticalanal- signalis believed to be dominatedby multiple scattering ysisofamplitudeandphasefluctuationsofdirectarrivals, and is highly coherent along the array [6]. We evaluate we refer the reader to e.g. Zheng and Wu [4]. the Hilbert transformof the verticaldisplacement which yields the imaginary part of the complex analytic signal ψ(t,r)=A(t,r)eiφ(t,r). Fromthecomplexfield,twodefi- I. PHASE DISTRIBUTIONS nitionsofthephasecanbegiven: (1)Thewrappedphase φisdefinedastheargumentofthecomplexfieldψ inthe We study data sets from a temporary experiment de- range( π :π]. (2)The unwrappedphaseφ is obtained u − ployed at Pinyon Flat Observatory (PFO), California, by correcting for the 2π jumps – occurring when φ goes in 1990 by an IRIS program. This site exhibits a high through the value π - to obtain a continuous function level of regional seismic activity. The array (see Fig. 1) with values in R. ±The φ distribution is flat [7]. How- contained 58 3-components L-22 sensors (2Hz) and was ever, more information can be extracted by considering 2 higher-order statistics of the phase. For this purpose we 101 considerthespatialderivativeofthephase,whichcanbe estimated in two different ways. 100 (1) The first measurement relies on the difference of ) the wrappedphases ∆φ betweentwoseismometerssepa- φ ∆ ratebyadistanceδ. Applyingthesimplefinitedifference (10−1 P ′ formula φ ∆φ/δ an estimate of the spatial derivative ≃ ≈ isobtained. Notethatthephasedifference∆φtakesval- δ )/10−2 uesbetween 2πand+2πwhichdoesnotallowaprecise ′φ estimate for−the distribution of the derivative for values P( P(∆φ)measured P(∆φ )measured roughly larger than π/δ. Beyond this value our mea- 10−3 P(φ′)/uδfit surementswillbe dominatedbyfinitedifferenceartifacts P(∆φ)fit and the distributionis biasedby the 2π jumps occurring 10−4 P(∆φu)fit within the distance δ. 10−2 10−1φ′δ ∆φ 100 101 (2)Thesecondmethodusesthedifferenceofthephases ≃ φ spatially unwrapped at each time step. This yields u ′ FIG. 2: Distribution of the first derivative of the phase nor- another estimate of the derivative: φ ∆φ /δ which is ≈ u malized by the inter-station distance δ and measured using expected to suppress finite difference artifact. In prac- finite difference formula. (◦): wrapped phase; (+): un- tice it is impossible to discriminate a rare but physical wrapped phase. The colour lines represent the fits with large phase jump from a small fluctuation that causes a Gaussian theory. Green: wrapped phase; blue: unwrapped 2π jump just within the range δ. The only possibility phase; red : phase derivative. The fitting parameters are along1-Darraysis to impose that the largestadmissible Q=2.774 10−4 m−2 and g(δ)=0.993204. phase difference between two stations be smaller than π. ′ Hencethisφ estimatetakesvaluesin( π/δ :π/δ]andis − biasedclose to π/δ by the unwrappingprocessing errors. Inthelimitδ 0,thetwodefinitionsoughttobeequiva- N ltewnot bsteactaiuonsestt→heendpsrotboab0.ilitByyofapvehraasgeinjugmopvserbethtwee8ensetihse- P(ψ1···ψN)= πN d1etCexp−Xi,j ψi∗C−ij1ψj, (1) mic records, the lag-time in the coda, the east-west and   north-southorientations,andthesensorpositionswithin whereC = ψ ψ∗ isthecovariancematrix[8]. Itiscon- the array’sgridatfixedδ =7m,we calculatethe twore- venient tijo ushe niorjmialized fields so that C = 1. Then, ii sultingphasederivativedistributionswhichareshownto theoff-diagonalelementsareequaltothefieldcorrelation be non-uniform in Figure 2. It is also instructive to con- functionC =C( r r ). Fromthe jointdistribution ij i j | − | siderthesecond(third)derivativesofthephasewhichare of two fields, the probability distribution of the phase governed by the 3 (4)-point statistics which are plotted difference at two points located δ apart can be obtained inFigure3. Higher-orderderivativesareobtainedbyap- by integrating over the amplitudes [9]: plyingstandardfinitedifferenceformulastothewrapped 1 g2 F cos−1( F) phase(thischoiceisexplainedinthenextsection). Since P(∆φ)=P˜ − 1+ − (2) thethreefirstderivativedistributionsareevenfunctions, (cid:20)1 F2(cid:21)(cid:20) √1 F2 (cid:21) − − we only represent the positive values. They have all where F =gcos(∆φ) and g = ψ(r δ/2)ψ (r+δ/2) ; similar properties. For small values of the random vari- h − ∗ i P˜ = (2π ∆φ)/4π2 if ∆φ is the difference of the ables,the distributions arenearly flat. For largervalues, − | | wrapped phase and P˜ = 1/2π if ∆φ is the difference the distributions are governed by a power-law decay ex- of the unwrapped phase ∆φ . In the limit δ of to- cept for some peaks which stem from the finite distance tally uncorrelatedfields P(∆uφ)=P˜. In the l→imi∞t δ 0, between the seismic stations. In the following section, → wegetthefollowingformulaforthephasederivative[10]: we will demonstrate that the transition between the two behaviours is governed by the wavelength and that the 1 Q ′ P(φ)= . (3) power-lawdecayisaveryaccuratesignatureoftheGaus- 2[Q+φ′2]3/2 sian nature of the vertical displacements. Seismic coda is believed to be composed of multi- where Q = C′′(0) = 2(1 g)/δ2 for δ 0. Figure − − → ply scattered waves. Upon scattering, the many partial 2 demonstrates that the phase difference between two waves - associated with different paths in the medium nearby seismometers follows the prediction of Gaussian - would achieve random and independent phase shifts. circular statistics with excellent accuracy: we observe a The Gaussian nature would then follow from the cen- good agreement between the theoretical distribution of ′ tral limit theorem. We will thus assert that the coda φ (Eq. 3)andthemeasurementsover3ordersofmagni- waves obey circular Gaussianstatistics (CGS) according tude in probability and 2 orders of magnitude in deriva- towhichthejointprobabilityofNcomplexfielddisplace- tive. A clear discrepancy occurs for large values - typi- ′ ments ψ , recorded at positions r , is written as cally when φ π/δ - which can be perfectly explained i i ≈ 3 by the finite distance between the seismic stations. This is demonstrated in the same plot which shows excellent agreement between formula (2) and the measured finite- 100 difference statistics of the phase. We observe that the formula for the derivative (3) agrees with the observa- n tions on a larger range when the estimate of φ′ is based /δ ) on the wrapped phase difference. As a consequence, we )n prefer this method to find the higher derivatives. ((φ10−2 P In the frequency band of interest, there is experimental evidence that the vertical component of the coda is dominated by scattered Rayleigh surface waves [6]. As P(φ′′)measured P(φ′′′)measured a consequence, we expect the correlation function of the 10−4 P(φ′′)fit field to be given by [11] P(φ′′′)fit 10−2 10−1 100 101 ∗ φ(n)δn C(r)= ψ(0)ψ (r) =J (kr)exp( r/2ℓ), (4) 0 h i − FIG.3: Comparison oftheGaussiantheoryandobservations which agrees well with observations [7]. Equation 4 con- forthephasederivativesdistribution P(φ(n)) wheren=1or tains two length scales: the wavelength 2π/k and the th 2 denotes the n derivative of φ with respect to the spatial scattering mean free path ℓ. We emphasize that formula coordinate. From the fit we find: Q=2.774 10−4 m−2, R= (4) is valid even when the medium is slightly inelastic. 2.75 10−2 m−2 and S =4.0 10−2 m−2. Weak absorption does not enter the correlation function obtained from the long-time-tail of a diffuse wavefield [12]. Theformofthecorrelationfunction(4)impliesthat Q k2/2,ifkℓ 1. UsingtheparameterQobtainedby ≃ ≫ fitting the data with equation (3), we infer a dominant wavelengthλoftheorderof267mwhichagreeswiththe Cφ′(r) φ′(0)φ′(r) =(k/πr)exp( r/ℓ) (5) ≡h i − one predicted on the basis of the vertical profile of the elastic constants below the PFO array [13]. Q offers an for r > λ. Formula (5) has one crucial property. Con- accurate way of estimating the wavelength, alternative trary to the field correlation function C, Cφ′ does not to SPAC measurements [11]. The use of a narrow band oscillate on the wavelength scale but decreases with the signaliscrucialbecausetheparametersg andQstrongly mean free path as the sole characteristic length scale. depend on frequency. Any determination of the mean free path based on for- From the joint Gaussian distribution of 3 and 4 fields, mula (4) is impossible because the exponential decay is we canderiveanalytically the joint probabilityfunctions masked by the rapid cyclic oscillations. P(φ′,φ′′), P(φ′,φ′′,φ′′′) , featuring two new constants R As above,weestimate the phasederivativecorrelation and S [9]. From these formulas, the marginal distribu- function in a finite difference approximation using the tpiroonbsaPbi(liφt′y′),diPst(rφib′′u′)ticoannsboef tehvealufiarstet,dsneucomnedr,icaanlldy.thTihrde fro′′rm)/uδl2a. CCφ′o(nr)tra≃ryh(t∆oφtuh(er′)p∆roφbua(bri′l′i)tiy/δd2is=triCbu∆tφioun(|,r′w−e | derivativesofthephaseexhibitanasymptoticpower-law found out that the unwrapped phase difference offers a decay with exponents 3, 2, 5/3,respectively. These better estimate for the phase derivative. The wrapped − − − universal exponents (i.e. independent of the medium phase difference correlation function is too much domi- properties) can be obtained analytically and provide a nated by large 2π jumps with no physical interest. sensitive fit-independent test of the Gaussian nature of TheunwrappedphasedifferencecorrelationC (r)is ∆φu the seismic coda. This is illustrated in Figure 3. The measured along the two orthogonal arms with an aper- agreement leaves no doubt that coda waves are in the ture of 252 m and δ = 21 m interstation distance. The multiple scattering regime. data are averaged over orientation, lag-time in the coda andseismicevents. Theresultispresentedinlinearscale in the main part of Figure 4 and shows a decay domi- II. PHASE DIFFERENCE CORRELATIONS nated by the 1/r factor along the arms of the array, as predicted by formula (5). This supports our 2D picture We have shown that the distributions of the phase of wave diffusion. Due to the finite difference Cφ′(r =0) derivatives provide accurate information on the short- achievesafinitevalueatr =0,whichwefoundtobecon- range correlation properties of the field. In the follow- sistent with the variance of the unwrapped phase differ- ing we use the phase difference correlations to put some ence calculated from h∆φ2ui= −ππd∆φu(∆φu)2P(∆φu) constraints on the degree of heterogeneity which is re- and g(δ) = 0.98. The parameRter g(δ) has been deter- sponsibleforlong-rangecorrelationsalongthearray. For minedindependently byfitting the observeddistribution GaussianstatisticsandforsurfacewavesobeyingEq.(4), of P(∆φ ) with formula (2) for δ =21 m. u the phase derivative correlationfunction is [12]: Different reasons exist why it is more difficult to mea- 4 0.2 2 geometry. The results for different values of ℓ at fixed k are shown in Figure 4 together with the experimental ) φui1 results. The simulated function log(r×C∆φu) exhibits φui0.15 φ∆u 0 ainlEinqe.a(r5d)e(csaeye iwnistehtainsFloipguer−e14/).ℓTashewraefsosreee,ncofrorrelδat→ion0s ∆ ∆ φu rh of finite-size phase difference still offer direct access to ∆ n(−1 the mean free path of the waves in the crust. It is quite h0.1 l = encouraging to see (inset Figure 4) that the asymptotic r) −20 100 200 300 400 500 exponential regime is already reached for r > λ/5. We ( φ r (m) find that the 1/ℓ slope could in principle be measured ∆ C0.05 if the apertur−e of the network would have been a few wavelengths in size (>500 m), which in general is much smaller than the mean free path. Unfortunately, the 0 experimental error bars of our data set are too large to 0 50 100 150 200 250 r (m) permit accurate estimates of ℓ at PFO although it can roughly be bounded between 1 km and 10 km, much FIG.4: Unwrappedphasedifferencecorrelationfunction. (◦): smaller than the typical path (40 km) taken by coda experimental results. Dashed line: 1/r fit. Inset: logarithm waves coda during 20 s. This has to be interpreted as a ofthecorrelation multipliedbyr fordata(◦witherrorbars) rough estimate for the mean free path of the Rayleigh and numerical simulations at fixed k and ℓ=10 km (green), waves. ℓ=1 km (red), ℓ=500 m (red) and ℓ=200 m (yellow). In conclusion, from their observed first three spatial derivatives of phase, seismic coda waves are proved to surethecorrelationfunctionofthephasederivativethan obey Gaussian statistics with high accuracy, with a lo- e.g. theprobabilitydistributions. First,4th-orderstatis- cal correlation on the scale of the wavelength. At longer tics require more averaging to suppress unwanted fluc- scales, we demonstrate that the correlation function of tuations in the data. Secondly, the interstation dis- the spatial derivative of phase offers a new, promising tance δ = 21 m along the arms reduces the correlation opportunity to measure directly the mean free path ℓ of between the fields at two nearby stations significantly, seismic waves. This measure is not affected by absorp- which favours systematic errors in the derivative. Fi- tion neither to scattering anisotropy. As opposed to the nally, due to frequent breakdowns of the sensors located coherentbackscatteringeffect [14], the proposedmethod near the ends of the arms, the data could not be aver- does not require the sensors to be located in the near aged over all sensor positions and all seismic events. As field of the source, but requires an array of seismic sta- a consequence the correlations for distances r > 180 m tions with sub-wavelength spacing between stations and had to be excluded. with a total aperture of a few wavelengths. Our work Since formula 5 is valid only in the limit δ 0, we → highlights the phase as a useful physical object to study have evaluated the impact of the finite inter-station seismic coda waves. distance inherent to our experimental set-up. To this end, we have simulated N correlated Gaussian random WethankJ.Page,M.Campillo,P.RouxandE.Larose field displacements on a virtual array with the PFO for useful discussions. [1] K.Aki, J. Geophys. Res. 74, 615 (1969). [8] J. W. Goodman, Statistical optics (New York: Wiley, [2] K.AkiandB.Chouet,J.Geophys.Res.80,3322(1975). 1985, 1985). [3] H. Sato and M. Fehler, Seismic wave propaga- [9] M. L. Cowan, D. Anache-Ménier, W. K. Hildebrand, tion and scattering in the heterogeneous Earth (AIP J. H. Page, and B. A. van Tiggelen, Phys. Rev. Lett. Press/Springer Verlag, New York,1998). 99, 094301 (2007). [4] Y. Zheng and R.-S. Wu, Geophys. Res. Lett. 32, 14314 [10] A. Z. Genack, P. Sebbah, M. Stoytchev, and B. A. van (2005). Tiggelen, Phys. Rev.Lett. 82, 715 (1999). [5] T.J.Owens,P.N.Owens,P.N.Anderson,andD.E.Mc- [11] K. Aki,Bull. Earthq. Res. Inst. 35, 415 (1957). Namara, Tech. Rep., Incorporated Research Institution [12] B. A. van Tiggelen, D. Anache, and A. Ghysels, Euro- in Seismology (1991). physics Letters 74, 999 (2006). [6] L.Margerin,M.Campillo,B.VanTiggelen,andR.Hen- [13] F. L. Vernon, G. L. Pavlis, T. J. Owens, D. E. Macna- nino, ArXive-prints803 (2008), 0803.1114. mara, and P. N. Anderson, Bull. Seism. Soc. Am. 88, [7] See epaps document no. [ ] for formulae. this document 1548 (1998). can be reached through a direct link in the online arti- [14] E. Larose, L. Margerin, B. A. van Tiggelen, and cle’s html reference section or via the epaps homepage M. Campillo, Phys.Rev.Lett. 93, 048501 (2004). (http://www.aip.org/pubservs/epaps.html).