Polarization of the Earth’s teleseismic wavefield

Summary

The average polarization of the Earth's teleseismic wavefield is imaged by applying a vector stacking technique to ∼32 000 three-component broad-band seismograms. Polarization record sections capture the average radial-vertical particle motion in two narrow frequency bands: one above the microseism peak (0.5–2.0 Hz) and one below (0.03–0.07 Hz). Within these record sections, colour indicates polarization angle and brightness denotes particle motion linearity. These parameters are estimated using the time-domain method of Jurkevics (1988), enhanced to allow the simultaneous analysis of multiple sources. The low- and high-frequency record sections highlight different features in the teleseismic wavefield. The low-frequency stack shows P, S and their surface multiples as well as more exotic phases originating from combinations of reflections, refractions and conversions. The high-frequency stack detects fewer arrivals. However, the attenuation of the multiple S waves allows the detection of later-arriving core phases. Both the high- and low-frequency stacks illuminate regions containing consistent waveforms and warn of contaminating arrivals. The coherence of lower-mantle-turning P waves is indicated by the agreement between the observed and predicted polarization angles. In contrast, a large discrepancy between the observed and predicted low-frequency S angles evidences strong SPL contamination between 20° and 60°. Both polarization stacking and conventional amplitude stacking techniques allow the detection of seismic phases that are buried within the noise on single seismograms, but polarization stacking additionally yields polarization angle and particle motion linearity.

body waves, polarization, traveltimes, whole Earth

1 Introduction

Fine details of deep earth structure lie buried within a few hundred thousand broad-band seismograms recorded by the Global Seismic Network (GSN). It is difficult to see the subtle signals contained in such an enormous waveform collection by examining the individual traces. Thus, images are produced from large subsets of this and other waveform databases to reveal details of Earth structure such as the fate of subducted slabs (Castle & Creager 1998), mantle layering (Revenaugh & Jordan 1989), deep reflectors (Kawakatsu & Niu 1994), discontinuity thickness (Benz & Vidale 1993), the extent of small-scale mantle heterogeneities (Hedlin 1997), long-wavelength discontinuity topography (Flanagan & Shearer 1998), the impedance contrast across the 660 km discontinuity (Estabrook & Kind 1996), and small-scale core–mantle boundary topography (Earle & Shearer 1997b).

The long-period (≥ 60 s) global seismic wavefield is captured in stacked record sections (Shearer 1994b), which are several hours long and cover the angular range from 0° to 180°, At periods close to 20 s, Shearer (1991) used an automatic gain control (AGC) algorithm to image the Earth’s main body wave arrivals. Astiz (1996) extended the AGC method to three-component seismograms for a range of frequency bands.

These global record sections clearly show the predominant body wave arrivals and are used in many ways. For example, they (1) identify robust features that can be scrutinized for Earth structure, (2) aid in earthquake detection and location by providing templates for matched filter processing (Shearer 1994a), (3) document phases not included in published theoretical traveltime curves, alerting seismologists to the existence of contaminating phase arrivals, and (4) are useful for education purposes.

Here I apply a vector stacking technique to GSN seismograms and capture the main polarization signature of the high- and low-frequency (0.5–2.0 and 0.03–0.07 Hz) global seismic wavefield in the radial–vertical plane. In addition to imaging the main body wave arrivals, this method determines the average particle motion linearity and direction. Knowledge of these polarization parameters allows the classification of arrival type (shear or compressional) and the determination of arrival angle.

Polarization stacking has several advantages over traditional stacking techniques such as automatic gain control and multiple-event coherent waveform stacking. Emergent arrivals such as SPL are sometimes not detected by AGC algorithms. These phases can be identified using polarization methods because polarization analysis relies on phase and amplitude differences between components, whereas AGC techniques are sensitive to rapid amplitude increases along a single trace. Additionally, the polarization angle and particle motion linearity can be used to determine the origin of scattered energy and characterize seismic coda. Multiple-event coherent waveform stacking techniques are less effective at high frequencies due to pulse distortion introduced by heterogeneities and source variations; however, a high-frequency phase arrival’s polarization is often consistent amongst events, allowing the application of polarization stacking.

2 Polarization of Teleseismic Arrivals

The ground motion of a teleseismic arrival can be quite complicated, being affected by seismic velocity heterogeneity, anisotropy, surface conversions, attenuation, source characteristics and wave type. However, for a limited time window and frequency band, it can be fitted using a simple elliptical polarization model.

When this model is applied to two-component data, the ground motion is fitted to an ellipse defined by a polarization angle (θ) and linearity (ρ) (Fig. 1). The definitions of θand ρ are study-dependent. Here θ is taken to be the angle from the vertical (positive away from the source) to the major axis of the best-fitting ellipse and ρ=1 − (a/b)ⁿ, where a and b are the lengths of the ellipse’s minor and major axes respectively and n is an empirically chosen parameter. The polarization angle θ ranges from − 90° to 90° and ρ varies from 0 for circularly polarized energy to 1 for linearly polarized energy. A proper choice of n ensures the values of ρ are well distributed between 0 and 1.

Figure 1

Determination of the polarization parameters. The dashed line represents the particle motion of a compressional arrival in the radial–vertical (R–V) plane. The polarization parameters θ and ρ are obtained from the best-fitting ellipse (solid line).

The polarization angle and linearity aid in arrival classification (e.g. P, S, Rayleigh or Love). Body waves display nearly linear polarization, whereas the surface wave particle motion is more elliptical. P and SV teleseismic arrivals can be distinguished because compressional wave particle motion is correlated in the vertical–radial plane, whereas shear wave motion is anti-correlated.

The P polarization angle approximates the arrival angle. Therefore, θ can aid in hypocentre location as well as arrival classification. Note that the arrival angle deviates slightly from θ because of surface mode conversions (Nuttli 1961; Aki & Richards 1980; Bokelmann 1995). The P polarization angle becomes closer to vertical as the source–receiver separation increases, while SV polarization approaches horizontal. Further details on seismic signal polarization can be found in seismological and geophysical time-sequence analysis texts (e.g. Aki & Richards 1980; Kanasewich 1981).

3 Data

The polarization stacks are generated from broad-band waveforms recorded by GSN stations from January 1988 to July 1997. These data were obtained from the online IRIS FARM archive, which contains seismograms from 2327 earthquakes for this time period. The FARM has a minimum-magnitude cut-off of 5.7 M_w for events shallower than 100 km and a 5.5 M_w cut-off for deeper events. To increase the signal-to-noise ratio and avoid complications resulting from depth phases, only the 973 earthquakes with magnitudes greater than 6.0 M_w and depths less than 50 km are used in this paper. Fig. 2maps the source–receiver midpoints that correspond to the 32 431 three-component seismograms used here. Note the good global coverage.

Figure 2

Data coverage. The pluses indicate source–receivermid-points for seismograms used in the polarization stacks.

4 Method

Both time- and frequency-domain algorithms have been developed to estimate a seismogram’s time-varying polarization signal. Time domain techniques generally utilize eigen-value decomposition of the particle motion covariance matrix within a moving time window (e.g. Flinn 1965; Kanasewich 1981; Phillips 1993). This method has been extended to analytic seismograms [a complex time-series defined as S(t)+iH(t), where S(t) is the seismogram and H(t) is its Hilbert transform] (Vidale 1986; Bataille & Chiu 1991; Morozov & Smithson 1996). Frequency-domain techniques produce polarization estimates as a function of time and frequency (Park 1987; Lilly & Park 1995). Relatively few studies (Jurkevics 1988; Bataille & Chiu 1991) have taken advantage of the signal enhancement afforded by stacking polarization estimates from multiple seismograms.

In this study the algorithm of Jurkevics (1988) is extended to allow the combination of seismic records from multiple earthquakes. As with previous time domain methods, the polarization parameters in a given time window are determined using eigenvalue decomposition of the particle motion covariance matrix Ψ. For the two-component vertical–radial case we have  

(1)

where  

(2)

where N is the number of points in the time window and the vertical (V) and radial (R) components are assumed to have a zero mean.

The eigenvalues and eigenvectors of Ψ define an ellipse that is the least-squares best fit to the particle motion. The orientation of the ellipse’s major and minor axes is defined by the eigenvectors, and the axes’ lengths are equal to the square root of the corresponding eigenvalues. Ψ is a positive semi-definite matrix, thus its eigenvalues are real and non-negative.

The polarization parameters shown in Fig. 1are obtained from the eigenvectors and eigenvalues. θ is the angle from the vertical to the eigenvector corresponding to the largest eigenvalue (λ₁) and ρ=1 − (λ₂/λ₁)^n/2. For these data, setting n=1.0 produces values for ρ that are well distributed between zero and one.

To combine the polarization estimates from multiple earthquakes, I average normalized covariance matrices for all radial and vertical trace pairs within a given time–distance window. The covariance matrices must be normalized before averaging to account for different magnitude earthquakes. They are normalized such that the sum of the major and minor axes of the best-fitting ellipse is equal to one. This ensures that the average polarization angle calculated from two ellipses with equal major- to minor-axis ratios will lie halfway between the polarization angles of the individual ellipses. The normalized covariance matrix, ~ Ψ, is defined as  

(3)

The procedure for data processing and stacking is as follows.

(1) The three-component seismograms are corrected for the differing frequency-dependent gains often present on the vertical, north–south and east–west components. The relative gains are obtained from the instrument response at the middle frequency of the pass band.

(2) The seismograms are filtered to a high- or low-frequency pass band (0.5–2.0 or 0.03–0.07 Hz). This increases the STN by avoiding the microseism peak, and the band-limited signals provide better polarization estimates.

(3) The seismograms are rotated into their radial (R) and vertical (V) components. This study documents polarization in the R–V plane, leaving the transverse component unused.

(4) The normalized covariance matrix ~ Ψ is calculated from the R–V particle motion in a time window of length ΔT centred on each sample point along the seismogram. The choice of ΔT depends on the predominant frequency of the seismogram. ΔT of 1.4 and 35 s were used for the high- and low-frequency band passes respectively.

(5) All ~ Ψ falling within a given time–distance bin are averaged and the polarization angle and rectilinearity are calculated from ~ Ψ. The bin sizes are 0.5°× 10 s for the high-frequency stack and 0.75°× 15 s for the low-frequency stack. The bin sizes were chosen by trial and error to maximize the benefit of noise reduction from averaging whilst providing the necessary resolution. The smaller bin size used for thehigh-frequency data reflects its more rapid time–distance variations.

(6) An image is constructed. The colour of each bin is determined by the polarization angle and its brightness is scaled by the rectilinearity (Figs 3and 4).

Figure 3

Time-distance variations of the Earth's average low-frequency (∼0.05 Hz) body wave polarization where the colour of each bin is determined by the polarization angle and its brightness is scaled by the rectilinearity. See text for details.

Figure 4

Time-distance variations of the Earth's average high-frequency (∼1 Hz) body wave polarization where the colour of each bin is determined by the polarization angle and its brightness is scaled by the rectilinearity. See text for details.

5 Results

The average polarization of the low- and high-frequency global seismic wavefield is shown in Figs 3and 4. The colour at a given range and time corresponds to the polarization angle θ and the brightness is proportional to the linearity ρ, Note that the angular colour scales differ between the high- and low-frequency stacks (lower right of Figs 3and 4). Although this complicates direct comparison between frequency bands, the different colour scales highlight the individual features appearing in each stack.

Only the polarization in the radial–vertical plane (P and SV energy) is considered. The polarization in the radial–horizontal plane (SH and SV energy) depends on the focal mechanism. Thus, the polarization is inconsistent between earthquakes and the stacking method will not work. The method is, however, independent of an arrival’s polarity.

In Figs 3and 4, regions with consistent body wave arrivals appear bright, a result of their linear particle motion. The average ground motion generated by seismic noise is approximately random and has a near-circular polarization that appears dark in the stacks.

Figs 5and 6show the IASP91 (Kennett & Engdahl 1991) ray-traced traveltime curves for arrivals visible in the two pass bands. Some of these phases find their way from hypocentre to seismograph through a complex combination of reflections, refractions and conversions.

Figure 5

The traveltimes for phases imaged in the low-frequency stack (Fig. 3) calculated using the IASP91 velocity model.

Figure 6

The traveltimes for phases imaged in the high-frequency stack (Fig. 4) calculated using the IASP91 velocity model.

Energy sampling deeper within the Earth arrives at angles closer to the vertical. This is evident for several phases imaged in Figs 3and 4. The polarization angle of the direct P wave decreases as range increases, and at identical ranges θ is smaller for P than for the shallower turning PP. Also, energy that turns in the outer core (PKPab) has a noticeably different θ than energy that turns in the inner core (PKPdf).

Seismic energy that travels a distance greater than 180° is recorded at a source–receiver range of less than 180°, For example, a station 90° away from an earthquake will record a PKKP arrival that has travelled 270° (see Fig. 6). Such ‘wraparound phases’ travelling between 180° and 360° arrive from a direction opposite to that defined by the receiver–source azimuth. For these distances, negative θcorresponds to compressional energy and positive θ corresponds to shear energy (opposite to the relation for distances less than 180°). The low-frequency PPP demonstrates this effect. PPP energy travelling less than 180° has a positive θ (seen as a yellow and orange streak in Fig. 3), whereas θ is negative for PPP arrivals travelling more than 180° (seen as a green streak in Fig. 3).

5.1 Low-frequency stacks

The low-frequency stack (Fig. 3) images the Earth’s polarization response at a dominant period of 20 s. In order to distinguish P- and S-polarized arrivals, the colour scale has sharp breaks at θ=0° and 90°, The compressional energy arriving between P and S is clearly imaged in magenta, red and yellow. This wedge consists of P, the surface multiples PP, PPP and PPPP, the core reflections PcP and PcS, and the core phases PKP and PKS. Also, the wraparound PKP, PP, PPP, SKSP and PKPPcP phases are seen near a range of 180° as green streaks with negative slopes. Another section of predominantly P-polarized energy lies between SKKS and SS. It consists of the surface S-to-P converted phases (e.g. SP, SPP).

The light blue streaks image shear energy arrivals, including S, SS, SSS, SSSS, ScS, SKS and SKKKS. The near-horizontal P-polarized energy immediately following the S and multiple S arrivals, known as shear-coupled PL waves, originates from coupling between the teleseismic SV energy and the fundamental leaking mode for the crustal waveguide (PL) (Oliver & Major 1960). The amplitude of the shear-coupled PL waves following direct S(SPL) diminishes at ranges beyond ~ 52°, where the P leg of the surface-converted SP arrival no longer undergoes a total internal reflection at the Moho. These arrivals have been used to determine uppermost mantle velocity structure (Zandt & Randall 1985; Zhang & Langston 1996).

5.2 High-frequency stacks

Fig. 4captures the polarization of the Earth’s seismic wavefield in the frequency band 0.5–2.0 Hz. Fewer phases are seen in these stacks, primarily a consequence of greater attenuation at high frequencies. The S and ScS arrivals are barely visible and the shear wave surface multiples are completely absent. On individual seismograms, the direct S phase is often linear and above the noise level. However, its polarization angle is often inconsistent from record to record, resulting in a near-circular average polarization. This is probably an effect of scattering or anisotropy.

The absence of multiple S waves and surface waves in the high-frequency stacks reduces the signal-generated noise later in the record section, uncovering the core phases SKKP, PKKP, PKKKP and PKPPKP. These arrivals are wraparound compressional arrivals and have a negative θ. Also, less contamination from P and PP coda produces a clearer PcP.

5.3 Predicted polarization

The observed and predicted P-polarization angles are compared in Fig. 7. The surface displacements of the incident, reflected and converted phases combine to give the observed ground motion. However, the particle motion of the reflected and converted phases is close to that of the incident wave. Thus, the polarization angle can be approximated from the definition of the ray parameter ( p=sini/V, where i is the angle of incidence and V is the velocity).

Figure 7

Difference between observed and predicted P-wave polarization at low and high frequencies.

In Fig. 7, the ray parameter is obtained from the IASP91 velocity model, and the best-fitting velocities of 7.5 and 5.2 km s^{− 1} for the low and high frequencies respectively were found by trial and error. The observed polarization angles were obtained using 1°× 15 s bins centred 10 s after the IASP91 theoretical P traveltime. These estimates were then smoothed using a 3° running average.

The observed low-frequency P-wave polarization angles are greater than the corresponding high-frequency measurements, a result of the deeper higher-velocity material sampled by the long-period energy. Also apparent is the well-fitted decrease in polarization angle as a function of epicentral range to 100° at both high and low frequencies. At ranges greater than 100° the scatter in the high-frequency observations increases significantly due to the rapid amplitude decay of P_diff.

This approach assumes an equivalence between polarization angle and angle of incidence. Thus, it does not account for the free-surface effect that alters the observed polarization angle from the true incidence angle. This correction can be up to a few degrees for the incidence angles considered here and depends on the Poisson Ratio (Bokelmann 1995). Corrections of this order are not needed here to show the general agree-ment between the observed and predicted polarization angles as a function of range and the difference between high- and low-frequency measurements.

The predicted SV polarization was calculated by the same technique used for P and is shown along with observations in Fig. 8. The observed and predicted polarization angles agree between 60° and 95°, At ranges less than 60° long-period SV waveforms are strongly distorted by SPL energy, causing a large deviation between the calculated and observed SV polarization. At ranges greater than 95°, contamination from SKS alters the observed polarization from expected values. This analysis clearly shows the distance range in which studies utilizing SV waveforms should be conducted.

Figure 8

Difference between the observed and predicted S-polarization angle. An S velocity of 4.1 km s⁻¹ was used to calculate the predicted angle, which is the average velocity in the top 85 km of the IASP91 earth model. Note the discrepancy between observed and predicted polarization at ranges of less than ∼ 60°. At these distances S-to-P converted energy is trapped in the crustal waveguide and the diference between the SV and SP arrival times is less than the period of the SV wave. At distances of less than 30° the observed polarization angle increases. However, the figure is truncated at 30° because the overlapping upper mantle triplicated arrivals complicate the calculation and measurement of the polarization angle.

6 Discussion

The record sections presented clearly show the average polarization of the Earth’s main body wave arrivals, demonstrating the effectiveness of multiple-event polarization stacking. The images provide a road map for detailed waveform studies by highlighting regions of constant arrivals and warning of contaminating phases.

This technique can be used to estimate the ray parameter of incoherent seismic wave trains, requiring only consistent polarization angles between seismograms rather than coherent waveforms. Thus, it may aid in the location and characterization of the Earth’s small-scale heterogeneities that generate high-frequency teleseismic wave trains (Haddon 1972; King 1975; Earle & Shearer 1997a).

When applying polarization stacking to small-amplitude arrivals, one must account for non-linear biases introduced by non-circularly polarized noise (Bataille & Chiu 1991). The back-ground noise in the high-frequency stack is vertically polarized, which will shift polarization angle estimates towards zero. This effect is negligible for large arrivals but must be corrected for in studies focusing on phases with low signal-to-noise ratios.

Polarization stacking offers an advantage over amplitude stacking by providing estimates of arrival angle and linearity; however, amplitude stacking methods remain superior at detecting extremely small phases. For example, upper-mantle discontinuity phases detected by aligning on reference phases prior to amplitude stacking (Shearer 1990) remain below the noise level in the polarization stacks. Also, the polarization stacks contain no information about relative amplitude, and thus cannot be used to obtain estimates of the impedance contrast across seismic discontinuities (e.g. Richards 1972).

Global body wave arrivals appear sharper in automatic gain control stacks (Astiz et al. 1996), yet these AGC stacks provide no information about the arrival angle and do not image the shear-coupled PL waves that are visible in Fig. 3. Polarization and amplitude stacking methods are clearly complementary.

Acknowledgements

Many colleagues have provided suggestions for the improvement of this study and I would particularly like to thank John Vidale, Vera Schulte-Pelkum, Steve Persh, Keith Richards-Dinger, Peter Shearer, Oz Pathare and two anonymous reviewers. I also thank the personnel of the Global Seismic Network for collecting the data, and the Incorporated Research Institutions for Seismology for archiving and distributing it. This work was conducted with funding from an NSF Earth Sciences Postdoctoral Research Fellowship.