Databases of surface wave dispersion

We derive phase velocity maps from two global databases of surface wave dispersion, and compare them to assess the relative quality of the observations. Although based upon similar sets of seismic records, the databases show some significant discrepancies, which we attempt to trace to the different automated measurement techniques employed by their authors. This exercise is relevant to our current understanding of the Earth’s upper mantle, whose elastic structure is most importantly constrained by surface wave observations like the ones in question.


Introduction
Seismic surface waves are dispersive; their speed of propagation is a function of their frequency, or, in other words, individual harmonic components of the surface wave seismogram (individual modes) propagate over the globe at different speeds.For any given earthquakeseismometer couple, one can isolate from the seismogram each harmonic component, and measure its average speed (phase velocity) between source and receiver.We call «dispersion curve» a plot of average phase velocity against frequency.From a large set of measured dispersion curves, with sources and stations providing as uniform a coverage of the Earth as possible, local phase velocity heterogeneities can be mapped as a function of longitude and latitude.At each location, perturbations in phase velocity with respect to a given reference model are weighted averages of seismic heterogeneities in the underlying mantle (Jordan, 1978).
Through the steps outlined above, measurements of surface wave dispersion provide the most important global seismological constraint on the nature of the Earth's upper mantle (e.g., Ekström, 2000).Over the last decades, two large intermediate-period (35-150 s) fundamental mode dispersion databases have been assembled by Trampert andWoodhouse (1995, 2001) and Ekström et al. (1997) from global networks of broadband seismic stations, and are now available to the public.Both are currently employed in global tomography studies (e.g., Boschi and Ekström, 2002;Spetzler et al., 2002;Antolik et al., 2003;Beghein, 2003;Boschi et al., 2004;Boschi, 2005); their accuracy becomes even more important when they are used to constrain observables that are inherently hard to resolve, like the azimuthal anisotropy of surface wave propagation (Laske and Masters, 1998;Becker et al., 2003;Trampert and Woodhouse, 2003).Boschi and Woodhouse (2005) point to the substantial disagreement between published global models of surface wave anisotropy, and the need to explain it partly motivates the present study.
The Utrecht and Harvard databases are in good, but not complete agreement; it has been observed by their authors, and we again illustrate below, that those of Trampert andWoodhouse (1995, 2001) require phase velocity maps with a whiter spectrum than those of Ekström et al. (1997), to be explained comparably well.This discrepancy might depend on the a priori selection of globally recorded seismic records in the two studies: Ekström et al. (1997) obtained phase velocities from digital seismograms from the Global Seismographic Network (GSN), using earthquake sources between 1989 and 1995, Trampert andWoodhouse (1995, 2001) used digital GDSN and GEOSCOPE seismograms recorded between 1980 and 1990.The discrepacy could also arise from the differences in the measurement algorithm designed and employed by the two groups.In both cases theoretical («synthetic») seismograms are calculated from theoretical dispersion curves, and compared to the measured one.An inverse problem is then solved to find the theoretical dispersion curve that minimizes the misfit between measured and synthetic seismograms.Such optimal dispersion curve is the dispersion measurement sought.As we shall illustrate below, the two groups implemented this procedure in different ways.
Thus far, 2-d phase velocity and 3-d shear velocity maps have been derived independently from the two databases by different authors, making use of different tomographic algorithms; Ritzwoller et al. (2001) and Shapiro and Ritzwoller (2002) merged all data in a single set that they inverted, but did not provide a comparative analysis of data quality.
We present two sets of phase velocity maps, derived independently from the two databases, with the same parameterization, inversion algorithm and regularization scheme.At the period where discrepancies are most severe (35 s Love waves), we refine the parameterization and explore the solution space, in an attempt to determine whether Trampert and Woodhouse's (2001) database can indeed resolve coherent isotropic heterogeneities of high spatial frequency.

Measuring surface wave dispersion
Following Ekström et al.'s (1997)  with X denoting propagation path length measured along the great circle.u(ω) can be calculated, e.g., by surface wave JWKB theory as described by Tromp andDahlen (1992, 1993).Given an observed seismogram, surface wave dispersion between the corresponding source and receiver is measured by finding the curves A(ω) and c(ω) (hence Φ(ω)) that minimize the misfit between the observed seismogram and u(ω).The best-fitting c(ω) is the seeked dispersion curve.Ekström et al. (1997) and Trampert andWoodhouse (1995, 2001) find synthetic seismograms in slightly different, but equivalent, ways.Both groups write the corrections to amplitude and phase, with respect to PREM (Dziewonski and Anderson, 1981), as linear combinations of a set of cubic splines fi(ω), with coefficients bi and di, respectively, (2.3) (2.4) and the dispersion measurement is reduced to the determination of bi, di (i=1, ..., N).
The outlined inverse problem is highly nonlinear.Ekström et al. (1997) accordingly chose to solve it through the non-linear SIMPLEX method (e.g., Press et al., 1994).They parameterized their dispersion curves in terms of only N = 6 splines; this number being relatively small, the curves tend to be smooth and no further regularization is needed.Conversely, Trampert andWoodhouse (1995, 2001) described dispersion curves as linear combinations of as many as N=36 splines.They accomodated the nonlinearity of the inverse problem by first measuring group velocities, which are «secondary observables easily calculated by the seismograms».They then derived phase velocity dispersion curves from them, through the nonlinear least-squares algorithm of Tarantola and Valette (1982), regularized by some roughness damping constraint.
The main reason for requiring dispersion curves to be smooth, is to avoid fictitious fluctuations due to «2π» ambiguity: notice from eq. (2.1) that the surface wave seismogram is a periodic function of Φ with period 2π ; when dispersion is measured, values of Φ different of exactly 2π will provide an equally good fit to the data; and perturbations δΦ larger than 2π cannot simply be discarded, because at short to intermediate periods (say, <100 s), heterogeneities in the Earth's structure can suffice to cause phase anomalies of more than one cycle.Both groups solved this ambiguity requiring the dispersion curve to be smooth (see, e.g., fig. 2 of Trampert and Woodhouse, 1995); Ekström et al. (1997) noted that at longer periods (>100 s) a phase measurement can be associated with a total cycle count with no difficulty; they therefore preferred to first measure dispersion in the longer period range, and extend the measured dispersion curve to shorter period in a subsequent step, with the requirement that it be smooth, as they illustrated in their fig.3.

Intermediate resolution maps at different periods
We derive Love and Rayleigh wave phase velocity maps at all periods (35 to 150 s) where measurements are available; a separate tomographic inversion is performed at each period to find an independent phase velocity map; we then compare maps obtained from Trampert and Woodhouse's (2001) database with those obtained from Ekström et al.'s (1997).Both databases have, at all periods, a relatively good coverage of the entire globe.We employ, in both cases, the same tomographic parameterization (approximately equal area pixels, of size 5°×5°degree at the equator), regularization (a combination of norm and roughness damping), and inversion algorithm (LSQR) (Boschi and Dziewonski, 1999;Boschi, 2001): differences between the maps must then arise from discrepancies in the data.The regularization scheme was chosen after a number of preliminary inversions (Carannante, 2004).
The phase velocity distributions that we found share the general character of already published maps.For any given regularization scheme, inversions of Harvard data lead to a systematically higher variance reduction than those of Utrecht ones, also confirming published results.
The correlation between maps improves with increasing period, and is systematically higher for Rayleigh, with respect to Love waves.Figure 1 allows a visual comparison between Love and Rayleigh wave velocity as imaged from the two databases, at short (35 s) and intermediate (∼75 s) periods.In fig.2, we show the harmonic spectra of the same images, derived as described, e.g., by Boschi and Dziewonski (1999).The agreement between the images is remarkably good at ∼75-80 s, both in the spatial and frequency domains (comparing the third and fourth panels of fig.2, one might notice that, as anticipated, the Rayleigh wave spectra are slightly more similar to each other than the Love wave ones); high correlation between images based on the Trampert and Woodhouse's (2001) and Ekström et al.'s (1997) databases is confirmed at all measured intermediate and long periods (Carannante, 2004).
As period decreases, a loss of correlation is evident; it is maximum at 35 s (measurements at shorter periods are not provided in the databases), both for Love and Rayleigh waves, but stronger in the case of Love waves.The 35 s Love wave maps disagree most evidently in the regions of thickest crust (Andes, Tibet), where Ekström et al.'s (1997) data provide a much lower phase velocity than those of Trampert and Woodhouse (2001).The top panels of fig. 2 show that the discrepancy has a global character, with the former images being characterized by a much redder spectrum than the latter.

Higher resolution, short-period phase velocity maps
The exercise described in the previous section confirms that a discrepancy exists between the databases in question.The discrepancy is practically negligible at intermediate to long periods, but becomes more severe with decreasing period.In a short-period regime, tomographic inversions of the two databases result in phase velocity maps that would have, at least in certain regions, quite different implications on the properties of the upper mantle.
To investigate the nature of the discrepancy, we focus our attention on Love wave observations at 35 s, where the most evident disagreement is found.We refine our grid, reducing the pixel size to 2°×2°at the equator; this parameterization, involving ∼10 000 free parameters, should be sufficient to resolve any short-wavelength structure «seen» by the data.This is an important issue: if either database contains information on the structure of very high harmonic degree, that a 5°×5°grid cannot represent, aliasing of this signal could in principle be a reason for discrepancies between images.
We must likewise rule out the possible inadequacy of regularization as a reason for discrepancy: we repeat our inversions with a number of different regularization schemes.A few examples are shown in fig.3; at each damping level, the phase velocity distributions derived from the two databases differ, with a pattern similar to that outlined by our lower resolution inversions.The spectral plots of fig. 4 are also qualitatively consistent with those in the top panel of fig.2, except for the obviously underdamped case (norm and roughness damping parameters equal to 0.005), where the large l harmonic coefficients clearly take non-physical values.
The most outstanding difference between maps derived from Ekström et al.'s (1997) and Trampert and Woodhouse's (2001) databases is found, again, in regions of anomalously thick crust.Regardless of regularization scheme and parameterization, Tibet and Andes are coherently slow according to the former database, locally fast according to the latter.We compare the results of our inversions with a theoretical 35 s Love wave phase velocity map (bottom of fig. 3), computed, following Boschi and Ekström (2002), from a spherically symmetric mantle model combined with the laterally heterogeneous crustal model Crust-2.0 (Bassin et al., 2000) derived from observations entirely independent from the ones discussed here.This simplified model favours uniformly low phase velocity at Tibet and Andes; for Trampert and Woodhouse's (2001) observations to be correct, shear velocity then has to be very high in the upper mantle underlying those regions (see, e.g., Boschi and Ekström, 2002, fig. 7), while a simpler mantle model would be sufficient to explain Ekström et al.'s (1997) data.
Given the systematic nature of the discrepancy, always stronger in regions of thick crust, we could infer that either shorter-period signal (more sensitive to crustal thickness) «leaks» into the 35 s Love wave measurements of Ekström et al. (1997), or longer-period signal (more sensitive to upper mantle structure) «leaks» into Trampert and Woodhouse's (2001).Figure 5, where we show correlation between maps at a set of neighbouring periods, is to some extent consistent with the latter hypothesis, while the former cannot be tested in the absence of global maps of phase ve-locity at shorter periods.In both cases, discrepancies would be ascribed to faults in the measurement algorithms.
The high spatial frequency character of maps derived, at this period, from Trampert and Wood- house's (2001) data suggests, alternatively, that their data might simply be noisier.This hypothesis is difficult to test.Carannante's (2004) histograms of both databases at all periods have not disclosed any clear systematic difference between the Utrecht and Harvard databases.et al.'s (1997) databases (solid and dashed lines, respectively).Subsets of the two databases at Rayleigh (top) and Love (bottom) wave periods approximately 35 s, 40 s, 50 s, 60 s, 75 s, 100 s and 150 s are inverted with the 5°×5°pixel parameterization of Section 3.1; the correlation of each of the resulting maps with the one, from the same database, at the neighbouring longer period is shown as a function of period.Correlation is computed according to Becker and Boschi (2002), and all harmonic degrees up to 40 are considered.For both Love and Rayleigh waves at short periods, maps from the Utrecht database are better correlated with each other than those from the Harvard one; the opposite is true at longer periods.
In principle, scattering and anisotropy in surface wave propagation are possible causes for aliasing and subsequent discrepancies between databases: an accurate dispersion measurement could contain, for example, a coherent «aniso-tropic signal» that might have been wrongly removed in a less accurate one.Scattering and anisotropic effects are also perioddependent.Again, we are not yet capable of testing this hypothesis; attempts at mapping the azimuthal anisotropy of surface wave phase velocity have not led to consistent results (Laske and Masters, 1998;Becker et al., 2003;Trampert and Woodhouse, 2003; see also Boschi and Woodhouse, 2005 for a comparison); Boschi (2005) notes that accounting for scattering effects in inversions of the Harvard database appears to lead to an effective improvement in resolution, but, at the present stage, his results are only preliminary.

Further exploration of the solution space
We find the «L-curves» (e.g., Hansen, 1998) associated with the Love wave 35 s, 2°×2°pixels inverse problem (Section 3.2): in fig.6 (top) data misfit (1-variance reduction) is plotted as a function of total roughness (eq.2.21 of Boschi, 2001), for a number of phase velocity maps derived from the Utrecht (solid lines) and Harvard (dashed lines) databases, with roughness minimization as the only damping constraint.Roughness minimization, with respect to other damping schemes, limits the effect of the starting model on the final solution (e.g., Inoue et al., 1990, Section 3.3.1 and fig.2).Values of the roughness damping parameter range from 0.005 (minimum misfit) to 0.5 (maximum misfit).We also take the derivative of misfit with respect to total image roughness (fig.6, bottom), to make the two L-curves comparable.
If a small increase in image complexity is sufficient to reduce the misfit significantly, one must infer that meaningful signal is being mapped, and that the image resolution is also being improved.Conversely, if misfit remains approximately constant even for large increases in total roughness, the regularization must be inadequate, so that the effect of noise in the data, and of subsequent nu-merical instabilities, become prominent in the solution.Ideally, the vertex of the «L»'s in curves like those of fig.6 should identify the optimal values of total roughness and damping parameter, leading to the lowest possible misfit while minimizing the fictitious mapping of noise.
The L-shape of the curves in fig.6 indicates that the data are of sufficiently high quality to resolve meaningful Earth structure; it is hard, however, to identify an optimal value for the damping parameter, as the high-resolution grid employed here is probably fine enough to allow a certain amount of noise to be mapped into fictitious phase velocity heterogeneity.The curves take a more pronounced L-shape in the case of  21)).Misfit is defined as 1 minus the variance reduction.
Harvard data, easier to explain in terms of isotropic phase velocity variations.
After numerous independent inversions of both 35 s Love wave databases, we show in figs.7 (Harvard) and 8 (Utrecht) variance reduction as a function of norm and roughness damping parameter.We achieve a systematically better fit of Harvard data, with respect to Utrecht ones, regardless of the regularization.Consistently with the simple shape of the curves in fig.6, we find that variance reduction decreases monotonically with increasing damping parameters.
We conclude that the analysis of Section 3.2 holds true independently of the chosen damping parameters.

Conclusions
We confirm that a discrepancy exists between Ekström et al.'s (1997) and Trampert andWoodhouse's (1995, 2001) measurements of teleseismic surface wave dispersion.At long periods, it only concerns the short spatial frequency component of imaged phase velocity, while at short periods the discrepancy is significant at all harmonic degrees (figs. 1 and 2 above).
We have investigated the nature of this discrepancy, deriving independent, isotropic, global phase velocity maps from both databases, for both Love and Rayleigh waves at all available periods.
A preliminary set of inversions, with a parameterization of intermediate nominal resolution (5°×5°), shows that the Harvard data are fit systematically better than the Utrecht ones.This effect could be explained in various ways.The Harvard measurements might simply be cleaner, either because Ekström et al.'s (1997) algorithm is more accurate than Trampert andWoodhouse's (1995, 2001), or because the Harvard group started with a set of seismograms of higher quality.Alternatively, the Utrecht database might have higher sensitivity to small-scale structure, damped out from Ekström et al.'s (1997) measurements, and not resolved by a grid of 5°×5°pixels.
Having repeated our analysis with a much finer (2°×2°) parameterization, and with a wide spectrum of regularization schemes (figs.3 and  4; figs.6 through 8), we can rule out the latter speculation: the increase in nominal resolution has not affected the nature of the discrepancy between images based on Utrecht versus Harvard data, nor the systematic difference in data-fit.
We have, however, neglected anisotropy and scattering; while it is very likely that their effect be only minor, it is possible that models of phase velocity propagation including azimuthal anisotropy (Laske and Masters, 1998;Becker et al., 2003;Trampert and Woodhouse, 2003), and accounting for the effects of scattering on data sensitivity (e.g., Spetzler et al., 2002) are significantly different, in their isotropic component from the ones presented here.This issue calls for further work.
An additional explanation for the discrepancy in question could be a leakage of signal when dispersion is measured between neighboring surface wave periods.Figure 2 suggests that this effect might be severe enough, in Trampert and Woodhouse's (2003) algorithm, to cause the Moho anomalies of Tibet and Andes (to which short period Love wave should be strongly sensitive) to be missed.This speculation will need to be substantiated by benchmarking the two algorithms with each other.

Fig. 2 .
Fig. 2. Spherical harmonic spectra of the maps from fig. 1.The spectra of maps derived from Utrecht data are plotted in grey; the black curves are the spectra of the corresponding maps, derived from Harvard data.

Fig. 3 .
Fig. 3. Love wave dispersion data at 35 s, from both Trampert and Woodhouse's (2001) (left) and Ekström et al.'s (1997) (right) databases, are inverted with four different regularization schemes.Roughness and Norm damping parameters, always coincident in this experiment, are given to the left of each couple of maps.The bottom panel shows a theoretical 35 s Love phase velocity map: at each pixel of Crust-2.0, local mode eigenfunctions were used to calculate the phase velocity resulting from a spherically symmetric mantle and heterogeneous crust.We do not expect tomographic and theoretical maps to concide, but they should share some prominent features.The colour scale is the same as in fig. 1. 951

Fig. 4 .
Fig. 4. Harmonic spectra of the images in fig.3, plus (top panel) an obviously underdamped one, shown for comparison.Again, spectra associated with Utrecht and Harvard databases are in grey and black, respectively.

Fig. 5 .
Fig. 5. Correlation between maps as a function of period, forTrampert and Woodhouse's (2001) andEkström et al.'s (1997)  databases (solid and dashed lines, respectively).Subsets of the two databases at Rayleigh (top) and Love (bottom) wave periods approximately 35 s, 40 s, 50 s, 60 s, 75 s, 100 s and 150 s are inverted with the 5°×5°pixel parameterization of Section 3.1; the correlation of each of the resulting maps with the one, from the same database, at the neighbouring longer period is shown as a function of period.Correlation is computed according to Becker andBoschi (2002), and all harmonic degrees up to 40 are considered.For both Love and Rayleigh waves at short periods, maps from the Utrecht database are better correlated with each other than those from the Harvard one; the opposite is true at longer periods.

Fig. 8 .
Fig. 8. Variance reduction of 35 s Love wave dispersion from the Utrecht database, as a function of norm (x-axis) and roughness (y-axis) damping parameters.

Fig. 7 .
Fig. 7. Inversions of the 35 s Love wave Harvard database: variance reduction as a function of norm (x-axis) and roughness (y-axis) damping parameters.