Upper mantle attenuation and velocity structure from measurements of differential S phases

Summary

Although much progress has been made in determining the 3 D distribution of seismic wave velocities in the Earth, substantially less is known about the 3 D distribution of intrinsic attenuation. In this study variations in attenuation and shear velocity of the Earth's mantle are constrained using measurements of differential traveltime and attenuation.

The data are broad band displacement SH seismograms filtered to have energy in the period range 8–20s. The seismograms are obtained from over 600 globally distributed earthquakes of magnitude, M_w, 5.5 or greater.

Differential traveltimes and differential t* values from multiple S phases are estimated by a waveform fitting method, resulting in approximately 4300 measurements for SS–S in the distance range 50°–105° and 1000 measurements for SSS–SS in the distance range 90°–179°. Each measurement consists of a differential traveltime and a corresponding differential t*.

The differential traveltimes and t* values are inverted to obtain models of the lateral variation of shear velocity and lateral variation of q_μ, where q_μ=1/Q_μ. The models explain the data well but have limited depth resolution. The velocity models show good correlation with previous studies; in particular, low velocities are observed underlying mid oceanic ridges and convergent margins and high velocities are found for continental regions. The q_μ models show shield regions to be less attenuating than PREM, with ridges appearing as highly attenuating features. The models have limited depth resolution and to address this problem we also present a shear velocity model obtained from the combination of body wave and surface wave data sets.

1 Introduction

An understanding of the magnitude and distribution of the elastic and anelastic properties and density in the Earth's interior is important for a wide variety of geophysical, geochemical and astronomical studies. Since the late 1970s, seismic tomography has been used to construct models of the 3 D structure of the interior of the Earth (e.g. Dziewonski 1977; Woodhouse & Dziewonski 1984, 1986, 1989; Dziewonski 1984; Grand 1994; Masters et al. 1996; van der Hilst et al. 1997). Over the last two decades there has been a dramatic increase in both the quantity and the quality of data available. This, along with advances in computation, has allowed increasingly more detailed models of Earth structure to be constructed. The velocity structure of the Earth is now quite well known at long wavelengths, although some areas such as the transition zone and some sections of the lower mantle lack adequate resolution. However, substantially less is known about the 3 D distribution of intrinsic attenuation. The only previous global, 3 D study of body wave attenuation is that of Bhattacharyya et al. (1996). Studies of attenuation have also been performed using surface waves; these include Durek et al. (1993), Romanowicz (1990, 1994, 1995); and Selby (1998). Studies of the attenuation structure are of interest as a knowledge of attenuation structure can provide additional constraints on Earth properties such as the temperature and state of the material through which a seismic wave has travelled. However, this will not be possible until models achieve a much more detailed level of spatial resolution. Furthermore, improved models of seismic wave velocity should be possible if the effects of attenuation are taken into account.

Attenuation studies on regional length scales (e.g. across a continent) have been performed in a number of geographical regions (e.g. Barazangi et al. 1975; Sipkin & Jordan 1980; Lay & Helmberger 1981; Schlue 1981; Lay & Wallace 1983, 1988; Taylor & Bonner 1986; Chan & Der 1988; Sheehan & Solomon 1992; Sipkin & Revenaugh 1994). These studies provide higher resolution information in particular geographic regions. By incorporating results from regional models into global models it should ultimately be possible to obtain more refined and better constrained global models of attenuation.

A large number of studies have used multiple ScS body wave phases to obtain the whole mantle averaged shear attenuation beneath a particular region (Sipkin & Jordan 1980; Lay & Wallace 1983, 1988; Chan & Der 1988; Sipkin & Revenaugh 1994).

Frequency domain methods and the technique of phase equalization and stacking were used by Sipkin & Jordan (1980) and Lay & Wallace (1983, 1988). Sipkin & Jordan (1980) obtained values of Q_ScS for a variety of regions including the Pacific, South America and China. Lay & Wallace (1983, 1988) determined both ScS traveltime and attenuation values beneath Mexico and Central America and the western United States respectively. A time domain method of obtaining Q_ScS was developed by Chan & Der (1988). They used a waveform and amplitude matching scheme to examine the regional variation of Q_ScS using data in the frequency band 0.02–0.1Hz. The areas covered by their study included the Pacific, South America, Eurasia and North America. Sipkin & Revenaugh (1994) used a combination of time (waveform inversion) and frequency (phase equalization and stacking) domain techniques to estimate both differential traveltimes and the multiple ScS attenuation operator under China.

All of these studies of Q_ScS found large variations of attenuation across regions, implying large lateral variations in mantle attenuation. A good correlation with present day levels of tectonic activity was found (active tectonic regions producing high values of attenuation, shield regions producing low values of attenuation).

In this study, a new technique for making differential traveltime and attenuation measurements from shear waves is presented. Using such measurements we construct models of the variations in shear wave velocity and intrinsic shear attenuation q in the upper mantle.

2 Data

The data are digital, broad band, transverse component, long period displacement SH seismograms recorded over the period 1989–1996. Transverse component data are used for two reasons: they are not subject to the effects of S–P phase conversions and they eliminate the effects of possible interference of SKS near 80° (Woodward & Masters 1991). Long period data are used because the wavelengths involved allow large scale structure to be studied and also because the effect of fine scale structure is averaged out (Woodward & Masters 1991).

The seismograms are obtained from over 600 globally distributed earthquakes of magnitude, M_w, 5.5 or greater. The seismograms lie in the epicentral distance ranges 50°–105° and 90°–179° for the SS–S and SSS–SS data sets respectively. For the SS–S data set this distance range is selected so as to avoid the core shadow region for shear waves, and also, for epicentral distances greater than 50°, both the direct and the reflected phases sample the source and receiver regions in a similar manner. For the SSS–SS data set, the distance range is chosen in order to avoid rays bottoming in the transition zone.

We have analysed over 2900 events, selecting approximately 600 for further processing. The criterion for event selection is that both the S and SS or SS and SSS phases can be clearly distinguished from noise. For each selected event, the S, SS and SSS phases are windowed out interactively. The windowing is performed subject to the following criteria:

• (i) for a given station, the real and synthetic seismograms are similar;

• (ii) both the S and SS or SS and SSS phases must be clearly visible, that is, they must have a high signal to noise ratio.

Fig.1 shows several transverse component seismograms from an event of 1996 March 4. The ‘a’ and ‘b’ intervals denote the windows selected for the S and SS phases respectively. Also included in Fig.1 are examples of seismograms (stations KBS and KIV) for which the SS phase could not be identified.

3 Methodology

In the frequency domain the S and SS phases of the seismogram can be represented as a product,

where Ω(ω) represents the source time function, I(ω) is the instrument response function, a_S and a_SS describe the excitation and geometric spreading of the S and SS phases respectively and e^−iωtS and e^−iωtSS account for the traveltime delay. T(ω,t*) represents the attenuation and dispersion operator and is defined as (e.g. Chapman et al. 1988)

Re arranging eq.(1) and substituting into eq.(2), the following expression is obtained:

Now setting t*_SS−t*_S=Δt*, a_SS/a_S=a and t_SS−t_S=Δt and introducing a Hilbert transform due to a caustic encountered by SS(Aki & Richards 1980a,b), the following expression relating the S and SS phases is obtained:

Δt is the differential traveltime between the SS and S phases and Δt* is the differential t*. The factor i accounts for the fact that the SS undergoes a Hilbert transform relative to the S phase (Choy & Richards 1975).

For the SSS–SS differential measurements, the corresponding expression is

Δt₁ is the SSS–SS differential traveltime, Δt*₁=t*_SSS−t*_SS is the SSS–SS differential t* and a₁ is an amplitude factor relating the SSS and SS phases.

These expressions are used in conjunction with a non linear parameter estimation procedure (Bevington 1969) to obtain values for the differential traveltime and differential t* for SS–S and SSS–SS.

The parameters Δt* and a trade off against each other and as a result a two step procedure is employed to fit the waveforms. The algorithm is initially applied to synthetic data to obtain a value for a, the amplitude ratio between the synthetic S and SS phases. It is assumed that the real data have the same amplitude ratio and the algorithm is subsequently applied to the real data holding the value of a fixed. The synthetics are calculated by the summation of 20000 toroidal normal modes to a period of 6s using the method of Gilbert (1971). The modes are calculated using the code of Woodhouse (1988). The results from the second run of the fitting algorithm give the values for the SS–S (or SSS–SS) differential traveltime and differential t* relative to the PREM values. The fitting algorithm provides error estimates for Δt and Δt* and a. We also calculate a ‘fit parameter’ that is used as an additional indication of the quality of each estimate. The fit parameter is defined as the sum of the squared difference between the two waveforms divided by the sum of the squares of the known SS (or SSS) waveform. If a fit parameter of zero is obtained then the two waveforms are identical, and the closer to zero the fit parameter is, the higher the quality of fit between the two waveforms.

Fig.2 gives an example of the waveforms obtained from fitting the synthetic data. The left hand plot shows the traces prior to execution of the fitting algorithm and the right hand plot shows the traces after applying the algorithm. The parameters α₁, α₂ and α₃ give the starting values before applying the fitting algorithm. The parameters β₁, β₂ and β₃ give the final values of Δt, Δt* and a after running the fitting algorithm.

Fig.3 shows an example of the waveforms obtained by fitting the real data. As previously, the left and right hand plots show the waveforms before and after running the fitting algorithm respectively. The fit parameter for fitting the real data is 0.02, which indicates a very high degree of agreement between the real and calculated SS phases.

4 Results from the waveform fitting method

4.1 Traveltime residuals

The waveform fitting method yields a total of 4207 SS–S and 992 SSS–SS measurements of differential traveltime and t*. The residuals are given relative to the PREM value of differential traveltime or t*.

The distributions of SS–S and SSS–SS traveltime residuals are given in Fig.4. Both the SS–S and SSS–SS traveltime residuals have distributions that are close to Gaussian. The SS–S traveltime residuals have a mean of 0.12s and a standard deviation of σ=4.54s. The SSS–SS residuals have a mean of −0.81s and a standard deviation of σ=5.26s.

We can also examine the geographical distribution of the SS–S residuals by plotting them at the SS phase bounce point. For epicentral distances in excess of 50°, both the direct and reflected phases sample the source and receiver regions in a similar manner. Also, both the direct and reflected phases bottom at depths of 700km or greater. The mantle is thought to be more heterogeneous at the surface than at such depths and therefore the source of any observed traveltime or attenuation residual can be attributed to structure in the mantle or crust in the vicinity of the reflected (SS) phase's bounce point (Woodward & Masters 1991). The SSS–SS residuals cannot be represented in such a way.

Fig.5 shows the SS–S traveltime residuals plotted in 5° spherical caps at the SS bounce point. The spherical caps are plotted using the method of Woodward & Masters (1991). Any two caps having centres within 10° arc distance potentially have data in common. Positive residuals are plotted as crosses and negative residuals are plotted as triangles. Each spherical cap consists of the average of all the residuals having bounce points within a 5° radius of the centre of the cap. Plotting the data in spherical caps has the effect of smoothing the small scale variability and uneven geographical distribution of measurements. Positive residuals are observed along the mid oceanic ridges of the Pacific and convergent margins. Negative residuals are observed across the western Pacific and across continental areas.

4.2 t* residuals

The waveform fitting method provides 4207 SS–S and 992 SSS–SS differential t* (i.e. Δt*) residuals. The residuals are calculated relative to the PREM differential t* values. A positive residual is obtained when the differential attenuation obtained from fitting waveforms is greater than the PREM value and vice versa.

The distributions of both the SS–S and SSS–SSt* residuals (shown in Fig.6) are close to Gaussian. The SS–St* residuals have a mean of −0.18s and a standard deviation of σ=2.69s. The SSS–SSt* residuals have a mean of 0.36s and a standard deviation of σ=3.01s.

Fig.7 shows the SS–St* residuals plotted in 5° spherical caps at the SS bounce point. Positive residuals are plotted as crosses and negative residuals are plotted as triangles. Each cap contains at least three points, which results in 50 per cent coverage.

4.3 Data coverage

Fig. 8 illustrates the bounce point coverage for the combined SS–S and SSS–SS data sets. Each SS–S datum contributes one bounce point; each SSS–SS datum contributes three bounce points, two from the SSS ray and one from the SS ray. The bounce point coordinates are calculated by ray tracing. In Fig.8 the darker the shading, the higher the density of bounce points. White areas indicate regions where no bounce points occur. The bounce point density is highest across the Pacific region, with in excess of 350 bounce points occurring in some caps. The regions with the sparsest bounce point density occur across South America and Africa, as is to be expected from the global distribution of sources and stations.

5 Modelling

5.1 Model parametrization

We parametrize our models using spherical harmonic basis functions to describe lateral variations and spline functions for the depth dependence. A model of seismic anomalies (either velocity or attenuation) is then given by

where C_klm are the model parameters to be solved for, Y_l^m(θ,φ) are the spherical harmonics and f_k(r) are the radial basis functions (depth splines). The depth splines are defined from the Moho to the core–mantle boundary with the spline spacing increasing with depth. The unequal spacing allows variations at shallow depths to be more finely resolved. In this study we expect the anomalous signal to be attributed to sources in the shallow mantle and therefore choose a depth parametrization that results in sensitivity to approximately 400km depth.

The velocity models are constructed using harmonics of degree and order 12 and depth splines over six knots spanning the depth interval from the Moho to 400km, which results in a total of (12+1)²×6=1014 unknowns. Similarly, the attenuation models are constructed using harmonics of degree and order 8 using the same depth parametrization, resulting in a total number of 484 unknowns to be solved for.

5.2 Modelling of seismic parameters using weighted damped least squares

3 D models of the variations in attenuation and shear velocity are obtained by inverting the traveltime and attenuation residuals. The method of inversion is that of damped weighted least squares, where we seek to find the model m by solving the equation

where d is the data vector, m is the model vector and A is the matrix of derivatives relating d and m. Explicitly, the data vector is a column vector with the differential traveltime or attenuation residuals as its elements, i.e.

where each d_i is a single differential traveltime or attenuation residual. The solution to eq.(8) is given by

where D represents the inverse data covariance matrix (see Section5.5), U is an n×n matrix whose columns are the eigenvectors of A′^TA′ and Λ is defined as

where the λ_i represent the eigenvalues of A′^TA′ and η is a ‘damping’ parameter to achieve a compromise between model norm and variance ratio (see Section5.4) . A′ is defined as A′=DA.

The resulting model, m, will consist of a total of (l_m+1)²N_s parameters, where l_m is the maximum spherical harmonic degree and N_s is the number of depth splines used in the parametrization of the model. The model parameters are ordered spline by spline in order of increasing depth, e.g.

where each m_i is a column vector containing the (l_m+1)² spherical harmonic coefficients associated with depth spline i.

5.3 Data selection

Given a noisy data set, the process of data selection is important. The event selection process allowed 597 out of approximately 2900 events to be selected for further processing. However, further selection is necessary before inverting the data. The least squares inversion procedure requires a normal distribution of measurements so we must check that the data possess such a distribution. Least squares is particularly sensitive to outliers and therefore these should also be removed. Although it is important to remove the outliers, it is also important not to restrict the data set so severely as to significantly alter the shape of the normal distribution.

For data obtained using the waveform fitting method, only data points with fit parameters of 0.2 or less are selected for use in the inversion process as these are considered to be the most reliable. The threshold of 0.2 was selected subsequent to extensive visual examination of the waveform fits. This reduces the SS–S data set from 4207 to 2537 points and the SSS–SS data set from 992 to 387 points. The distributions of these reduced data sets are observed to be Gaussian with no removal of outliers required.

5.4 Damping

Damping is achieved by selecting a value of η that defines an eigenvalue cut off. This eigenvalue cut off can be chosen either from the value of η directly (λ_i=ηλ_max) or by examining the relationship between λ_i and the resulting variance ratio. λ_i denotes the ith eigenvalue and λ_max denotes the maximum eigenvalue. The effects of damping are best illustrated by considering the trade off between the variance ratio and model size. Fig.9 shows the relationship between the model size and the number of effective eigenvalues. The model size, |m|, is defined as

where m_i is the ith element of the model vector and n is the total number of model parameters. The number of effective eigenvalues, λ_e, is calculated from

where λ_k is the kth eigenvalue, η is the cut off parameter used to define where to begin tapering the eigenvalues, and λ_max is the maximum eigenvalue and is equal to λ₁.

Fig.10 shows the relationship between the variance ratio and the number of effective eigenvalues. The variance ratio, v, provides a measure of how well a given model fits the input data and is defined as follows:

These two figures illustrate the trade off between the variance ratio and model size. As the number of effective eigenvalues is increased the model size increases gradually until around η₁₂, where it increases rapidly. Conversely, the variance ratio decreases as the number of effective eigenvalues increases. Ideally the ‘best’ model would have as small a variance ratio as possible. However, from an inspection of Figs 9 and 10 it is clear that a compromise must be reached. In order to minimize the variance ratio whilst keeping the model size constrained we use a value of η=0.005 in the inversion procedure.

5.5 Data weighting

In a least squares inversion the data error for all measurements should be equal to ensure that the resulting model is not distorted by poor measurements. The data errors may be set to the standard deviations obtained from the measurement algorithm. Following the inversion the data error can be estimated a posteriori.

If we use these error estimates to define the data weighting for the inversion procedure, the value of χ² obtained for each data point is much greater than 1.0. We conclude that the errors provided by the waveform fitting method are greatly underestimated. Therefore, we weight all data points equally (equivalent to setting all values of σ=1.0), which reduces the inversion problem to that of damped least squares. This means that the matrix D in eq.(10) is proportional to the identity matrix.

5.6 Corrections to the traveltime residuals

The traveltime residuals must be corrected for the effects of the Earth's ellipticity and varying crustal structure. These corrections are applied prior to inversion. The ellipticity corrections are made using the methods described by Kennett & Gudmundsson (1996). They used the approach of Dziewonski & Gilbert (1976) to provide tabulated ellipticity coefficients for a wide range of seismic phases including S and SS. They also provided a theory for combining different phases, which allows the ellipticity correction for SSS to be computed.

The crustal corrections are made using model CRUST5.1 of Mooney et al. (1998). Since we are concerned with differential measurements, we need only account for the effects of the crust at the reflected phase's bounce points. For example, for the SS–S measurements the crustal correction is calculated only at the SS bounce point. In general, the crustal corrections cause a reduction in differential traveltime for paths with bounce points in continental regions and an increase in differential traveltime for paths with oceanic bounce points. The effect is most pronounced for the SS–S data set as only a single bounce point is involved. For the SSS–SS data set the effects are more complex as the crustal correction is calculated from the two SSS bounce points and a single SS bounce point.

5.7 The models

Fig.11 shows velocity models MVSS and MVSSS plotted at a representative depth of 190km. Model MVSS is the model obtained from inverting the SS–S traveltime residuals and model MVSSS is the model obtained from inverting the SSS–SS traveltime residuals. Each map shows the lateral variation in δv/v_s given as a percentage variation from PREM. Hotspots are plotted as green circles.

Fig.12 shows attenuation models MQSS and MQSSS plotted at a representative depth of 190km. MQSS is obtained from inverting the SS–St* residuals and model MQSSS is obtained from inversion of the SSS–SSt* residuals. Each map shows the lateral variation in Δq. The blue colours correspond to areas in which the attenuation is less than in PREM and the red colours correspond to areas in which the attenuation is greater than in PREM.

5.7.1 Combining the SS–S and SSS–SS models

We construct a joint SS–S and SSS–SS model by combining the derivative matrices for the SS–S with those for the SSS–SS data set. This combined derivative matrix is inverted, resulting in a velocity (or attenuation) model that contains information from both the SS–S and SSS–SS data sets. The number of SS–S residuals outnumbers the number of SSS–SS residuals by a factor of 6 and therefore we should expect that the joint model will be dominated by the SS–S data. A consequence of this is that the joint model will generally fit the SS–S data set far better than it fits the SSS–SS data set. If the number of SS–S and SSS–SS residuals were equal, we should expect the resulting joint model to fit the two data sets equally well. However, since this is not the case a compromise is sought. Essentially we weight the SS–S and SSS–SS data sets in such a way as to produce a joint model that fits both data sets equally well. Fig.13 illustrates the weighting procedure. Fig.13 shows the variance ratios obtained for the SS–S and SSS–SS attenuation data sets (computed against the joint model) for different weightings. The weightings w₁, w₂ are used as multipliers for the derivative matrices for the SS–S data and SSS–SS data respectively. For Fig.13, w₁ and w₂ take values ranging between 1000.0 and 1.0. Point C (Fig.13) corresponds to a high relative weight for SS–S (where w₁/w₂=1000.0) so that the resulting model is essentially the same as a model constructed from SS–S data alone. Point D corresponds to w₁/w₂=0.001, so that the resulting model is essentially the same as a model obtained using SSS–SS data alone. Point A corresponds to equal weighting (i.e. w₁/w₂=1.0) of the two data sets. Point B is chosen by calculating the ratios v₁/v_1min and v₂/v_2min, where v₁ and v₂ are the variance ratios for the SS–S and SSS–SS data sets for a given choice of weighting and v_1min and v_2min are the minimum variance ratios for the SS–S and SSS–SS data sets (e.g. for Fig.13v_1min=0.81 and v_2min=0.59). We take point B to represent the best compromise between the two data sets.

Figs 14 and 15 show the resulting velocity and attenuation models denoted MVCOMB and MQCOMB, respectively, obtained from combining the SS–S and SSS–SS data sets. The variance ratios obtained for models MVCOMB and MQCOMB are respectively 0.33 and 0.77.

6 Resolution tests

We test the resolution of our models by performing a simple resolution test. Essentially the procedure is as follows:

• (i) generate the predicted data from the test model using the derivative matrices calculated for the real data set (using the combined SS–S and SSS–SS data set);

• (ii) invert the predicted data using the same damping and weighting as for the real data set to obtain the output model.

The test model is a checkerboard pattern generated for l=12 and m=6 at a depth of 256km. The input or ‘test’ model is given by Fig.16.

Fig.17 shows the output model from the resolution test. The same checkerboard pattern is now observed across all depths with approximately equivalent amplitude. This leads to the conclusion that while the lateral resolution of our models is good (the checkerboard pattern is reproduced successfully), the depth resolution is poor (the pattern is spread across all depths). This means that we cannot attribute the lateral variations of δv/v_s or Δq observed in the models to a particular depth. Results for different target depths are essentially the same as those shown in Fig.17.

7 Velocity model from combining body wave and surface wave data sets

The resolution testing illustrates that although our models have good lateral resolution, the depth resolution is somewhat limited. In order to obtain improved depth resolution we combine the body wave data set of this study with the surface wave data set of van Heijst & Woodhouse (1999). Combining the two data sets also allows greater confidence to be assigned to the resulting model, providing the data sets are compatible. An earlier version of this model that included only the SS–S body wave data is presented in van Heijst et al. (1998).

The body wave data set is identical to that used for constructing model MVCOMB and consists of a total of 2921 measurements of differential traveltime, 2535 from SS–S measurements and 386 from SSS–SS measurements. No weighting of the individual SS–S and SSS–SS data sets is performed. Both ellipticity and crustal corrections from CRUST 5.1 are applied to the traveltime residuals.

The surface wave data set consists of close to 1000000 Rayleigh wave fundamental mode and overtone dispersion measurements that are made for minor arcs (van Heijst & Woodhouse 1999) and major arcs. A full description of the measurement technique is given by van Heijst & Woodhouse (1997). The surface wave data set is corrected for the effects of the crust using CRUST 5.1. Corrections for the Earth's ellipticity are also applied to the surface wave data before inversion. The relative weighting of overtone and traveltime data was chosen to achieve satisfactory variance reductions for each data set.

The body wave data are mainly sensitive to structure in the upper mantle but they also have some sensitivity to structure in the lower mantle due to the turning of the S rays in this region. Because of the sparse ray coverage in the lower mantle we restrict discussions to the top 1200km of the mantle. The surface wave and overtone data have sensitivity to the upper 1200km with reduced sensitivity to deeper structure. Fig.18 shows the resulting velocity model, denoted MVBS, over a range of depths from 50–1094km.

8 Discussion

The velocity models MVSS and MVSSS (Fig.11) exhibit many similarities. The main difference between them is that the mid Atlantic ridge is observed as a low velocity feature in model MVSSS and is neither high or low in model MVSS. Also, model MVSSS shows the East African rift valley to be a low velocity feature, whereas model MVSS shows this to be a higher than average velocity feature.

Model MVCOMB is obtained by inversion of the combined SS–S and SSS–SS data sets. The velocities along the faster spreading ridges such as the East Pacific Rise are lower than the velocities along slower spreading ridges such as the mid Atlantic ridge. Areas of low velocity are also observed to the North of the African rift valley and in a region across Asia to the northeast of India. Areas of high velocity correlate with the continents; for example, Australia, South America, much of Eurasia, Africa and the northeast of America and Canada are all observed to have higher than average velocities. Much of the Atlantic Ocean is also observed to have higher than average velocity. The velocities across the Pacific ocean are observed to increase westwards with distance from a spreading centre.

The lack of depth resolution means that it is impossible to attribute the regions of high or low velocity to particular depths. Many other studies of lateral variations in upper mantle shear velocity have observed such correlation with surface tectonics (e.g. Woodhouse & Dziewonski 1984; Su et al. 1992; Grand 1994; Masters et al. 1996).

Model MQSSS has many features in common with model MQSS; however, there are several differences. The region of extremely low attenuation to the west of Africa is considerably more intense in model MQSSS and South America appears as a strong low attenuation feature in model MQSS whereas it appears as a weak high attenuation feature in model MQSSS. Model MQSSS has an area of very high attenuation to the east of southern Africa; this feature also appears in model MQSS but it is not as intense. Both models show a region of high attenuation to the northeast of India; however, this region lies further south in model MQSSS. Each of these discrepancies can be attributed to data coverage. The bounce point coverage across these areas as depicted in Fig.8 is particularly sparse and therefore caution should be used when interpreting the model in these regions.

Model MQCOMB is the preferred q_μ model obtained from inverting the SS–S and SSS–SS data sets simultaneously. It exhibits features from both model MQSS and model MQSSS. Areas of high attenuation, with the exception of South America, correlate with spreading centres and convergent margins. For example, high attenuation is observed along the East Pacific Rise, Pacific–Antarctic Rise, Macquarie Ridge and Kermadec–Tonga Trench. Continental areas generally correlate with regions of lower than average attenuation, with the exception of South America, which appears as a region of higher than average attenuation. Low attenuation is observed across Australia, Antarctica and much of Eurasia, Africa and North America. The amplitude of the variations in q lies between ±0.020.

Both the SS–S and SSS–SSt* residuals (with respect to PREM) have a non zero mean. This means that the models of q_μ have a significant degree 0 term. Since the PREM Q structure is based on very long period modal data, we believe that the degree 0 term may represent evidence of the frequency dependence of q_μ in the frequency range 15–150s, which should be the subject of future study. The amplitude of the q_μ variations and the fact that physically Q cannot be negative means that the probable source of the q_μ anomalies is in the low velocity zone (depths 80–220km). This is consistent with the observations of Durek et al. (1993), who suggested that the source region of anelastic heterogeneity lies largely in the shallow mantle between depths of 100 and 300km.

Unlike the models obtained from inversion of the body wave data sets, model MVBS has good depth resolution to depths of at least 1000km. This is due to the inclusion of the surface wave and overtone data set. At 50km depth, model MVBS is characterized by low velocities along mid ocean ridges and convergent margins. These low velocities are also observed at 129km but the intensity is lower. The continents are observed as high velocity features to depths of around 256km and possibly deeper, with the maximum intensity observed around 129km. The African rift region is observed as a low velocity feature to depths of around 329km. Over the depth range 409–703km the model is dominated by the high velocity anomalies associated with subduction, and the central Pacific remains an area of low velocity. Between depths of 951 and 1251km the largest anomalies are high velocities in northern Africa, eastern Asia and North America.

8 .1 Comparison of the velocity and attenuation models MVCOMB and MQCOMB

From a physical point of view, if both velocity variations and q variations are due to temperature variations, it would be expected that regions of high attenuation would correlate with regions of low velocity. Comparison of models MVCOMB and MQCOMB leads to the following conclusions:

• (i) areas of low attenuation tend to correlate with areas of high velocity, for example, across Australia, Eurasia, Eastern North America, the Middle East, Antarctica and the Eastern Atlantic;

• (ii) areas of high attenuation tend to correlate with areas of low velocity, for example, along the ridges surrounding the Pacific Ocean, along the Atlantic–Indian ridge and in a small region to the northwest of India.

Global correlation coefficients for spherical harmonic degrees 1–8 are greater than 0.25 at all depths, significant at the 97.7 per cent level. We also observe several areas where the two models are not in agreement, for example, South America, an area to the east of Africa and the central Pacific. Substantial differences between the velocity and attenuation maps are an indication that velocity and attenuation variations reflect both temperature and compositional variations.

8.2 Comparison of models MVCOMB, MQCOMB and MVBS with existing models

Although only SS–S and SSS–SS data are used in the construction of model MVCOMB, some striking similarities are observed between it and existing models. Upper mantle velocity models are relatively well determined and differ little from one another. Here we make comparisons with models M84C (Woodhouse & Dziewonski 1984), SAW12D (Li & Romanowicz 1996), S16B30 (Masters et al. 1996) and S20RTS (Ritsema et al. 1999). Comparison of model MVCOMB (Fig.14) with depths 50 and 200km of model M84C of Woodhouse & Dziewonski (1984) shows many similarities. The correlation of high and low velocity anomalies with tectonic features is observed by each model. The relative intensity of the high velocity anomaly across Africa is less in model MVCOMB than in model M84C at 200km. However, the amplitude of this high velocity feature is approximately the same in each model. Comparison of model MVCOMB with depths 150 and 250km of model SAW12D of Li & Romanowicz (1996) again shows many similar features. The most significant difference between the two models lies in the region around the western Pacific rim. This feature is observed as a high velocity anomaly in model SAW12D but as a low velocity anomaly in model MVCOMB. Comparison with depths 70 and 170km of model S16B30 of Masters et al. (1996) also shows many similarities.

Comparison of model MVBS (Fig.18) with existing models of shear velocity perturbation shows many similarities. In the upper mantle, model MVBS shows correlation with tectonic features: high velocities underlying continental regions and low velocities underlying spreading ridges. Similar features have been observed by many studies, in particular Woodhouse & Dziewonski (1984), Su et al. (1994), Li & Romanowicz (1996), Masters et al. (1996) and Ritsema et al. (1999). In particular, it is closely in agreement with S20RTS (with which it shares the overtone data).

We now make brief comparisons of the attenuation model MQCOMB with existing models. Comparison of model MQCOMB with model QR19 of Romanowicz (1995) and the shear attenuation model of Bhattacharyya et al. (1996) show limited agreement. Comparison with model QMU3b of Selby (1998) gives rise to many similarities. The compatibility of models QMU3b and MQCOMB has already been demonstrated by Reid et al. (1999).

We also make a comparison of model MQCOMB with several regional studies. The studies of Lay & Helmberger (1981), Taylor & Bonner (1986) and Chan & Der (1988) all find the North American region to have high attenuation across tectonically active areas (e.g. the Basin and Range Province) with low attenuation across shield areas. Model MQCOMB is in general agreement with these conclusions.

Sheehan & Solomon (1992) studied the North Atlantic region using measurements of differential SS–S attenuation and traveltime. They found SS–S attenuation measurements to be positively correlated with differential traveltime residuals. They also found that both the SS–S traveltime and the attenuation residuals decrease with increasing seafloor age. Model MQCOMB shows limited agreement with these observations; however, the data coverage across this region is very sparse (see Fig.8).

Sipkin & Revenaugh (1994) investigated attenuation and traveltime variations across China. They found regions of low attenuation across the Tarim platform and the fold systems of northeastern China, with high attenuation elsewhere on the continent. Model MQCOMB shows China to be a predominantly low attenuation feature but does not have sufficient resolution to distinguish between the regions discussed here.

Finally, we compare model MQCOMB with the work of Barazangi et al. (1975). They used pP phases to map variations in the attenuation of high frequency (0.5–2Hz) compressional waves in the wedge of mantle lying above nearly all the inclined seismic zones on Earth. They found high attenuation behind the Tonga, New Hebrides, Mariana and Japanese island arcs and also across the Andean Altiplano of western South America. Low attenuation features were observed across the Indonesian, Philippine, Izu Bonin and New Britain–Solomon regions. With the exception of the Japanese island arcs, model MQCOMB is in excellent agreement with the results of this study.

8.3 Geophysical implications

In summary, it is found that regions of low velocity tend to correlate with the spreading ridges and convergent margins. Regions of high velocity tend to correlate with continental regions. The velocities observed across the Pacific increase westwards, related to the aging and therefore cooling of the oceanic lithosphere with distance from a spreading ridge. The ridges and convergent margins are characterized by high attenuation features with low attenuation generally observed underlying continental areas.

9 Conclusions

In this study, a new method for measuring the differential properties of multiple S phases has been presented. The ‘waveform fitting method’ enables measurements of both SS–S and SSS–SS differential traveltime and attenuation to be made. Using the results from the waveform fitting method, models of the lateral variations of upper mantle shear velocity and intrinsic attenuation were calculated. To improve the depth resolution of our velocity model we also inverted a combined body wave and surface wave data set.

Model MVCOMB is one of the first models of upper mantle shear velocity perturbations to include differential SSS–SS measurements. Despite the small number of measurements, the SSS–SS measurements show great promise. They provide very good sampling of the upper mantle as there are three bounce points for each source and receiver pair: one arising from the SS ray and two arising from the SSS ray. Although the depth resolution of model MVCOMB is limited, the features observed are consistent with earlier observations: low velocities are observed underlying mid oceanic ridges and convergent margins, with high velocities underlying continental regions. Comparisons with existing models of shear velocity heterogeneity showed many similarities.

Model MQCOMB is the fourth truly 3 D model of upper mantle intrinsic attenuation after the models of Romanowicz (1995), Bhattacharyya et al. (1996) and Selby (1998). The body wave model of Bhattacharyya et al. (1996) was found to be somewhat incompatible with earlier surface wave results. The results presented here, on the contrary, display a significant degree of similarity to recent results from surface wave modelling (Selby & Woodhouse 1999; Selby 1998). Results obtained by combining surface wave data and the body wave measurements show that models can be found that reconcile the two kinds of data; these models will be reported elsewhere. The amplitude of the degree 0 term provides, we believe, evidence of the frequency dependence of q in the frequency range of this study (≈15s) and that of the mantle wave data of previous global studies (≈150s). This is a significant finding that should be further researched.

Model MVBS, obtained from combining body wave and surface wave data sets has good depth resolution due to the surface wave data with sensitivity in the lowermost mantle arising from the turning of the deepest body wave S rays. In this region there is excellent agreement with existing studies (Woodhouse & Dziewonski 1989).

Acknowledgments

We gratefully acknowledge the research studentship awarded to FJLR by the Natural Environment Research Council. We wish to thank Neil Selby for numerous useful discussions.

References

Author notes