Long-Range Correlations in the Diffuse Seismic Coda

See allHide authors and affiliations

Science  24 Jan 2003:
Vol. 299, Issue 5606, pp. 547-549
DOI: 10.1126/science.1078551


The late seismic coda may contain coherent information about the elastic response of Earth. We computed the correlations of the seismic codas of 101 distant earthquakes recorded at stations that were tens of kilometers apart. By stacking cross-correlation functions of codas, we found a low-frequency coherent part in the diffuse field. The extracted pulses have the polarization characteristics and group velocities expected for Rayleigh and Love waves. The set of cross-correlations has the symmetries of the surface-wave part of the Green tensor. This seismological example shows that diffuse waves produced by distant sources are sufficient to retrieve direct waves between two perfectly located points of observation. Because it relies on general properties of diffuse waves, this result has potential applications in other fields.

Seismologists have used coherent seismic waves to image the structure of Earth's interior. Velocity variations of seismic waves with depth can be derived from arrival times (with the use of ray theory) or from the dispersion properties of coherent surface waves. One difficulty with these imaging techniques is that they require energetic sources such as large explosions or earthquakes that can be located with a specified accuracy. Any approach that would help to evaluate the response of Earth to a perfectly known point source—that is, the Green function—would be most welcome.

This problem of imaging in strongly diffractive media is also a challenge for acoustics or optics. It was recently demonstrated in laboratory experiments with ultrasonic and thermal noise that the Green function can be measured from the correlation properties of diffuse fields (1, 2). Here, we show that the use of field-to-field correlation to retrieve the Green function is a valid approach not only in the extremely controlled and favorable conditions of the laboratory but also with natural signals such as seismograms produced by earthquakes. In seismology, it has been recognized that coda waves, which make up the late part of seismic signals (Fig. 1), are the result of scattering from small-scale heterogeneities in the lithosphere (3–5). The physics of coda waves cannot be fully understood with classical ray theory. Multiple scattering plays a prominent part in the seismic coda, and seismologists have made use of the radiative transfer theory to model the coda intensity (6–9). Recently, the diffusive character of the coda was revealed (10) by investigating the property of mode equipartition (11). This phenomenon is a property of diffuse elastic waves and shows up as a stabilization of the ratio of S- and P-wave energies in time, independent of the source. Radiative transfer and diffusion are concepts that apply only to the evolution of the average energy of waves in random elastic media (12–14). They disregard the phase of the diffuse field despite experimental evidence of the importance of phase information in optics (15) and acoustics (16).

Figure 1

Location map of the broadband stations CUIG, YAIG, and PLIG of the Mexican National Seismological Network (black squares) and epicenters of 30 earthquakes of the data set (white circles). Inset: An example of a record of one of these events at station PLIG (vertical component).

Here, we use the coherence of diffuse waves to retrieve direct waves between two points at Earth's surface. In this approach, the cross-correlation function between the wave fields produced by a single source at two points is averaged over the source location. Assuming a modal representation of the wave field, this spatially averaged correlation is an approximation of the Green function between the two points of observation (1, 2, 17,18). However, because we cannot expect a homogeneous distribution of earthquakes, we have to rely on the distribution of scatterers responsible for the diffusion to perform a sufficient averaging. An alternative argument can be found in the property of modal equipartition of the diffuse field. Equipartition occurs because multiple scattering tends to homogenize the phase space. For direct arrivals, the energy is distributed in the phase space in a manner that depends on the nature and position of the source. In contrast, energy becomes uniform in phase space when entering the diffusive regime. This property is independent of the details of the heterogeneities that produce the scattering. Considering a time window delayed enough from the first arrival for the waves to have become diffuse as a result of multiple scattering, we can write the displacement u at location R and time t in the form of its expansion in the eigenfunctionsΦ n of the elastic medium:Embedded Imagewhere Ωn are the eigenfrequencies and ɛn are statistically identical independent random variables. The expression of the cross-correlation of the displacements at two different locations is, on average, close to the Green function between these two locations.

We applied this approach to a seismic data set from central Mexico that fits several basic criteria: (i) The region is seismically active, with numerous earthquakes of magnitude large enough to excite late coda (M > 4.5), (ii) good-quality broadband records of these events are available, (iii) the Green function is already known and displays features striking enough to be recognized easily in a noisy time series, and (iv) high-frequency coda waves exhibit a diffuse behavior in this region (10), a property that we expect to be verified in a broad frequency range.

The stations PLIG and YAIG were selected because of the availability of good-quality records of 101 regional events (Fig. 1). The horizontal components of the seismograms were rotated assuming the interstation great circle path to be the radial direction. No filter was applied to the broadband seismograms. For most records, the signal-to-noise ratio was good enough to process the seismograms over coda windows of a few hundred seconds (Fig. 1, inset). As a result of the exponential decrease of coda amplitude with time, a simple cross-correlation between the coda signals recorded at the two stations would strongly overweight the earliest part of the coda. To compensate for amplitude attenuation with time, we divided the coda windows into 100-s-long segments with an overlap of 25 s. Then we computed the cross-correlations between these truncated signals and normalized the amplitudes of each correlation to unit maximum. The resulting 595 normalized correlations were then stacked to give the signals (Fig. 2A). Similar results are obtained by disregarding the amplitude completely and considering one-bit signals (19, 20).

Figure 2

(A) Stacked correlation functions computed from vertical-component records at station PLIG and vertical (Z), radial (R), and transverse (T) component records at station YAIG. Because the individual cross-correlations were normalized before the stack, the amplitudes are arbitrary. (B) Particle motion plots in the time window of the pulse (15 to 35 s).

The resulting average cross-correlations computed from vertical-component records (Z) at station PLIG and vertical, radial, and transverse component records (respectively Z, R, and T) at station YAIG (Fig. 2A) show a clear 8-s-period pulse between 20 and 30 s in both the vertical and radial traces, whereas no coherent signal is visible on the transverse component. During the stacking process, the amplitude of the pulse increases linearly with the summation order, showing that the pulse results from the stack of coherent signals present in most of the individual correlation traces. In contrast, the average amplitude of the noise varies as the square root of the number of individual windows, as expected for a summation of incoherent signals. The coherent signal is not visible in individual cross-correlations because the signal-to-noise ratio is only ∼0.2 for the case of a 100-s window. The particle motion (Fig. 2B) is restricted to the propagation plane with elliptic polarization typical of a Rayleigh wave pulse, a disturbance propagating at the free surface of an elastic body.

We computed the average correlations between all components of the ground motion at the two stations (Fig. 3A). The theoretical Green tensor has been computed in a three-layer crustal model (21) and convolved with an 8-s-wide Ricker wavelet. It displays only a few distinct features. It is dominated by a strong Rayleigh pulse at 25 s for both the vertical and radial point-force sources, and a Love pulse—a horizontally polarized guided shear wave—for the transverse point-force source. For the radial force, the Rayleigh pulse on the radial component (R/R) is preceded by a strong-amplitude body wave. The tensor obtained by stacking correlations between coda records of regional earthquakes at stations PLIG and YAIG displays the same symmetries as the theoretical Green tensor (Fig. 3). Moreover, the arrival times of the pulses in the Z/Z, Z/R, R/Z, R/R, and T/T components of the stacked correlations coincide with those of the Rayleigh and Love signals in the Green tensor. This coincidence in arrival time, as well as the clear Rayleigh and Love polarization of the correlation pulses (Fig. 2B), proves that the observed signals are identified as the Rayleigh and Love pulses of the Green tensor and, most important, that the coda correlation technique does indeed retrieve the surface-wave part of the actual Green tensor between the two stations.

Figure 3

Comparison between the nine stacks of correlation traces at stations PLIG and YAIG (A) and the nine components of the theoretical Green tensor (B) computed for a 69-km source-station distance. The 1-D average shear wave velocity model used here was measured for the crust of Central Mexico from inversion of group velocity dispersion curves estimated for paths between the Guerrero-Michoacàn subduction zone and Mexico City (20).

To make sure that the pulse is not simply a surface wave that is generated repeatedly at the coast by the conversion of oceanic waves and that propagates in the direction defined by the two stations, we performed the same test with another pair of stations, YAIG and CUIG, oriented in a different azimuth (Fig. 1). The stacked correlation signals also display pulses with arrival times and polarizations close to the Rayleigh and Love modes of the theoretical Green function (fig. S1), excluding the alternative interpretation of induced surface waves.

So far we have not been able to extract either the high-frequency part of the Green function or the body waves. The lack of high frequencies is most probably a result of the absence of high-frequency waves in the late coda because of anelastic absorption, which acts as a low-pass filter. Another explanation could be that the fundamental modes of Rayleigh and Love waves at low frequency are the part of the field with the simplest modal representation. Retrieving the Green function relies on the orthogonality of the set of eigenfunctions that constitutes the total random field. All cross-products vanish in the averaging, assuming a distribution of sources, or scatterers, that spans the whole space. However, the volume where the spatial source averaging is performed is in practice limited by the number of earthquakes and the locations of scatterers. We speculate that only eigenfunctions with amplitudes concentrated in a zone where inhomogeneities are densely distributed can be adequately extracted. This is the case with the Rayleigh and Love waves, the eigenfunctions of which have a limited penetration in the upper part of the crust where the distribution of scatterers is likely to be dense.

We expect to retrieve both the Green function and its time reciprocal if the diffuse field is perfectly random. This could be the case with an isotropic distribution of sources around the stations or in a finite body. Because all earthquakes are located south of both station PLIG and station YAIG, there is a preferential direction of transport of diffuse energy. This results in a better reconstruction of the Green function in one of the time directions. We also considered a couple of stations along the coast (fig. S2A) for which the distribution of epicenters is more symmetric. The wave propagation is much more complex there (22) than in central Mexico, but some features of the Green function emerge from the noisy correlation stacks, such as a clear dispersed Love wave that can be seen in the two directions of time (fig. S2B).

Digital seismic networks provide a large number of coda records, which can be used to compute impulse response between perfectly located positions. This new kind of seismogram could help to produce images of the inner Earth structures without the uncertainties of origin time and source location encountered with traditional earthquake data. A similar approach is applicable in other domains where time series of diffuse waves are available.

Supporting Online Material

Figs. S1 and S2

  • * To whom correspondence should be addressed. E-mail: Michel.Campillo{at}


View Abstract

Navigate This Article