Observation of Equipartition of Seismic Waves

R. Hennino,1 N. Trégourès,2 N. M. Shapiro,3 L. Margerin,1,4 M. Campillo,1 B. A. van Tiggelen,2,* and R. L. Weaver5 1Laboratoire de Géophysique Interne et Tectonophysique, Observatoire de Grenoble, Université Joseph Fourier, B.P. 53, 38041 Grenoble Cedex, France 2Laboratoire de Physique et Modélisation des Milieux Condensés, CNRS/Université Joseph Fourier, B.P. 166, 38042 Grenoble Cedex 09, France 3Instituto de Geofisica, Universidad Nacional Autonoma de México, CU 04510 México, D.F., Mexico 4Department of Geosciences, Princeton University, Guyot Hall, Princeton, New Jersey 08544 5Department of Theoretical and Applied Mechanics, University of Illinois at Urbana–Champaign, Urbana, Illinois 61801 (Received 29 September 2000)

Scattering of classical waves in controlled laboratory experiments has extended our knowledge of wave propagation considerably [1][2][3].One new challenge is to apply this knowledge to waves in natural environments, whose properties are complex and not well known, for example, medical imaging [4], remote sensing [5], and seismology.
The seismic case is particularly difficult because several specific complications, such as near-field detection and mode conversions at the Earth's surface, rule out the application of methods developed in optics and acoustics.In addition, seismic records are intrinsically time resolved.Sheng et al. [6] considered the possibility of Anderson localization in time-dependent signals that have propagated through a layered Earth crust.
The seismic coda refers to the pronounced exponential time tail observed in the seismograms of regional earthquakes in the frequency band 1-10 Hz, reaching the level of seismic noise after sometimes more than 10 times the travel time of direct waves [7].One important observational fact is that the time decay coefficient of the energy in the seismic coda, the so-called coda Q factor, is a regional constant [8], dependent on frequency, but independent of distance, depth, and magnitude of the seismic source.
Many studies have made an attempt to relate physical properties of the Earth's lithosphere to the observed regional seismic coda.The first pioneering studies by Aki and Chouet [7,8], followed by others [9,10], interpreted seismic coda as elastic waves that are singly scattered from inhomogeneities in the Earth's crust.It was shown that a combination of single scattering and dissipation can account for some aspects of the seismic coda [11].Recent numerical studies suggested multiple scattering as the origin of seismic coda [12][13][14][15][16].Our own study estimated mean free paths ᐉ ഠ 20 70 km in the frequency band 1-10 Hz [16], comparable to the thickness H of the Earth's crust.For a shear wave velocity b ഠ 3.5 km͞s in the crust this implies a scattering mean free time of only 15 s, much smaller than the observed duration (180 s) of the coda.
A direct and model-independent confirmation of the pertinence of multiple scattering is important for two reasons.First, it is important for future studies in extracting information correctly from the seismic coda, i.e., to do the inverse problem.Second, it would facilitate the application of many powerful mesoscopic tools developed in acoustics and optics [2,3].However, several complications arise.One "elementary" test for multiple scattering is the randomization of the propagation direction of the waves [2].Unfortunately, seismic observations measure local elastic displacements u͑r, t͒, and not "specific intensities," i.e., angularly resolved fluxes, as in optics.Another elementary test would be an observation of the randomization of polarization [5] as indicating multiple scattering, since single scattering should be highly polarized for a polarized source such as an earthquake.The essential problem in seismic studies is that transverse S and longitudinal P waves cannot be separated in a ground displacement at one point.Only recently, Shapiro et al. [17] presented an approach to separate P from S waves using a small-aperture array.A second complication is the presence of the Earth's free surface that coherently mixes transverse SV waves with longitudinal P waves [18].
In this Letter, observational data of seismic coda are confronted with a first principle in multiple scattering to resolve the controversy "single or multiple scattering in the seismic coda."We have investigated the principle of equipartition, first discussed by Weaver [19,20] for elastic waves in a solid body, more recently by Papanicolaou et al. [21] and Turner [22] using a transport equation, and finally for light in liquid crystals [23].It is based on the fact that multiple scattering tends to homogenize phase space: the spectral energy density around a specified frequency, originally distributed in phase space in a fashion that largely depends on the nature of the source, eventually becomes uniform.Equipartition has the remarkable property to be independent of the unknown fluctuations that cause the scattering.It would imply any energy ratio to become time independent, and independent of the magnitude, polarization, and distance of the seismic source.In single scattering from a polarized source these energy ratios would vary significantly with time and among different sources.
The total elastic wave energy density is given by [24] with u͑r, t͒ the local displacement vector, and r the mass density.Deep-reflection and surface-wave studies have shown that the Lamé coefficients l and m of the strainstress tensor are roughly constant and equal in the Earth's crust, with fluctuations less than 5%.The different terms in Eq. ( 1) represent kinetic (K), compressional (P), and shear (S) energy density.The last term (I) is an interference term involving cross terms of the kind ≠ i u j ≠ k u l whose averages vanish except near boundaries.We have chosen to study S͞P, K͑͞S 1 P͒, and I͑͞S 1 P͒, and the ratio H 2 ͞V 2 of kinetic energies for elastic displacements in different directions for which no spatial derivatives have to be carried out.
The most direct way to facilitate the evaluation of the three partial derivatives would be to install seismic receivers close to each other and at different depths.However, installing numerous seismic receivers, especially at depth, is extremely expensive.An array of four receivers at the corners of a 50 m side square [17] was temporarily set up close to the city of Chilpancingo (Mexico).The high seismic rate in Mexico made it possible to record, during the three months of the experiment, a series of local earthquakes with magnitudes between 4 and 5, and a spread of roughly 300 km in epicentral distance.The sensors are CMG-40T seismometers connected to Reftek digitizers, whose absolute time was provided by GPS satellites.The three vertical derivatives ≠ z u i were deduced from the measured two horizontal partial derivatives by imposing the stress-free condition l͑div u͒d iz 1 m͑≠ i u z 1 ≠ z u i ͒ 0 for i x, y, z at the local Earth's surface.
We selected 12 earthquakes that exhibit a pronounced seismic coda in the frequency band 1 3 Hz.Figure 1 shows the observed seismograms of one selected event of magnitude 4.3 at a distance of 35 km.They have been averaged over time windows of Dt 15 s separated by 0.3 s.The first arrivals of S and P waves, as well as the exponential coda at long times, can easily be identified.In the regime of seismic coda all ratios fluctuate around a constant value until the signal-to-noise level is reached.The horizontal lines in Figs.1c and 1d locate the time average.We emphasize that in this regime, the local elastic energy density decays by 4 orders in magnitude.
We interpret the observed time-independent energy ratios in the coda as a sign of equipartition.The time- dependent fluctuations around the mean for one event are attributed to the Gaussian fluctuations that ought to be uncorrelated in time and of the order of 1͞ p DtDf͞2 for a bandwidth Df.In spite of the large variations of the sources in seismic magnitude and distance, Fig. 2 shows that the measured energy ratios are roughly the same for all sources.Especially the uncertainties in the ratios of potential energies (Table I) are small.The nonzero value of ͗I͘ indicates the presence of mode conversions at the Earth's surface.
We now attempt to understand the observed ratios by assuming that the coda is an equipartitioned random field.If u͑r, t͒ is the local, time-dependent displacement vector, it can be expanded into the eigenfunctions u n of the elastic medium with eigenfrequencies v n , u͑r, t͒ X n ´ne 2iv n t u n ͑r͒ . ( In the presence of disorder, all modes get mixed and the ´n become time-dependent random variables.For small disorder the eigenfrequencies v n are not significantly altered.The simplest way of expressing equipartition is that the average ͗´n͘ 0 and that ͗´n´‫ء‬ m ͘ s 2 ͑t͒d nm .This means that the mode amplitudes are independent random variables with random phase but with equal variance [25].This statement implies an "equipartition" of the elastic energy (1) among the different modes.The time dependence s 2 ͑t͒ cancels in any energy ratio.
In a homogeneous unbounded medium the modes are plane waves with either S or P polarization.In that case ͗I͘ 0, ͗K͗͘͞S 1 P͘ 1, and ͗S͗͘͞P͘ 2͑a͞b͒ 3 [19,21,26], in terms of the P velocity a p ͑l 1 2m͒͞r and the S velocity b p m͞r, including a factor of 2 to acknowledge the twofold degeneration of the transverse S waves.This ratio is recognized as the ratio of the density of states (DOS) for both modes.For l m we get ͗S͗͘͞P͘ 10.39.The observed energy ratios do not agree with these quantitative predictions.
A more realistic model should take into account the mode conversions at the Earth's free surface, and at the interface between crust and mantle (the Moho) at roughly 30 km depth.The equipartition process within a distance of a few shear wavelengths (l s ഠ 2 km) from the surface is dominated by the first.We have calculated the exact elastic eigenmodes of a homogeneous plate, bounded by two free surfaces [25,27].For a thickness H . 10l S , the equipartition values near the surface, calculated from Eq. ( 2), are independent of H, with surface values that agreed with previous work for the elastic half-space [25].Figure 3 shows how the ratios ͗S͗͘͞P͘ and ͗K͗͘͞P 1 S͘ depend on depth if equipartition between all modes is assumed.The oscillations in Fig. 3, associated with the local density of states (LDOS) near the boundary, are well known in optics [28,29].Our present study has access to the surface values only (Table I).The agreement between model and observations is remarkable.Table I also shows the theoretical values when either only surface waves or only bulk waves are considered.In both cases they fall outside the estimated error bar of the observational values.This leads to a picture of the coda consisting of an equipartitioned mixture of both Rayleigh and bulk waves, which is expected when the absorption time of Rayleigh waves exceeds the time to mode convert into bulk waves.
The observed value H 2 ͞V 2 ഠ 2.5 agrees less well (Table I).One possible source of discrepancy could be the presence of a "true" effective surface different from both the local geographical free surface and the local surface of the Earth, due to fluctuations in the slope of the Earth's surface.To investigate this option, the full kinetic covariance tensor ͗ 1 2 r≠ t u i ≠ t u j ͗͘͞K͘ was calculated from the data.Its off-diagonal elements are less than 5% 6 5%, i.e., consistent with zero and thus with our equipartition model.However, the observed eigenvalues ͗X 2 ͘:͗Y 2 ͘:͗Z 2 ͘ 0.27:0.45:0.28(Y is north) disagree with the theoretical diagonal elements ͗X 2 ͘:͗Y 2 ͘:͗Z 2 ͘ 0.32:0.32:0.36.No frame can be found that makes the entire covariance matrix consistent with our equipartition model.Some other source of discrepancy must exist.
Alternatively, the discrepancy could be attributed to excess absorption of surface waves.It is readily shown that if the absorption time of a surface wave is less than the time to mode convert to bulk waves, its contribution to the steady state partition will be diminished.Table I shows that a small deficit of Rayleigh waves could bring H 2 ͞V 2 into accord with theory without seriously compromising the others.In addition, an anisotropic absorption of Rayleigh waves might even explain the anomalous value of X 2 ͞Y 2 0.6.In view of the geographical location of the array, a mountain area stretched in the direction E-W, it is quite possible that the component Y is privileged over the X direction.
In conclusion, we have observed a first principle of multiple scattering -energy equipartition -in seismograms recorded in Mexico.This strongly supports the interpretation of seismic coda as a genuine multiple scattering process, and excludes single scattering as an alternative explanation.

FIG. 1 .
FIG. 1. Observed seismogram, bandpassed between 1 and 3 Hz, for event 11 at an epicentral distance of 35 km from the detection array and a magnitude of 4.3.(a) Linear plot of the bandpassed displacement measured as a function of time.(b) Semilogarithmic plot of the energy density.A distinction is made between kinetic energy (K, dashed line), shear energy (S), and compressional energy (P).(c) Linear plot of the energy ratio S͞P.(d) Linear plot of the energy ratios K͑͞S 1 P͒ (solid line) and H 2 ͞V 2 (dashed line).The horizontal lines denote the estimated time average.

FIG. 2 .
FIG. 2. Observed time-averaged energy ratios S͞P (top left), K͑͞S 1 P͒ (top right), H 2 ͞V 2 (bottom left), and I͑͞P 1 S͒ for all 12 events.The shadowed bar denotes the mean value with the standard deviation.The error bars denote the observed time-dependent fluctuations.

TABLE I .
Observed averaged energy ratios (in the frequency range 1-3 Hz, with standard deviation) compared to the theory near the free surface, with l m independent of depth.The third and fourth columns are calculated by considering both bulk and Rayleigh waves.The last columns are calculated by considering only surface Rayleigh waves or only bulk waves. .3. Theoretical prediction for the depth dependence of the average energy ratios S͞P (normalized to bulk value), K͑͞P 1 S͒ and I͑͞P 1 S͒, assuming that l m. FIG