Global reference data sets: surface waves 1809 from various techniques can be attributed to discrepant theoretical approximations, reference Earth models and processing schemes. Phase-velocity variations yielded by the inversion of the summary data set are highly correlated (R ≥ 0.8) with those from the quality-controlled contributing data sets. Long-wavelength variations in fundamental-mode dispersion (50-100 s) are largely independent of the measurement technique with high correlations extending up to degree ∼25. Agreement degrades with increasing branch number and period; highly correlated structure is found only up to degree ∼10 at longer periods (T > 150 s) and up to degree ∼8 for overtones. Only 2ζ azimuthal variations in phase velocity of fundamental-mode Rayleigh waves were required by the reference data set; maps of 2ζ azimuthal variations are highly consistent between catalogues ( R = 0.6-0.8). Reference data with uncertainties are useful for improving existing measurement techniques, validating models of interior structure, calculating teleseismic data corrections in local or multiscale investigations and developing a 3-D reference Earth model.

Mantle processes; Computational seismology; Seismic anisotropy; Seismic tomography; Surface waves and free oscillations.

A fundamental goal in seismology is to accurately determine the elastic structure of the Earth's interior. Elastic reference models are widely used in the geosciences for the modelling and interpretation of seismic sources and planetary interiors. Earth's interior has traditionally been described in terms of spherically symmetric (1-D) structure with physical properties varying radially within concentric shells such as the upper mantle and outer core. It is now well established that there are substantial three-dimensional (3-D) variations in the mantle and that 3-D tomographic models are useful as starting models for more detailed imaging studies and in constraining physical parameters such as temperature, grain size, fabric and composition. While radial (1-D) reference Earth models have been available for several decades (e.g. Dziewoński & Anderson 1981;Kennett et al. 1995), only recently has global seismic imaging reached a point where the development of a 3-D reference Earth model (REM3D) can be envisaged. A key component in this endeavor is to accurately characterize the arrival times of various phases observed on broad-band seismograms.

Surface waves are the most prominent phases recorded at teleseismic distances at periods longer than 30 s, especially from shallow-focus earthquakes. Two types of surface waves are observed, distinguished by their polarization during propagation through the Earth: Love (SH) and Rayleigh (P-SV) waves, recorded on the transverse and vertical/longitudinal components, respectively. Surface wave arrivals are denoted by the orbit number (e.g. N o = 1 for minor-arc L1 or R1 waves), a proxy for the number of times the wave circles around the Earth (N c = [N o -1]/2 for odd N o , N o /2 otherwise). The wave trains excited by large mega-thrust earthquakes (M w ≥ 7.5) circle the Earth multiple times (N c ≥ 1) for many hours and manifest as discernible higher-orbit arrivals (e.g. L3-L5, R3-R5). Generation and propagation of surface waves can also be classified based on the properties of the corresponding normal modes. Fundamental-mode surface wave trains are excited more strongly by shallow and intermediate-depth earthquakes (h < 250 km) and appear well separated from other phases at teleseismic distances ( > 30 • ). Higher-mode or overtone vibrations are excited by deeper earthquakes and appear as faster propagating, compact wave packets that contribute to the long-period body waveforms (e.g. Takeuchi & Saito 1972). Characterizing surface waves and overtones is critical for the construction of elastic reference Earth models.

In addition to their large amplitudes, surface waves are characterized by a frequency-dependence of velocity (i.e. dispersion). In cohort with other complementary data sets, laterally variable dispersion resulting from structural heterogeneity is crucial for mapping the upper mantle, transition zone and mid mantle (e.g. Masters et al. 2000;Ritsema et al. 2004;Moulik & Ekström 2014). Accounting for dispersion is also useful for locating earthquakes (e.g. Ekström 2006b;Howe et al. 2019), signal enhancement through phase-coherent stacking of seismograms (e.g. Ma et al. 2014), and inverting the centroid-moment tensors (CMTs) of seismic sources (Dziewoński et al. 1981;Ekström et al. 2005). Several techniques have been used to directly or indirectly measure dispersion of fundamental-mode surface waves and overtones (e.g. Dziewoński et al. 1972;Herrin & Goforth 1977;Lerner-Lam & Jordan 1983;Cara & Lévêque 1987;Stutzmann & Montagner 1993;Trampert & Woodhouse 1995;Ekström et al. 1997;van Heijst & Woodhouse 1997;Debayle 1999;Yoshizawa & Kennett 2002a;Beucler et al. 2003;Lebedev et al. 2005;Visser et al. 2007;Ma et al. 2014). These techniques use various processing and fitting procedures with different assumptions on crustal structure, mode coupling, reference model, geodetic parameters and attenuation. To date, no systematic assessment of the consistency in the resulting measurements has been performed. Such comparisons can help identify outliers or systematic biases and provide method-agnostic estimates of measurement uncertainty. Reference data sets with uncertainties are crucial for testing hypotheses about mantle structure, such as those concerning the depth and lateral variations of radial and azimuthal anisotropy (e.g. Trampert & Woodhouse 2003;Visser & Trampert 2008;Ekström 2011;Ma et al. 2014;Schaeffer et al. 2016). Additionally, inversions based on the reconciled reference data set can inform parametrization and regularization choices that strongly impact the inferences of mantle structure (e.g. Spetzler et al. 2002;Sieminski et al. 2004;van der Hilst & de Hoop 2005;Boschi et al. 2006;Trampert & Spetzler 2006). This is the first in a series of papers that describe a community effort to construct a 3-D reference model (REM3D) for the Earth's mantle. A major objective is to provide quality-controlled, comprehensive and publicly available seismological data sets with corresponding Downloaded from https://academic.oup.com/gji/article/228/3/1808/6408466 by University of California School of Law (Boalt Hall) user on 03 January 2022

The propagation of seismic surface waves was first developed in the framework of a flat, layered model of the Earth, to which corrections for sphericity are applied at far regional and teleseismic distances. Later, the equivalence between a propagating surface wave formalism and a normal mode formalism in a spherically symmetric Earth model, was established (e.g. Gilbert 1976;Aki & Richards 1980). Given the ensemble of Rayleigh (or Love) wave trains propagating along a great circle path, we denote the successive surface wave trains propagating in the direction of the minor arc from the source to the receiver as R1, R3, R5 (or L1, L3, L5) and those propagating in the opposite direction as R2, R4, R6 (or L2, L4, L6). Surface wave trains can be interpreted in terms of Rayleigh-wave equivalent spheroidal modes n S l or Love-wave Downloaded from https://academic.oup.com/gji/article/228/3/1808/6408466 by University of California School of Law (Boalt Hall) user on 03 January 2022 Table 1. Summary of the raw data used in this study. Measurements comprise either eigenfrequency perturbations, phase velocity curves or phase anomalies (in seconds) at various discrepant ranges of periods (Section 4.4). Both minor-arc (L1, R1) and major-arc data (L2, R2) are analysed from two sources (GDM52, MBS11), while higher-orbit measurements (L3-L5, R3-R5) are included from GDM52. Mode branches represent either the fundamental-mode measurements (n = 0) or those of the higher overtones (n = 1-18). N contrib is the number of raw dispersion measurements at discrete frequencies reported by the authors.

Also provided are the percentages of raw measurements whose station and source metadata can be retrieved or cross-validated (Section 4.4). N raw is the number of raw measurements after interpolation to reference frequencies, which is typically less than N contrib for catalogues that over-sample the dispersion curves but can be higher when limited frequencies are contributed (e.g. Cambridge19). N grid is a measure of the spatial sampling and refers to the number of knot pairs on an evenly spaced grid of 2562 points where the summary data from the catalogue is available. N eq , N st and N path are the number of earthquakes, stations and paths in N contrib , respectively. A total of N contrib = ∼227 million (49.65 million Love, 177.66 million Rayleigh) dispersion measurements reported at discrete frequencies were analysed towards the reference data set. Source equivalent toroidal modes n T l with radial order n and angular order l. For a particular mode type (spheroidal S or toroidal T), the displacement time series recorded at the receiver r can be expressed as a sum over normal modes as follows:

where n ω l is the complex eigenfrequency of the mode, and n A l is its excitation amplitude for the particular source-station configuration.

Alternatively, the displacement time series can be expressed, in the high frequency limit, as a sum of propagating surface waves as:

where A n (ω) and n (ω) are the amplitude and phase of the nth surface wave overtone as a function of angular frequency ω (Aki & Richards 1980).

Dropping the overtone index n, the phase for a particular source receiver pair comprises four contributions in a spherically symmetric Earth model

where S is the contribution from the source, R is the receiver phase shift, P is the contribution to the phase due to propagation from the source to the receiver and M is a signed integer, which represents the number of passages through the source antipode (e.g. M = 0 for R1 and L1, M = 1 for R3 and L3, M = -1 for R2 and L2). Similarly, the amplitude term A n (ω) can be decomposed into a product of contributions: Flow chart of various processing steps in this study. Overall, seven major steps are adopted in the construction of the reference data set resulting in contributed, raw, homogenized, clean and summary catalogues. The corresponding section in the manuscript where the step is discussed is specified in bold red. RSDF files in both ASCII and HDF5 formats (Table A1) store the measurements on homogenized and original paths.

where A S is the magnitude of the excitation at the source, A R is the receiver amplification, A F is the geometrical spreading factor and A Q is the decay factor due to anelastic attenuation along the ray path. The propagation phase is defined in terms of a phase velocity C(ω) as where X is the distance travelled by the wave. When we assume propagation to follow the great circle path, X equals , the great-circle distance that is calculated using a distance factor ( F = 111.31948) after applying the geocentric conversion factor (W = 0.9933056) to the locations in geographic coordinates (e.g. Seidelmann 1992;Moulik & Ekström 2021). In a spherically symmetric Earth model, C(ω) does not depend on the source-station geometry. For a given mode branch n, phase velocity is related to the corresponding mode eigenfrequency as

where R = 6371 km is the mean radius of the Earth (Jeans 1927).

In the 3-D Earth, the phase velocity measured on a given source-station path depends on the location of the source and station, and on the source radiation pattern to account for potential off-great circle propagation. A common assumption made in interpreting phase velocity measurements is the 'path average' approximation (PAVA, Woodhouse & Dziewoński 1984). Propagation is assumed to occur along the great circle path, and the propagating phase only depends on the distance between the source and receiver. Given a reference spherically symmetric Earth model, in which the phase velocity for a given surface wave branch is C 0 (ω) [or its inverse, the phase slowness P 0 (ω)], the propagation phase in the reference model is:

and the measured propagation phase P between two points distant by can then be written as

where δC and δP are the perturbations in phase velocity and slowness, respectively, due to variations in velocity away from the reference spherically symmetric model. The integer S accounts for indeterminacy due to the definition of phase modulo 2π . The varying structure along the path s is most conveniently described by 'local' phase slowness perturbations δp(ω, s), such that

Note that care must be taken to avoid 'cycle skipping' when inferring the slowness perturbation δP. Since the differences in dispersion between the predictions from a reference 1-D model and the real observations are small at long periods (>100 s), there is usually no ambiguity in the selection of S. Most surface wave dispersion studies start processing at longer periods, so that continuous dispersion curves can be anchored, and the total phase perturbation at shorter periods can be inferred with less ambiguity. Of particular interest to this study is the distribution of local phase slowness δp(ω, s) and its inverse, phase velocity δc(ω, s) at a given frequency and mode branch. Such two-dimensional (2-D) maps can be derived by the inversion of measured phase slowness perturbations δP(ω) over many source-station paths, while potentially including measurements on higher orbits. The resulting phase dispersion curves obtained over a band of frequencies at each point on the Earth's surface can be inverted in turn for elastic structure as a function of depth, using sensitivity kernels derived from normal mode perturbation theory or fully numerical approaches. The perturbation in phase velocity at a fixed frequency ω is related to the perturbation in local eigenfrequency at a fixed wavenumber k as

where U is the group velocity (e.g. Dahlen & Tromp 1998). While it is straightforward to derive eq. ( 9) for surface waves in the frequency domain, relating it to normal mode perturbation theory took some theoretical development. Simply perturbing the eigenfrequency of a mode only allows us to represent the effect of heterogeneity integrated over the entire great circle path (e.g. Jordan 1978), and therefore sensitivity to structure that is symmetric with respect to the centre of the Earth ('even order' heterogeneity). It can be shown that the PAVA approximation for surface waves is equivalent to along-branch mode coupling in the asymptotic limit of large angular orders of first order perturbation theory applied to normal modes (Mochizuki 1986;Park 1987;Romanowicz 1987). The along-branch coupling brings out the sensitivity of the modes to odd-order heterogeneity. Most surface wave and overtone phase dispersion measurement techniques implicitly utilize the PAVA approximation to relate 3-D structural heterogeneity at depth to observed slownesses. In the rest of the paper, we will use the greatcircle ray approximation (GCRA), which is a related infinite-frequency approximation that predicts phase delays from 2-D slowness maps without accounting for finite-frequency (e.g. FFT, Wang & Dahlen 1995b;Yoshizawa & Kennett 2002b;Zhou et al. 2004) or off-great-circle propagation effects adopted in exact ray theory (e.g. ERT, Woodhouse & Wong 1986;Larson et al. 1998). Similar structures can be obtained using GCRA, FFT and ERT theory depending on the choices of parametrization and regularization (Spetzler et al. 2002;Sieminski et al. 2004;Boschi et al. 2006;Trampert & Spetzler 2006). Based on synthetic experiments, GCRA accurately predicts phase anomalies of minor-arc phases and matches input phase slowness maps at global scales (e.g. Godfrey et al. 2019). The basic assumption in GCRA that rays travel along the great circle connecting the source and receiver may become less valid with increasing path length (e.g. Woodhouse & Wong 1986), such as in the case of higher-orbit measurements (L3-L5, R3-R5). A detailed comparison of theoretical assumptions for all wave types is beyond the scope of this study. However, we note that Wang & Dahlen (1995b) found little dependence of errors in the ERT approximation on the orbit number of surface waves. (Lebedev et al. 2005;Schaeffer & Lebedev 2013), GDM52 (Ekström et al. 1997;Ekström 2011), IPGP03 (Stutzmann & Montagner 1994;Beucler et al. 2003;Beucler & Montagner 2006), Lyon18 (Debayle 1999;Debayle & Ricard 2012), MBS11 (van Heijst & Woodhouse 1997;Ritsema et al. 2011), Scripps14 (Ma et al. 2014) and Utrecht08 (Yoshizawa & Kennett 2002a;Visser et al. 2007). Fundamental-mode (n = 0), minor-arc (L1, R1) measurements were the most common type of data across the eight contributions. Additional constraints on major-arc arrivals (L2, R2) were available from two sources (GDM52, MBS11), while higher-orbit measurements (L3-L5, R3-R5) were included from GDM52. All groups contributed measurements in terms of path-dependent dispersion curves for various overtone branches (n = 0-18) sampled unevenly at different sets of discrete frequencies. The contributions included measurements from recent analyses and unpublished updates in formats that accounted for the processing guidelines in this study (Section 4.3, the Appendix).

Rayleigh-wave dispersion data on the vertical component are more widely available (>3 times) than Love-wave measurements due to the inherently noisier records on the horizontal component seismograms. The contributed compilation represents the largest and most diverse set of surface wave arrival times assembled to date. Fig. 1 shows the reported locations of 44 871 sources and 12 222 receivers that contributed at least one observation to this study. Waveform data for the majority of catalogues were available from the Incorporated Research Institutions for Seismology (IRIS). All catalogues adopted in their measurement procedure the source mechanisms (CMTs) from the Global CMT project (Dziewoński et al. 1981;Ekström et al. 2005) in their measurement procedure. Recent implementations of the Global CMT algorithm minimizes the difference between observed and synthetic seismograms in three frequency bands and time windows. These include the body waves (40-150 s), long-period mantle waves (125-350 s) and surface waves with bandpass varying with event size (50-150 s for M W = 6). After salvaging and validating metadata (Section 4.4), our compilation included 40 122 earthquakes (moment magnitude, M W = 4.6-9.1) between 1976 and 2016 recorded on 10 469 stations and 310 networks accessible through the open IRIS data centre. The measurements were made on seismograms recorded on globally distributed permanent stations as well as temporary deployments. Some common permanent stations included the Global Seismographic Network (network codes II and IU), the Chinese Digital Seismograph Network (CD and IC), the Mednet (MN), Geoscope (G), Geofon (GE) and Caribbean (CU) Networks, the Global Telemetered Seismograph Network (GT), Brazilian Lithospheric Seismic Project (BL), United States National Seismic Network (US), Southern California Seismic Network (CI) and selected stations of the Canadian National Seismograph Network (CN). Temporary deployments included the Hawaiian PLUME experiment (ZF), the POLARIS array in northern Canada, Earthscope USArray transportable array (TA, 1693 locations), SKIPPY array in Australia (7B) and those of the United States Geological Survey (GS).

Fig. 3 shows the ray coverage of fundamental-mode Rayleigh waves (R1) at 100 s from various catalogues. Hit count is defined as the number of rays traversing every 2 • pixel, normalized by relative area to account for smaller pixels at higher latitudes. Global averages of hit counts for these waves were the highest (>3000) for the Cambridge19, MBS11, Dublin13, and Lyon18 catalogues. The inclusion of temporary PASSCAL array deployments helped improve hit counts in the Pacific Ocean Basin and Southern Hemisphere, especially in Africa, Antarctica and South America. Nevertheless, large areas in the southern oceans still lack good station coverage and hit counts differ laterally by up to 3-4 orders of magnitude. Several catalogues provide uneven coverage by repeatedly sampling paths from the numerous earthquakes in the Tonga-Kermadec subduction zone to the large number of stations in North America.

Several pioneering efforts since the 1960s have led to sophisticated techniques for measuring surface wave dispersion. In the interest of brevity, we will discuss only those aspects of measurement techniques that are relevant for data reconciliation. Several procedures are common across various techniques for measuring surface wave dispersion. For example, Rayleigh wave dispersion is determined from the vertical component seismograms and Love wave dispersion from transverse seismograms after rotation of the horizontal components using the great-circle back azimuth. Most measurement techniques compare synthetic and observed seismograms either in the frequency (IPGP03) or time domain (Cambridge 19, Dublin13, GDM52, Lyon18, MBS11, Utrecht08) while Scripps14 compares pairs of observed seismograms. Dispersion and amplitude of the synthetic waveform are adjusted to minimize the residual dispersion and the associated misfit between seismograms. The end product of interest is a smoothly varying perturbation in apparent phase velocity c 0 + δc valid for a range of periods, as well as parameters quantifying the quality of the measurement typically based on measures of waveform fit. Most dispersion measurement techniques proceed one record at a time using semi-automated schemes that use filter and processing criteria informed by domain experts.

While the overarching goal of the techniques are similar, details of the processing scheme can lead to inconsistencies in the measurements and inferences on Earth structure. Techniques for measuring surface wave dispersion can differ in their choices of: (1) fundamental-mode only versus multimode schemes; (2) methods for computing synthetic waveforms, and, when necessary, sensitivity kernels; (3) data processing Downloaded from https://academic.oup.com/gji/article/228/3/1808/6408466 by University of California School of Law (Boalt Hall) user on 03 January 2022 choices such as misfit criteria, windowing, filtering, and the use of cross-correlation; (4) framework for interpolation or parametrization across frequencies; (5) extent of automation and (6) criteria for selecting sources and stations.

Many techniques are designed to either process fundamental modes and overtones jointly in multimode schemes or consider fundamental modes in isolation. Measurement of fundamental mode phase velocities are considered relatively more straightforward if certain data processing criteria are adopted. Two of the contributed data sets, GDM52 and Scripps14, are fundamental mode-only catalogues, and restrict their analysis to shallow earthquakes (h < 50-250 km) that excite strongly fundamental mode surface waves and ensure that these wave trains are the dominant long-period phase in the seismograms. GDM52 measures dispersion using synthetic seismograms that do not account for the Downloaded from https://academic.oup.com/gji/article/228/3/1808/6408466 by University of California School of Law (Boalt Hall) user on 03 January 2022 contribution from overtones. Fundamental-mode energy is isolated by suppressing the contributions from the interfering overtones based on ideas from residual dispersion (e.g. Dziewoński et al. 1972), phase-matched filtering (e.g. Herrin & Goforth 1977) and optimally windowing the cross-correlation function (e.g. Ekström et al. 1997). Scripps14 uses dispersion predictions from GDM52 to calculate an 'undisperse' term that helps with aligning observed seismograms for clustering, especially at frequencies higher than 20 mHz where the procedure is more susceptible to cycle skipping.

The extraction of overtone information from surface wave seismograms is an underdetermined problem due to the similar group velocities and associated simultaneous arrivals of various branches. This process has a much wider range of quasi-linearity for Rayleigh waves than for Love waves, due to the clear separation of the fundamental mode. The choice of the starting 1-D or 3-D model for guiding dispersion measurements could therefore be more significant for Love waves and influence the results. Nevertheless, several multimode studies have been developed to extract the overtone signal in the data. Dublin13 and Utrecht08 determine time-frequency windows in which synthetic seismograms fit the data closely, identify the modes that contribute significantly to these waveforms, and measure their phase velocities. Cambridge19 and Lyon18 cross-correlate the complete observed waveform with pure-mode synthetics for different overtone branches to monitor along-branch dispersion and residual fits to the observed cross-correlograms. MBS11 uses an iterative mode-branch stripping technique in which the cross-correlogram between the observed waveform and the single most energetic mode branch is fit to determine phase velocity and amplitude perturbations, and the waveforms predicted for that branch are iteratively subtracted from the observed waveform. IPGP03 uses non-linear optimization to fit dispersion curves simultaneously to groups of waveforms from multiple nearby sources with different depths and source mechanisms to potentially make the extraction of overtone information less underdetermined.

A common source of discrepancy among surface wave studies lies in the theoretical procedure for calculating synthetic predictions. These could either involve corrections for undispersed waveforms to enable stable cross-correlation comparisons (Scripps14), or synthetic waveforms for comparison with observations in other catalogues. In most dispersion studies, synthetics are initially computed in a reference spherically symmetric (1-D) Earth model, though different choices of both the elastic (e.g. isotropic vs. anisotropic PREM) and anelastic structure in the reference models are common. Cambridge19 and Lyon18 use path-specific reference 1-D models that capture the average crustal structure along each path, while Dublin13 uses reference phase velocities c 0 (ω) that account for off-great-circle-path sensitivity in a 3-D crustal model. The non-linear optimization scheme used in IPGP03 could make the resulting phase measurements insensitive to the reference model used. These choices can affect the reference propagation phase ( 0 P , eq. 7) systematically with great-circle distance ( ), either directly through different reference phase velocities c 0 (ω) or through anelastic dispersion with frequency (Kanamori & Anderson 1977).

The differences in data processing across various techniques may be grouped into two categories. First, the techniques differ in the way the misfit is evaluated. In Cambridge19, GDM52, Lyon18 and MBS11, misfit is calculated on the cross-correlograms between observed and synthetic seismograms, which highlights sensitivity to a particular branch (Lerner-Lam & Jordan 1983) and enables precise measurements of dispersion (Dziewoński et al. 1972). Scripps14 cross-correlates pairs of observed seismograms to measure relative traveltime differences. Dublin13 and Utrecht08 calculate the misfit in the time domain within multiple time-frequency windows (Yoshizawa & Kennett 2002a;Lebedev et al. 2005). Secondly, a major difference among the techniques pertains to the construction of dispersion curves. Some groups minimize misfit between waveforms by parametrizing smoothly varying dispersion curves in terms of spline coefficients (GDM52, MBS11, Scripps14) or by imposing smoothness through a priori covariance (IPGP03). Alternatively, path-average 1-D models that are perturbations to a global or regionalized reference 1-D model are inverted using the PAVA approximation. These path-average models are then used to predict the dispersion curves for branches and frequencies that contribute substantially to the misfit (Cambridge19, Lyon18, Utrecht08).

While most dispersion data sets provide good geographic coverage due to the proliferation of seismographic networks, details of the measurement technique can place limitations on the number of available paths. In order to obtain reliable multimode dispersion measurements, IPGP03 requires waveforms from multiple nearby sources, which somewhat limits the geographic coverage of that data set. Since most surface wave techniques evaluate a single record at a time, various subjective criteria are used to quality control the data, automate the processing scheme and estimate uncertainty. However, IPGP03 quantifies uncertainty on phase dispersion parameters from the simultaneous inversion of waveforms from multiple, nearby sources, and Utrecht08 samples the full probability density function. Spurious measurements may be due to instrument polarity reversals, response function errors, timing problems, and dead channels (e.g. Ekström et al. 2007). Scripps14 and Dublin13 account explicitly for polarity reversals on the current Global Seismographic Network (GSN) based on a manual list of known issues. If unaccounted for, these polarity reversals can cause a half-cycle ambiguity (π ) in an isolated residual phase measurement. Dublin13 uses outlier analysis and removes from the data set the least mutually consistent measurements, which are likely to be contaminated by instrumental and event-location errors (Lebedev & Van Der Hilst 2008;Schaeffer & Lebedev 2013).

Our basic observation is the arrival time of a dispersed surface wave at a broad-band seismometer from an earthquake source. Contributed dispersion curves are typically provided either as a propagation phase anomaly (δ ; GDM52, Scripps14) or the inferred average phasevelocity perturbation ( δc; IPGP03, Dublin13, Lyon18, Cambridge19). Other studies (i.e. MBS11, Utrecht08) report fractional perturbation in eigenfrequency (δω/ω 0 ) to the normal mode nearest to the frequency of interest. These choices of how measurements are tabulated are associated with differences in measurement techniques (Section 4.2). Eigenfrequency perturbations are common in waveform approaches where large number of differential waveforms need to be evaluated (e.g. MBS11). Propagation phase anomalies are easily retrieved with Downloaded from https://academic.oup.com/gji/article/228/3/1808/6408466 by University of California School of Law (Boalt Hall) user on 03 January 2022 cross-correlation of narrow-band seismograms (e.g. GDM52) while phase-velocity perturbations are reported when inversion of a pathaverage 1-D model is part of the processing (e.g. Cambridge19, Lyon18). All contributed measurements are converted to propagation phase anomalies (δ , in seconds) of a surface wave component (e.g. R1) relative to a reference phase from the reference model ( 0P ). We will discuss propagation phase in either seconds or radians interchangeably in the rest of the paper.

The contributed source and station locations provided in geocentric coordinates are converted to geographic coordinates. Reference phase ( 0 P ) is calculated based on the radial reference Earth model reported in the study and the reported great-circle distance (eq. 7). In case of catalogues that report eigenfrequency perturbations, measurements are converted to propagation phase anomalies following eq. ( 10) as

We account for the discrepant values of the geodetic constants used in contributed data sets during these conversions whenever available (e.g. 111.1949 in GDM52). Due to the use of different geodetic constants, reported reference phases ( 0 ) from various catalogues have baseline differences that lead to discrepancies in phase anomalies (δ ). Geodetic contribution to the discrepancies is typically small, around 3-5 s for minor-arc Rayleigh (R1) waves at a period of 150 s. The uncertainties of propagation phase anomalies are also converted to seconds while preserving the relative uncertainty in reported data. All contributed data are stored in ASCII versions of reference seismic data formats (RSDF, the Appendix), where the columns represent measurements while metadata and original processing notes are preserved as headers (Table A1 and Fig. 2).

The availability of all relevant metadata is a critical requirement for reconciling contributed catalogues. Several contributions had missing or incomplete source and station information that needed to be salvaged or cross-validated against relevant sources. A persistent issue with the contributed catalogues was the lack of earthquake source information. When moment tensors from the Global CMT catalogue were used, appropriate indexing that would facilitate cross-validation (e.g. cmtname from Table A1) was sometimes not preserved. When the event names were provided, they were arbitrarily defined (e.g. custom timestamps) and did not correspond to those in the Global CMT catalogue. Some catalogues only provided the epicentre coordinates with no corresponding source depth or centroid time information. In active source regions with several hundreds of M W > 5.5 earthquakes every year, it became impossible to easily attribute the measurement to the correct CMT source mechanism. Since conservative quality control criteria were used in this study that could potentially exclude substantial portions of the catalogues (Section 5.3), we implemented a standard procedure to retrieve as much of the missing metadata as possible. The procedure for retrieving the source information was informed by: (1) range of source magnitudes and depths; (2) year of the study and duration of data analysed, if provided; (3) reported epicentre information and (4) date of the earthquake, if provided. After filtering the Global CMT catalogue based on criteria 1 and 2, the source nearest to the reported hypocentre ( = 0.01 • , depth h = 1 km) was found. If a unique source was not retrieved, an additional search was performed based on available date information (4) often codified as a timestamp in the contributed catalogues. We were able to cross-validate a substantial portion of the source mechanisms for all catalogues (73-99 per cent). Note that this procedure was not applied to the IPGP03 catalogue, whose measurements refer to source regions rather than specific earthquakes (Section 4.2). The most complete source metadata (≥99 per cent) were found for the GDM52, Dublin13 and Lyon18 catalogues.

A standard approach was also adopted to retrieve and validate missing station metadata against published databases. For every measurement, all stations active on the day of the CMT source event were filtered from the IRIS database. One of two procedures was selected based on the type of reported metadata. If no network and station code were available, stations within a threshold great-circle distance ( = 0.01 • ) were identified. We cross-validated reported network and station codes against any available codes if no station coordinates were available. In case of conflicts between network codes for the same station, we preferred IRIS network codes in a prescribed order (IU, II, CD, IC, MN, G, GE, CU, GT and CN). Location codes were preserved only when reported by the catalogue (e.g. MBS11), and no attempts were made to identify these during processing. These steps were repeated until a unique station code was found, which sometimes required manual intervention. A majority of the measurements in all catalogues cleared the metadata analysis for both sources and stations (Table 1). More than 96 per cent contributed measurements were cross-validated in several catalogues that preserve detailed information on their processing schemes (e.g. GDM52, Lyon18, Scripps14 and Dublin13). For the MBS11 catalogue, minor-arc measurements cleared both analyses at a substantially higher rate (≥99 per cent) than major-arc measurements. Almost all the source metadata for the Cambridge19 catalogue were found, but only 83-89 per cent of the stations could be validated with our choice to restrict analyses to stations accessible through to the open IRIS data centre. A substantial portion of Cambridge19 measurements are from stations whose waveforms are either available from other open repositories (e.g. European Integrated Data Archive, EIDA) or closed networks, thereby limiting their utilization towards this reference data set. Since Europe is already represented well by stations from the IRIS networks, the loss of information is not severe for this continent. After the retrieval of metadata, the great-circle distance was re-calculated and the phase anomalies were updated to account for any changes to the distance and the related reference propagation phase ( 0 P ). Differences between the reported and calculated distances are typically small (<0.003 per cent) but can lead to discrepancies of a few seconds for higher-orbit waves (R3-R5 and L3-L5).

Downloaded from https://academic.oup.com/gji/article/228/3/1808/6408466 by University of California School of Law (Boalt Hall) user on 03 January 2022 Table 2. Types of filter criteria used in the quality control and outlier removal to construct the clean data set. Quality refers to the filters applied to account for data availability in distance ranges with strong intra-catalogue consistency. Interference filters are from bands and distances where substantial discrepancies are seen owing to overtone effects. Source filters account for the expected strong excitation of the waves. In addition, threshold filters for intra-catalogue deviations and inter-catalogue discrepancies are also used (Section 5.3). Surface wave dispersion is reported across a wide variety of frequency ranges and sampling. Cambridge19 reports the dispersion curves coarsely sampled in frequency while Dublin13 has the finest sampling. For each catalogue, we calculated dispersion curves for every sourcereceiver path and overtone branch at a discrete set of reference periods roughly equally spaced in frequency (25s, 27s, 30s, 32s, 35s, 40s, 45s, 50s, 60s, 75s, 100s, 125s, 150s, 175s, 200s and 250s) using cubic spline interpolation. The implicit assumption of smoothly varying phase within the same nth-overtone branch is physically justified due to similarities in the corresponding depth sensitivities to radial structure. The interpolation procedure accounts for the intersection of Stoneley wave and core-mode branches with spheroidal overtones (e.g. Okal 1978;Dahlen & Tromp 1998); constant n therefore corresponds to a smooth overtone branch in which modes with neighbouring l have similar physical characteristics. However, the retrieval of smooth dispersion curves from real data can be made infeasible at near-nodal take-off angles and along paths that generate multipathing. We assumed that the various measurement techniques naturally exclude paths with such complications since they tend to provide poor fits with synthetic waveforms (e.g. Ekström 2011).

A large variation in the details of dispersion curves are noted in the contributed catalogues. To ensure reliable interpolation of the raw catalogues, we only included paths with dispersion measurements reported for at least three discrete frequencies. The minimum number of dispersion measurements was reduced to two for Cambridge19 data since only a narrow band of frequencies is available for higher overtone branches. For all catalogues except IPGP03, only measurements that have cross-validated earthquake sources and stations (Section 4.4) were included. Table 1 provides the resulting number of raw measurements (N raw ). The number of dispersion measurements was reduced substantially through this procedure; three catalogues with the largest number of resulting measurements are Cambridge19, MBS11 and Lyon18.

The raw data were then quality-controlled based on various criteria that facilitate inter-catalogue comparisons (Table 2). For the initial analysis, we selected earthquakes M W ≥ 5.5 that were likely to excite the relevant intermediate-period surface waves for various overtone branches. Shallow sources were used for fundamental modes (h = 0-250 km) while deeper sources are permitted for overtone data (h = 0-650 km). The analysis was restricted to period ranges where at least two catalogues provide independent constraints. Observations between 25 and 200 s were analysed for fundamental modes and narrower bands were considered for higher overtones (e.g. 40-50 s for the 6th overtone). In addition, we considered paths in the teleseismic distance range (30 • ≤ ≤ 150 • ) in order to avoid complexities at short distances and near the antipode. There are other methodological reasons to exclude measurements based on epicentral distance. Clustering of different events in IPGP03 render the average measurements along common ray paths unsuitable at short epicentral distances ( ≤ 55 • ), where discrepancies in the path-specific corrections for various events become comparable in size to the signal. The mode-branch stripping technique is more effective on longer paths where there is lesser overlap in the arrivals of higher-mode branches (van Heijst & Woodhouse 1999). The raw data obtained using these selection criteria were stored in HDF5 RSDF files to facilitate rapid processing and inter-catalogue comparisons (Fig. 2).

Fig. 4 shows the root-mean-square (RMS) strength of the phase anomalies in the raw catalogue after subtracting a global average phase-velocity contribution at each frequency. For fundamental modes, the phase-anomaly RMS variations increase roughly as the square of the frequency and reach up to 3 full cycles (6π ) for both Love and Rayleigh waves between 25 and 35 s in the GDM52 and Scripps14 catalogues. This trend can largely be explained by the greater sensitivity of higher frequency waves to the strong heterogeneity in the crust and upper mantle (i.e. the heterosphere, Dziewonski et al. 2010). RMS variations of fundamental-mode Love waves are substantially higher (by 0.5-0.6 wavelengths) than those of Rayleigh waves, especially at frequencies higher than 20 mHz, potentially due to the shallower sensitivity of Love waves at these periods (e.g. Takeuchi & Saito 1972). Other raw catalogues such as Dublin13 that provide data at these frequencies also show similar but less dramatic trends with frequency. The trends observed for fundamental-mode data result from the sensitivity of high-frequency waves to strong lateral variations in crustal thickness and velocities in the heterosphere. Nevertheless, some differences across phase anomalies in the raw catalogues can potentially be attributed to the details of the measurement techniques rather than Earth structure. For example, RMS variations for Dublin13 measurements (dark cyan symbols) are substantially lower than GDM52 and Scripps by up to 2.5 cycles at frequencies higher than 25 mHz. Large residuals in Dublin13 measurements may get removed as outliers during a conservative Downloaded from https://academic.oup.com/gji/article/228/3/1808/6408466 by University of California School of Law (Boalt Hall) user on 03 January 2022 Values for raw data are provided as symbols while those for clean data are specified as lines (refer Fig. 2). Note the improvement in consistency between RMS strength of fundamental-mode dispersion at periods shorter than 35 s from removal of outliers (Section 5.3). Downloaded from https://academic.oup.com/gji/article/228/3/1808/6408466 by University of California School of Law (Boalt Hall) user on 03 January 2022 Figure 5. Scatter density plot of 66 787 raw propagation phase anomaly measurements for 100 s R1 waves common to both GDM52 (δφ 1 ) and Lyon18 catalogues (δφ 2 ). Scatter points are coloured based on the spatial density of nearby points, with high density in red and low density in blue. Histogram of the differences between the reported measurements (δφ 2 -δφ 1 ) on a logarithmic scale. Note the full-cycle (2π ) band of discrepancies between the catalogues for a small subset of common paths (∼1000) that also manifest as minor peaks in the histogram. GDM52 reports slower velocities with arrival times that are 1.4 s longer on average than Lyon18. Both mean and median of the absolute differences in measurements binned every 2 • show a linear increase with great-circle distance between the source and receiver. Such minor discrepancies may arise from discrepant theoretical approximations, reference Earth models and processing schemes.

analysis procedure that selects only the most mutually consistent data based on the final tomographic model (e.g. Lebedev & Van Der Hilst 2008;Schaeffer & Lebedev 2013). For the available overtone data, the RMS variation increases with frequency but never exceeds ∼0.6 cycles. Dublin13 and Utrecht08 are the catalogues with the lowest RMS variations. No clear and systematic differences in RMS strengths are seen between Love and Rayleigh wave overtones. Lower RMS strength in overtone data is likely due to a peak in sensitivity in the transition zone and mid mantle, where strength of heterogeneity is known to be weaker than in the uppermost mantle.

The consistency across catalogues can be evaluated by comparing directly the raw phase anomalies at reference periods on identical paths. Successful metadata retrieval of a vast majority of measurements (73-100 per cent, Table 1) permits direct comparison on reported sourcestation paths. Overall, more than 95 per cent raw measurements of fundamental-mode dispersion are consistent across catalogues and can be readily reconciled. Fig. 5 shows the graphical output generated to summarize comparisons for 100 s minor-arc Rayleigh wave measurements from the GDM52 and Lyon18 catalogues; similar analysis is conducted for every pair of catalogues with overlapping constraints. Most of the 66,787 raw measurements common to both catalogues fall near the 1-to-1 line on the scatter density plot, reflecting a high level of inter-catalogue consistency. However, there are a few common paths (∼1000) on which the discrepancies between catalogues are large (≥ π /4 radians). Most of these paths correspond to full-cycle differences in phase (±2π), which can be a cycle-skipping issue in either catalogue. Moreover, GDM52 reports slower velocities with arrival times that are 1.4 s longer on average than Lyon18. Both mean and median absolute differences between the two catalogues binned every 2 • show a steady increase with great-circle distance. These variations may be due to methodological assumptions, such as differences in the dispersion correction, geodetic constants or attenuation in the reference Earth model. Figs 6 and S2-S13 summarize the catalogue scatter comparisons for minor-and major-arc waves at 50, 150 and 250 s. Mean differences in minor-arc data are usually low (< 2.54 s for R1 and L1 at 50 s) and rarely exceed π /4 radians for data between 25 and 250 s. Constraints on major-arc data afforded by MBS11 and GDM52 catalogues are also highly consistent with very low mean differences (<1.72 s) for 150-250 s data. While the level of agreement between catalogues is generally high, some inconsistencies are seen across a wide band of frequencies. Full-cycle differences likely related to cycle-skipping issues are observed for 50 s waves between several pairs of catalogues (e.g. Lyon18-GDM52) and are more evident in minor-arc Rayleigh-wave data. Cycle skipping problems are more acute at shorter periods where it is harder to resolve the ambiguity in the number of cycles (Section 3). Such issues are also noted for station pairs at large epicentral distances along strongly heterogeneous paths for which the accrued phase delay may approach or exceed the period. The use of GDM52 as the starting model in Scripps14 helps alleviate some of the cycle-skipping issues (Section 4.2), producing greater overall consistency between the two catalogues. Much of the scatter in several catalogue-pairs are from measurements at epicentral distances outside the distance range 30 • -150 • used in the construction of the reference data set.

In the interest of brevity, we summarize the agreement for all types of fundamental-mode measurements in Fig. 7. For every pair of catalogues, median absolute deviations in reported measurements of Love and Rayleigh waves are plotted against reference period. Median differences in fundamental-mode Love waves are uniformly low (<5 s) for all combinations of catalogues except for Cambridge19 where discrepancies with GDM52, MBS11 and Scripps14 can exceed 6 s at the longest periods (>150 s). In contrast, consistency deteriorates for Rayleigh waves with median absolute differences exceeding 6 s for some pairs of catalogues (e.g. Cambridge19-Scripps14) at the longest periods (>200 s). Dublin13 and Utrecht08 have the most consistent Love-and Rayleigh-wave measurements with median differences not exceeding 1.5 s across 40-150 s period band. This reflects the central role that automated multimode inversion (Lebedev et al. 2005) Kustowski et al. 2008), and from oceanic and continental profiles constructed using the GTR1 tectonic regionalization (Jordan 1981). All synthetic waveforms are narrow-pass filtered between 8 and 12 mHz. both catalogues (Visser & Trampert 2008). Median values for Rayleigh waves binned every 2 • in epicentral distance show a clear deterioration in agreement outside the distance range 30-150 • .

Despite the low median differences between raw catalogues, our comparisons reveal a systematic trend in the deviations of fundamental mode Love waves. Fig. 8(a) shows that discrepancies between data sets that explicitly account for contamination by overtones (Cambridge 16, Dublin13, MBS11 and Utrecht08), and those that do not (GDM52 and Scripps14), oscillate with epicentral distance, peaking at ∼20 • , 45 • , 65 • , 90 • , 110 • and 130 • . This effect is most prominent for Love waves between 60 and 125 s period with slight indications at 150 s, albeit at somewhat longer distances. Notably, phase delay discrepancies between catalogues within each group do not exhibit clear epicentral distance trends. No such periodicity is observed in discrepancies of Rayleigh wave data (Fig. 7). A potential source of this discrepancy between measurement techniques may be due to overtone interference, which is known to influence the measurements of fundamental-mode surface waves (e.g. Foster et al. 2014b;Hariharan et al. 2020). We investigated the effect of overtone interference on fundamental-mode Love waves with synthetic seismograms computed for the 1-D Earth model STW105 (Kustowski et al. 2008). We also used slightly modified versions of STW105 with oceanic and continental-type structures from the global tectonic regionalization (e.g. GTR1, Jordan 1981). For each model, we computed three sets of synthetics: fundamental mode only, fundamental mode and first overtone, and all mode branches that contribute to our band of interest (8-12 mHz). We then applied a narrow bandpass filter around a target frequency, computed the amplitude ratio and the difference in the instantaneous phase between pairs of synthetics at the predicted arrival time of fundamental-mode energy for the target frequency. Figs 8(b) and (c) shows the amplitude anomaly due to the first and all overtones. A clear oscillation is noted with epicentral distance, broadly consistent with the behaviour noted in raw catalogues. The precise periodicity of this interference at larger epicentral distances depends on the details of shallow structure. This effect may therefore appear more clearly at shorter epicentral distances in global aggregate comparisons. Based on this analysis, we concluded that measurements of fundamental-mode dispersion are influenced by how overtones are accounted for in the measurement technique. Isolating the fundamental-mode signal with windowing in the time domain may not completely alleviate the problem for distances (30-44 • , 48-68 • , 72-150 • ) and period ranges (60-125 s) where there is substantial contamination from overtone arrivals. We adopted interference filter criteria (Table 2) to exclude data within these ranges from catalogues that measure only fundamental-mode dispersion (GDM52, Scripps14).

Inter-catalogue consistency of overtone measurements deteriorates compared to that of fundamental modes, but remains sufficiently high to permit reconciliation amongst a subset of catalogues. Figs 9 and S14-S25 summarize catalogue scatter comparisons at 50 s from different overtones branches n = 1-6, demonstrating that the mean discrepancies remain low across branches. For the 1st overtone branch, mean differences can exceed 3 s for 50 s Love (e.g. Cambridge19-Utrecht08) and Rayleigh waves (Lyon18-Utrecht08); mean absolute Downloaded from https://academic.oup.com/gji/article/228/3/1808/6408466 by University of California School of Law (Boalt Hall) user on 03 January 2022 differences increase with great-circle distance in a non-monotonic fashion reaching ∼8 s at 120 • . Comparisons for higher overtone branches reveal a similarly high level of agreement (e.g. 50 s Dublin13-Utrecht08) and a clear trend with epicentral distance, although fewer common measurements are available. The median absolute differences in phase anomalies from various catalogues show a less clear trend as a function of period for 1st and 2nd overtone branches (Fig. 10) than is observed for the fundamental modes. Overall, the discrepancies across overtone branches 1-6 are uniformly low (<5 s) for the period bands specified in Table 2. Utrecht08 and Dublin13 catalogues show consistently higher levels of agreement than other catalogues across all mode branches except the first overtone. Discrepancies increase with great-circle distance for all combinations of catalogues with no clear relation to the overtone branch, indicating differences in the reference 1-D Earth models or geodetic constants.

Construction of a reference data set includes calculation of summary rays and removal of outliers for clean catalogues. These analyses provide further insights into the sources and estimates of uncertainties in the reported measurements.

The raw catalogues contain quality-controlled dispersion measurements interpolated on a set of reference frequencies along cross-validated source-station paths. Starting with the raw catalogues, we used a homogenization procedure to construct summary rays and evaluate consistency in measurements traversing similar paths. Previous studies on summary travel times have evaluated path similarity by finding the locations nearest to the source and receiver from a prescribed set of basis regions (e.g. Engdahl et al. 1998). Instead, we used K knot locations that are spaced nearly uniformly on the Earth's surface and assumed no prior knowledge of tectonics or data coverage (Fig. S1). The knot locations are given by the n-fold tesselation of a spherical icosahedron (Wang & Dahlen 1995a) For each wave type, frequency and catalogue of raw dispersion measurements, we evaluated summary rays between pairs of knot points. We chose K = 2562 splines with an average knot spacing 0 i of 4.33 • as the underlying grid for the homogenization process (Fig. S1). When assigning the original rays to pairs of knot locations, k and k , we corrected the observed propagation jth phase anomaly (δ j , where j = 1, 2,..., N in Table 1) by applying a multiplicative factor that accounts for the differences in path lengths, according to

where D is the minor-arc distance between knot pairs, j is the minor-arc distance of the original path, N o the orbit number and N c is the number of times circled around the Earth. The correction factor above accounts for the total distance traversed by the wave and is therefore larger for higher-orbits waves (R3-R5, L3-L5). We excluded subsets of catalogues when none of the knot pairs satisfied the minimum number Downloaded from https://academic.oup.com/gji/article/228/3/1808/6408466 by University of California School of Law (Boalt Hall) user on 03 January 2022

Figure 11. Intra-catalogue deviations in the raw fundamental-mode Love and Rayleigh phase-anomaly measurements. Median deviations (eq. 12, in seconds) along similar paths are calculated for fundamental-mode measurements at all periods for each catalogue (top row). For most catalogues, substantial variation in the median deviations with epicentral distance is not detected (bottom row). Dashed vertical lines denote the Quality filters adopted during the construction of the reference data set (Table 2).

of two contributing original paths, for example Rayleigh wave overtone data from Dublin13 (2nd Over. : 250-350s; 3rd Over. : 250-300s, 4th Over. : 150-200s) and MBS11 catalogues (2nd Over. : 250s, 4th Over. : 75s). Homogenized data and deviations/uncertainties across all knot pairs (k, k ) are defined, respectively, as the mean and standard deviation of all associated corrected measurements ( kk j ).

Deviations in measurements along similar paths, defined as those sharing the same knots for the source and receiver locations, can provide an empirical measure of consistency within each catalogue. We summarized intra-catalogue deviations by the median value of the standard deviations across all available pairs of knots in the homogenized data (eq. 12). This procedure was repeated for each catalogue, wave type, period and within bins of epicentral distance, and reported in seconds. Phase deviations in a self-consistent, high-quality catalogue are expected to be small relative to a cycle (≤π /4 radians). Fig. 11 shows the variation in intra-catalogue deviations of the raw fundamental-mode Love and Rayleigh phase-anomaly measurements as a function of frequency and epicentral distance. Love waves show a median variation in the range of 2-4.5 s and Rayleigh waves in the range of 2.5-6 s for all catalogues at these periods. Some catalogues exhibit larger deviations in measurements; Cambridge19 shows values higher by up to 2 s relative to the other catalogues, suggesting numerous outliers. For all catalogues except Cambridge19 and IPGP03, no substantial variation in the median deviations with epicentral distance was detected between 30 • and 120 • . Overall, intra-catalogue deviations along similar paths are small (≤5 s) and similar in the period range 25-250 s across several catalogues, and there are no systematic trends with epicentral distance. Catalogue deviations along similar paths are substantially higher for the overtone data than for fundamental modes (Fig. 12). The median deviations are the highest for the first (n = 1) overtone branch and decrease with overtone number n. Median deviations of IPGP03 are substantially higher than in the other catalogues (up to 15 s), with median deviations in the first and second overtone branch in excess of ∼10 s at all periods and distances. The large deviations are likely due to the limited geographic coverage of IPGP03 catalogue (Fig. 3), which has fewer overtone measurements than other catalogues (Table 1). Moreover, the starting models for the inversion of phase velocities in IPGP03 are typically the large-scale best solutions from non-linear optimization, which can differ substantially from PREM (Beucler & Montagner 2006). MBS11 also exhibits lower consistency in its overtone dispersion data, with median deviations similar to IPGP03 for the first overtone branch (∼5-10 s). The overlapping measurement periods cover a narrower period band for the third and higher overtone branch (50-60 s), where the median deviations are also uniformly low (±5 s) for all catalogues except the IPGP03 catalogue. Except for Lyon18 and Cambridge19, none of the larger catalogues show a clear trend of median deviation with great-circle distance for overtones in agreement with the fundamental-mode data.

Outlier identification and removal is critical to ensuring consistency across catalogues and robustness of phase-slowness inversions. It is particularly important to assess the relative quality of the measurements and flag inconsistencies because (semi)automated methods can be improved by our filter criteria. While these methods enable fast data processing and generation of large sets of dispersion data, they lack the detailed oversight of a domain expert inherent in fully supervised techniques. We obtained a clean data set on original paths and an associated clean homogenized data set for each catalogue after the removal of outliers. Our definition of an outlier is based on, (1) large intra-catalogue Downloaded from https://academic.oup.com/gji/article/228/3/1808/6408466 by University of California School of Law (Boalt Hall) user on 03 January 2022 Figure 12. Intra-catalogue deviations in the raw overtone phase-anomaly measurements. Median variations (eq. 12, in seconds) along similar paths do not vary substantially with period (left-hand column) or epicentral distance (left-hand column) for most catalogues. Note that this trend corresponds to an increase in phase error (in radians) with frequency in Fig. 4. Dashed vertical lines denote the Quality filters adopted during the construction of the reference data set (Table 2). deviations on similar paths, (2) clear half (±0.9-1.1 • π ) or full-cycle discrepancies (±0.9-1.1 • 2π) on original paths, (3) large inter-catalogue inconsistencies on similar paths (|δφ kk 2 -δφ kk 1 | > 0.4 • 2π), (4) quality criteria of distance and period ranges where the signal is most easily measured across techniques, (5) source depths and magnitudes where excitation is strongest and (4) criteria based on overtone interference (Section 4.5). Table 2 summarizes the configuration of some of these quality control criteria in our processing scheme.

Because measurements associated to a knot pair (k, k ) correspond to very similar source-station paths, we excluded all knot pairs and contributing paths if the standard deviations exceeded 10 times the median of corrected phase anomalies (δ kk j , eq. 12) across all catalogues. We adopted the median as the preferred value for catalogue comparisons since it is less affected by outliers. Half-and full-cycle discrepancies indicating polarity reversals and cycle skips that were identified in Section 4.5 were also excluded. Finally, we excluded the knot pairs and their contributing paths as outliers when there were large inconsistencies (>0.4 • 2π) between pairs of homogenized catalogues. Fig. 13 shows the effect of quality-control criteria (Section 4.4) and the outlier analysis for 100 s minor-arc Rayleigh wave measurements. Only a small fraction of paths (<0.01 per cent) from the contributed catalogues were excluded as outliers during the creation of the clean data set. A substantial fraction (>25 per cent, pie chart) of the removed outliers in IPGP03 are due to intra-catalogue inconsistencies while the bulk of outliers in other catalogues are due to source and quality filter criteria. A slight improvement in consistency within catalogues Downloaded from https://academic.oup.com/gji/article/228/3/1808/6408466 by University of California School of Law (Boalt Hall) user on 03 January 2022 Figure 13. Effect of quality control and outlier analysis on the various catalogues. Median and interquartile range of uncertainties remain roughly similar between the homogenized versions of raw and clean data (top row). Fraction of paths removed from filters corresponding to Quality control (red), CMT (green), large intra-catalogue uncertainties (blue), overtone interference (pink) and inter-catalogue comparisons (cyan), are provided as pie charts (Fig. 2, Table 2). The range of homogenized phase anomalies and their uncertainties decreases from processed to summary data (middle row). However, only limited original paths are removed between the raw and clean catalogues (bottom row), explaining the minor changes in median uncertainties (top row).

is observed after the removal of outliers based on catalogue median uncertainties. Interquartile ranges of homogenized phase anomalies and uncertainties also decrease from raw to clean catalogues. Since a limited number of original paths are removed between the two catalogues, median deviations do not change substantially between the two catalogues. Removal of outliers has a detectable effect on the RMS variations of phase anomalies in fundamental-mode Rayleigh waves, though the values remain similar to within ±0.5 radians (Fig. 4). The RMS variation in the clean catalogue at periods shorter than 35 s increases by up to ∼1 wavelength in case of Dublin13 data, which has the desirable effect Downloaded from https://academic.oup.com/gji/article/228/3/1808/6408466 by University of California School of Law (Boalt Hall) user on 03 January 2022 Figure 14. Half-cycle (π ) discrepancies between pairs of catalogues, attributed to polarity reversals in waveform data. These issues are preferentially observed at epicentral distances ∼77-92 • and on certain dates and stations.

of improving its consistency with GDM52 and Scripps14 catalogues. This is due to the substantial number of short paths with small phase anomalies in the Dublin13 catalogue, which gets removed during our procedure (Table 2).

Analysis of outliers provides interesting insights on the waveforms from which the measurements were derived. For example, half-cycle discrepancies (π ) between GDM52 and Scripps14 catalogues at some periods (Fig. 14, 100 s) are typically associated with specific stations. Geofon stations KSDI-GE, JER-GE and a few others (VSL-MN, AIS-G) show the most number of discrepancies (10-40 paths). Such outliers are likely due to reversed polarities of stations such as KSDI-GE and AIS-G in certain time periods, which may not be reflected in the instrument response history. The time history of potential reversals can be identified from the comparison of catalogues and the retrieved CMT source mechanism. We found some qualitative agreement between our list of possible polarity reversals and those used in the processing of Scripps14 data. Automatic detection of reversed polarities may require accurate propagation phase predictions from a 3-D reference Earth model.

The large amounts of surface wave measurements analysed in this study presented a major computational burden. At the same time, both inter-and intra-catalogue consistency suggests a level of redundancy in the contributed measurements that justifies constructing a summary data set. First, neighbouring modes within the same nth-overtone branch have very similar theoretical sensitivities to radial structure. Very fine sampling of the dispersion curve may not therefore provide independent structural constraints. Second, dispersion is often contributed in terms of propagation phase anomalies at arbitrarily sampled frequencies (N contrib , Table 1), which may not accurately represent the information content in a catalogue. Due to the approximations inherent to the measurement technique, Dublin13 provides finely sampled dispersion curves resulting in a very large number of contributed data (3-8 × N contrib , Table 1). However, the number of unique source-station paths (N path ) in Dublin13 is comparable or even lower than in other recent catalogues (e.g. Scripps14, Cambridge19, MBS11 and Lyon18).

Our processing scheme results in a summary data set that represents the best estimates and uncertainties of phase anomalies between knot pairs (Fig. 2). Phase anomalies in the summary data set (δ ) were calculated from the median estimates of all clean homogenized catalogues. The knot locations (K = 2562) specified in the homogenization of raw and clean catalogues are therefore preserved in the summary catalogues. Our procedure also accounts for overlapping coverage from various catalogues and inter-catalogue consistencies. Systematic inconsistencies due to geodetic constants and reference 1-D models were averaged out. All reference phases ( 0 P ) reported in the summary data were derived using updated geodetic constants (Section 4.3). Based on our experiments, such baseline issues do not influence substantially the lateral phase-velocity variations that are the focus of this study and most potential applications of the reference data set.

We estimated observational uncertainties of the phase anomalies empirically by comparing measurements for similar paths (Section 5.1). The knot pairs in the summary data set were determined to be of quality A, B or C, depending on the number and consistency of independent constraints from various catalogues. If only a single clean catalogue provided constraints, the knot pair was assigned the poorest C grade when a single source-station path was available, or a B grade when multiple paths were available with highly consistent measurements, that Downloaded from https://academic.oup.com/gji/article/228/3/1808/6408466 by University of California School of Law (Boalt Hall) user on 03 January 2022 is phase anomaly deviations did not exceed the median of catalogue uncertainties. When multiple clean catalogues provided constraints, knot pairs were assigned an A grade if they exhibited high levels of both intra and inter-catalogue consistency. The inter-catalogue consistency was considered high if all catalogues had uncertainties within the 2σ level of their median value. Due to the strict criteria adopted here, only a limited number of knot-pair paths in the summary data set were assigned the A grade.

Observed variations in the measured phase anomalies along similar paths reflect errors in the (1) source location, (2) CMT focal mechanism, (3) source excitation, due to incorrect source depth or local structure, (4) interference of the fundamental mode with overtones, (5) instrument-response correction and ( 6) errors due to seismic noise. Fig. 15 and Table 3 show the estimated uncertainty for the most numerous quality-B observations as a function of frequency, indicating a roughly linear increase in the phase uncertainty with increasing frequency. This trend is observed both in the summary catalogue and the contributing clean catalogues, although some catalogues such as MBS11 and Utrecht08 tend to have steeper gradients. A linear increase in phase uncertainty with frequency is also found in the case of overtone catalogues (Table 3). These observations may be related to possible source mislocation of 16-25 km in the direction to or from the station (e.g. Smith & Ekström 1996;Ekström 2011). It is difficult to reliably disentangle the combined effects of the error sources, especially because the empirical uncertainty reported here may underestimate the true uncertainty of the measurements.

Robustness in the features of mantle heterogeneity may be assessed with the phase-velocity variations inverted from the clean and summary data sets. Outstanding questions include identifying redundancies in the raw catalogues and structural complexities that are largely independent of the measurement technique.

We wish to determine the 2-D variations in local phase slowness or velocity as a function of colatitude θ and longitude ϕ. Such phase-velocity maps are optimized to fit a set of observed propagation phase anomalies at a given period, while accounting for azimuthal variations due to anisotropy. In a weakly anisotropic Earth, azimuthal variations in Love and Rayleigh wave velocities can be described by patterns with simple twofold (2ζ ) and fourfold (4ζ ) azimuthal symmetry (Smith & Dahlen 1973). Such variations are defined in terms of the propagation azimuth ζ with respect to the local meridian, following:

where p (ζ ) is the full, azimuthally varying, phase slowness and A, B, C and D are coefficients describing azimuthal variations in phase slowness with respect to the reference isotropic slowness p 0 . We can invert for lateral variations in anisotropic properties,

where laterally varying coefficients A(θ , ϕ), B(θ , ϕ), C(θ , ϕ) and D(θ, ϕ) are expanded in terms of a finite set of basis functions f i (θ, ϕ) on the surface of the sphere, for example

where a i are the model coefficients, and K is the total number of basis functions used in this representation. Since local phase slowness perturbations (δp/p 0 ) in the period range analysed in this paper can be larger than 20 per cent, we avoid the approximation 1/(1 + δc/c 0 ) ≈ 1 -δc/c 0 , which is commonly used to linearize the tomographic problem for small perturbations δc(θ, ϕ) in local phase velocity. When the results below are discussed in terms of velocity variations, (δc/c 0 ), these variations are calculated from the slowness variations by

In order to represent a smoothly varying direction of azimuthal variations in regions close to the poles, we reference the changing azimuth along a ray path to the 'local parallel azimuth' of a parallel line passing through the nearest spline knot location (Ekström 2006a).

As our set of basis functions, we choose spherical splines defined on an equispaced set of knot locations, given by the n-fold tesselation of a spherical icosahedron (Wang & Dahlen 1995a). The ith basis function f i (θ, φ) depends on the distance from the ith knot point as follows:

where 2 0 i is the full range of the ith spline basis function, and 0 i is set equal to the average distance between knot points. Specifically, we adopt a basis set containing K = 1442 splines with an average knot spacing 0 i of 5.77 • (Fig. S1). 3. The solid (and dash-dotted) lines show the phase error for the reference Rayleigh (and Love) wave data set, defined as the median of uncertainties across all catalogues.

Starting from eq. ( 9), we obtain a set of linear equations to predict the N number of observed phase anomalies for each period and wave type [δ obs. j ; j = 1, 2, ..., N]. These observations are related to the local slowness variations described by the set of spherical spline coefficients as Here, the integration is performed over the source-station path {θ j , φ j } corresponding to the jth observation, summation is over all non-zero splines along the path and p i are the coefficients of the isotropic part of the model. Based on prior studies of azimuthal anisotropy variations (Ekström 2011;Ma et al. 2014) and on expectations from mineralogy (e.g. Montagner & Nataf 1986;Montagner & Anderson 1989), we restrict our inversion to the 2ζ terms (A and B) for Rayleigh waves, and the 4ζ terms (C and D) for Love waves. This imposes the same number of free parameters in the Love and Rayleigh inversions and facilitates evaluation of the appropriate level of model complexity. For Rayleigh waves, neglecting 4ζ variations has previously been shown to not substantially influence the retrieved maps of 2ζ variations (Maggi et al. 2006).

The chi-squared misfit χ 2 which is minimized in solving for optimal p i , a i , b i , c i , d i is expressed as:

where j is the index of the datum, σ j is the observational uncertainty, δ pred. j is the predicted phase anomaly for a given set of spline coefficients and w j is a weight applied to the jth observation.

Weighting can be introduced to lower the contribution from highly sampled paths or increase the importance of paths important for coverage. To facilitate the use of data catalogues presented in this study, we calculate the number of observations (N P ) that share the same starting and ending point on 2562 evenly distributed points (Section 5.1). Each observation corresponding to this path is assigned an intra-catalogue weight

which downweights the contribution to the χ 2 from a path sampled repeatedly in a catalogue (e.g. Ekström 2011), such as those from the active Tonga-Kermadec subduction zone (Fig. 3, Section 4.1). Summary data sets on the evenly spaced grid have more uniform coverage by construction; we therefore set all weights to unity in inversions that utilize these data.

Several corrections are often applied before interpreting the propagation phase anomalies in terms of interior structure. δ ref, ellip denotes correction due to different reference models and Earth's hydrostatic ellipticity, δ ζ represents the correction for azimuthal variations in surface wave phase slowness, while δ crust is the correction due to the strong crustal heterogeneity (Table A1; Moulik & Ekström 2016, the Appendix). The corrected phase anomalies can then be attributed to variations in isotropic phase velocity (c) or slowness (p = 1/c) via eq. ( 9). Our forward operator is linear, and the inverse problem can therefore be expressed as d obs = Gm, where G is the sensitivity matrix derived using GCRA. Since our objective here is not to infer Earth structure, but rather to obtain surface wave phase slowness maps, we ignore all structural corrections in the remainder of this study. The RSDF HDF5 container files include fields that store the data corrections (Table A1).

Regularization is introduced to stabilize our inversion by minimizing the sum of χ 2 and additional terms that quantify the amplitude or roughness of the phase slowness variations. We define the isotropic roughness R as the RMS second-derivative gradient of the global isotropic phase-slowness variations

implemented by a discrete Laplacian on the knot points. Since the 1442 quasi-equispaced knot points are constructed via dyadic subdivision of a spherical icosahedron (Fig. S1), the discrete Laplacian is evaluated for the six nearest neighbours of all but the 12 knots that have only five neighbours. We similarly define the anisotropic roughnesses R 2ζ and R 4ζ as

where s nζ and c nζ are the spatially varying sine and cosine coefficients of 2ζ and 4ζ anisotropy. In a full inversion for isotropic, 2ζ and 4ζ -anisotropic phase-slowness variations, we then choose to minimize the quantity χ 2 , where

γ defines the relative preference assigned to fitting the observations and obtaining a smooth model and λ controls the relative smoothing of anisotropic compared to isotropic maps.

To account for the different sizes of catalogues contributed to this study, we scale the relative weights by the trace of the data sensitivity matrix tr(G T G) and explore the trade-off between data misfit and model roughness (i.e. L-curve analysis) across logarithmically equispaced values of γ between 0.001 and 1000. In Fig. 16, we systematically explore the effect of λ on the power spectra of isotropic and anisotropic variations obtained from the summary data set for 100 s fundamental mode Rayleigh waves. Due to uneven data coverage, λ values less than one lead to unstable results. By varying λ between 1 and 100, we find that values below ∼8 yield spectra of anisotropic variations that increase in power at shorter wavelengths (degrees 6), while those in the 5-13 range produce fairly flat spectra. Larger values of λ suppress much of the anisotropic variations above degree 5. The precise choice of λ does not influence the maps of long-wavelength isotropic variations below degree 12. Since values of λ ≥ 8 yield very similar power spectra of isotropic variations up to degree 20, we adopt this value (λ = 8) in all anisotropic inversions.

The complexity in our inverse problem is controlled both by the number of basis functions K (1442 spherical splines, Fig. S1) and the a priori regularization (e.g. Buja et al. 1989;Hastie & Tibshirani 1990). The amount of information extracted from observations or the effective degrees of freedom for model (N res ) can be approximated as

Downloaded from https://academic.oup.com/gji/article/228/3/1808/6408466 by University of California School of Law (Boalt Hall) user on 03 January 2022 where H is the 'hat' matrix used in a variety of data mining applications and tr( • ) denotes its trace (e.g. Cardinali et al. 2004;Ye 2012;Ruppert 2012). In least-squares inverse problems that are the subject of this study, H is equivalent to the resolution matrix R (Menke 1989). We systematically monitor how variance reduction (VR), χ 2 fits and reduced χ 2 red [defined as χ 2 /(N d -N res )] to a specific number of data constraints (N d , based on Table 1) vary with the strength of smoothing (γ , γ 2ζ and γ 4ζ ) and with the introduction of azimuthal anisotropy. We also calculate the Akaike information criterion (Akaike 1974, AIC = χ 2 +2• N res +c 2 ) and the Bayesian information criterion [Schwarz 1978, BIC = χ 2 + ln(N d ) • N res +c 1 ]. The constant terms c 1 and c 2 cancel out when comparisons ( AIC c , BIC) are made between candidate models using the same subsets of data. In our study, BIC penalizes model complexity more strongly than AIC, as it accounts explicitly for the large number of phase anomaly observations, which are considered independent. In isotropic-only inversions, the strength of smoothing γ is selected to minimize BIC, representing an optimal trade-off between model complexity and data fit. In anisotropic inversions, the damping weight is chosen to match the N res of the corresponding isotropic-only inversion. We do not construct phase slowness maps using data sets with fewer than 1442 measurements.

By analysing phase slowness maps, we assess the extent to which long-wavelength variations can be adequately resolved by a smaller homogenized data set that aims to remove the redundancies in the clean catalogues. Fig. 17 compares isotropic phase-velocity variations for R1 fundamental-mode arrivals at 25 s from the GDM52 catalogue and L1 arrivals at 50 s from the Cambridge19 catalogue. These data types Downloaded from https://academic.oup.com/gji/article/228/3/1808/6408466 by University of California School of Law (Boalt Hall) user on 03 January 2022 have good global coverage and are highly sensitive to the strong structural heterogeneity in the crust and uppermost mantle. Therefore, they are well suited for detecting the limits of consistency in derived structure. Phase-velocity maps were constructed from the larger clean catalogue on original paths and the corresponding homogenized data set. In general, we find excellent agreement between the maps, with similar power spectra and degree-wise correlations exceeding ∼0.9 up to spherical-harmonic degree 30. Additionally, agreement does not decrease substantially at short periods where wavelengths are comparable to the knot spacing used in summarizing data. Therefore, homogenized data sets can be used in lieu of clean catalogues for carrying out phase-velocity inversions across the 25-250 s period band. A substantial portion of the phase anomalies on the original paths carry redundant information, especially when making inferences on long-wavelength velocity variations. Based on the agreement between maps constructed with either data set, we conclude that most of the differences in the retrieved maps can be attributed to the improved uniformity of sampling achieved by the homogenized data sets.

The summary catalogues are similar to the individual homogenized catalogues in their knot locations but represent the super set of consistent measurements. Robust summary data sets would ideally capture most of the data variance in the clean homogenized catalogues and provide compatible structural constraints. First, we find excellent agreement between the isotropic phase-velocity variations (R ∼ 0.9) derived from the summary data set and all contributing clean homogenized catalogues, as discussed in detail in Section 6.3. Secondly, we compare the fits to every clean homogenized catalogue provided by phase-velocity maps inverted (1) from the catalogue itself, and (2) from the corresponding summary data set. For most types of surface waves, variance reductions (VR) to the clean homogenized measurements are very similar from the two sets of maps ( VR = ±5 per cent), suggesting an ability of the summary data set to capture most of the data variance in each catalogue. Some notable exceptions are the overtone measurements where there are large differences in how well the summary data captures the data variance of each catalogue. Potential reasons include substantial differences in the geographical coverage among catalogues (e.g. IPGP versus MBS11) or other systematic inconsistencies in the measurements (Section 4.5).

In addition to the direct comparisons of phase anomaly measurements presented in Section 4.5, the degree of consistency across catalogues can be assessed by comparing the related inferences on interior structure. We constructed maps of (an)isotropic phase-velocity variations from the summary data set of each type of surface wave. Fig. 18 show isotropic phase-velocity variations in Love waves at 100 s constructed using phase-anomaly measurements from the clean homogenized catalogues and the summary data set. For Rayleigh waves, isotropic phase-velocity variations from anisotropic inversions are preferred (Fig. 19, right-hand columns, Section 6.4). There is strong agreement between maps from various individual data sets, and the differences carry nearly indistinguishable structural implications. The maps delineate shallow tectonic features as observed in upper mantle tomographic studies. Geological features in the Southern Hemisphere are more prominent when the homogenized data are used in the inversion due to better coverage and uniform weighting of ray paths. For example, the fast cratonic features in Southern Africa are distinctly seen in Rayleigh wave maps constructed with every clean catalogue and the summary data set (Fig. 19, right-hand columns). Both the power of heterogeneity and level of consistency are strongest at the longest wavelengths (degrees <20). Moreover, data fit considerations also support the predominance of large-scale structure for the RMS variations in the catalogues. The RMS misfits to the homogenized data for degrees <20 are within 10 per cent of those corresponding to the full phase velocity map. Such large-scale features of global heterogeneity can adequately explain ∼90 per cent of the signal in the phase delay measurements.

In order to assess the level of consistency across length scales of heterogeneity, we compared the power spectra and degree-by-degree correlation between maps from various catalogues and the summary data set. Reductions in correlations can result from inconsistencies in the measurements comprising each data set or due to differences in coverage. Relatedly, differences in the power spectra can also reflect differences in the amount of smoothing that minimizes BIC for each catalogue. We found that the agreement between phase-velocity variations persists for fundamental mode surface waves across all frequencies. Fig. 20 shows the power spectra and degree-by-degree correlations of phase-velocity variations inferred for 50, 100 and 200 s Love and Rayleigh waves. For Love waves, high degree-by-degree correlations (R > 0.8) and similar power spectra across clean catalogues persist across degrees 1-15, degrading somewhat at the longest periods (200 s, L1). For Rayleigh waves, very high degree-by-degree correlations and similar power spectra persist across degrees 1-25 for most data sets, with the notable exception of IPGP03, whose small size tends to destabilize the retrieval of structure above degree ∼8. Inversions with the summary data set results in features of 25-250 s Love-wave maps that are highly consistent (up to degrees 15-25) with certain catalogues (e.g. GDM52, Scripps14 and MBS11) more than others (e.g. Dublin13 and Utrecht08). In case of Rayleigh waves, summary data sets are highly consistent up to degree 25 for other catalogues as well (e.g. Cambridge19, Lyon18 and Dublin13). This strong agreement on the patterns of heterogeneity (R > 0.8) for both Love and Rayleigh waves informs the target resolution appropriate for a consensus model of upper mantle structure. Agreement on phase-velocity variations tends to deteriorate with overtone number, reflecting the decreasing sizes and reduced consistencies in catalogues (Sections 4.5 and 5.2). As coverage and data availability decreases, smoothing becomes more important, which results in the differences among power spectra beyond a threshold degree. At 50 s, only the first three Love overtone branches have very high correlations (R > 0.8) at degrees up to 8 (Fig. 21). Given the low numbers and limited coverage of available overtone measurements for the higher branches, this degradation in degree-by-degree correlations is not altogether surprising. Rayleigh overtone maps derived from larger data sets show stronger agreements in inferred structure. Phase-velocity maps constructed from four of the data sets-Cambridge19, Lyon18, MBS11 and Utrecht08-exhibit high degree-by-degree correlations across degrees for branches 1-4 at 50 s and 1-2 at 100 s. The highest correlations are noted between data sets that are derived using closely related methodologies (e.g. Cambridge19 and Lyon18, Section 4.2). For a reference data Downloaded from https://academic.oup.com/gji/article/228/3/1808/6408466 by University of California School of Law (Boalt Hall) user on 03 January 2022 set, however, it is more important to ascertain whether different measurement techniques provide comparable results. As with Love overtones, map consistency degrades with overtone branch, though very high correlations persist up to degree 8 even for the 5th overtone. At 200 s, the 5th overtone maps show very little correlation between any of the data sets. Based on these comparisons, we adopted filter criteria that restrict attention to a subset of overtone branches and period bands during the construction of the summary data set (Table 2).

The systematic behaviour of degree-by-degree correlations reflects a number of structural and measurement factors. Phase velocity maps from the summary and the larger contributed data sets correlate highly (R ≥ 0.8) up to degree 20-25 for both fundamental-mode Love and Rayleigh waves at periods between 50 and 100 s. These data sets have good global ray coverage (Fig. 3) and dominant sensitivity to the heterosphere, a region where shear-velocity variations are the strongest (Dziewonski et al. 2010). Long-period waves (T > 150) sample a distinctly different pattern of heterogeneity below the heterosphere and accrue a smaller overall phase dispersion signal. Phase velocity maps of these long-period Love and Rayleigh waves correlate highly only up to degree ∼10. Overtones afford sensitivity to deeper regions of the highly correlated variations in isotropic phase velocity (Section 6.3). Inversions based on the summary catalogue are somewhat agnostic to the measurement technique and other subjective choices used in the construction of the contributed data sets. We evaluate the patterns of azimuthal anisotropy obtained with joint inversions of isotropic phase-velocity variations. Direct comparisons of isotropic and anisotropic inversions, however, are made difficult by the compounding effects of regularization and parametrization. Regularization through smoothing of the phase-velocity variations is necessary to stabilize the inversion. Smoothing may not always be justified based on geological considerations (e.g. strong chemical boundaries), but is often necessary when the total number of model parameter is large (e.g. including anisotropic terms). The imposed smoothing complicates the direct comparison of isotropic-only and joint anisotropic inversions; changes in data fits cannot be attributed directly to an actual signal of anisotropy or to differences in regularization. As the strength of smoothing increases, the effective number of model parameters decreases and the data fits degrade monotonically. We enabled a direct comparison of data fits of the isotropic-only and joint anisotropic inversions through the following procedure. First, we objectively chose an 'optimal' amount of smoothing (i.e. value of γ ) for the isotropic-only inversion that minimized BIC, which represents an optimal trade-off between model complexity and fit to data. Then, we systematically searched for the value of γ to be applied to the joint anisotropic inversion in order to produce a model with approximately the same number of resolved model parameters, N res . Because the number of resolved model parameters does not exactly match, we compared data fits using the reduced χ 2 red statistic, which partially accounts for the changing number of free parameters in the inversion (N d -N res ).

Fig. 23 shows the fast axes and strengths of 2ζ anisotropy in 100 s fundamental-mode Rayleigh waves from clean catalogues and summary data sets. The magnitude and direction of fast axes are in good agreement between various recent catalogues. Qualitatively, the inversion results from this experiment are similar to those of previous studies (e.g. Trampert & Woodhouse 2003;Ekström 2011;Ma et al. 2014). Prominent large-scale features include east-west fast axes in the central Pacific Ocean and the weakest anisotropy across Eurasia. Only the relatively small IPGP03 catalogue has coverage and resolution that is inadequate to robustly constrain the anisotropic patterns. We also investigated potential trade-offs between isotropic and anisotropic variations in phase velocity. Fig. 19 determined from an isotropic-only inversion alongside those from a joint anisotropic inversion of 100 s Rayleigh wave data sets. Isotropic phase-velocity variations are smoother in the joint inversions since they need to be described by fewer parameters. The most notable feature in the isotropic-only maps is the short-wavelength ripple or streak crossing the Pacific from the northeast to the southwest. This ripple is a stable feature of all high-resolution isotropic inversions of Rayleigh waves between 50 and 100 s period from several clean catalogues (Utrecht08, Dublin13 and GDM52) and summary data sets. In joint anisotropic inversions, the isotropic velocity ripple disappears and the phase velocity anomalies display a smooth increase in the Pacific Ocean basic from east to west, consistent with the pattern expected from increasing age and cooling of the lithosphere (Stein & Stein 1992). An outstanding question is whether the additional model complexity of anisotropy variations may be justified based on data. The improvement in data fit resulting from the introduction of anisotropic terms in the inverse problem is shown in Table 4. The introduction of anisotropic terms reduced χ 2 red for 25-250 s fundamental mode Rayleigh waves and most strongly at 100 s period ( χ 2 red > 0.5). The signal of anisotropy is therefore strong for these types of waves and an 2ζ anisotropic model can explain the observed data significantly better than an isotropic model of similar complexity. Addition of 4ζ azimuthal terms to inversions of Love waves and their overtones did not decrease the reduced χ 2 red at any frequency for models with similar numbers of resolved parameters. It is not clear based on current catalogues if azimuthal anisotropy can be reliably constrained from other wave types on the basis of data fit and parsimony arguments. For example, decreases in reduced χ 2 red fits attributable to the inclusion of anisotropy were less than 1 per cent for all Love overtones at 50 and 100 s period. In case of Rayleigh wave overtones, the reduced χ 2 red decreased less than 2 per cent at 50 s period and only 2-3 per cent at 100 s period. Notably, the

The main result of this study is a global reference data set of multimode surface wave dispersion analysed in collaboration with members of the seismological community and archived in a scalable data format. Our procedure for summary data sets effectively reduces redundancy, homogenizes geographic coverage, and averages out measurement, station and source errors (e.g. Pulliam & Stark 1993). The reference data set comprises phase delays and uncertainties of fundamental-mode Love and Rayleigh waves between 25 and 250 s and progressively narrower period bands for overtones up to the 6th branch. After accounting for modelling approximations and salvaging missing metadata, we demonstrate a high level of consistency across most contributed measurements. High quality of the fundamental-mode measurements is evidenced by their uniformly low deviations (< 5 s) along similar paths for all periods and epicentral distances. A few inconsistent outliers can be attributed to cycle skipping, station polarity issues or overtone interference at specific epicentral distances. While deviations are larger for overtones than for fundamental modes, they remain small compared to the wavelength (<π /4 radians) and show no systematic trend with epicentral distance. Despite complications with measuring overtone dispersion, consistency in the types of waves considered here remains high enough to permit reconciliation. Empirical uncertainties of the summary data set are low but increase systematically with frequency for all overtone branches. Future work on multimode dispersion will benefit from precise measurements at higher overtone branches across wider period bands. Future surface wave studies may converge towards greater agreement if certain choices in various processing schemes are made consistent. Reconciliation of large catalogues from diverse measurement techniques suggest potential guidelines. Open data-sharing policies across seismic networks will permit even larger portions of measurements and metadata to be cross-validated (e.g. Cambridge19), thereby extending the geographic coverage of current catalogues and promoting reproducible research. Strong agreement is noted within and across catalogues derived using semi-automated or supervised techniques. Lack of quality control in largely automated methods make issues like cycle skipping difficult to remove (e.g. MBS11), resulting in phase anomalies that can sometimes be abnormally large compared to the wavelengths (> 10π ). Baseline discrepancies between techniques may be easily avoided or reconciled if, (i) detailed information such as source and station metadata are preserved in scalable data formats (see the Appendix), (ii) standard geodetic constants and dispersion corrections from (an)elastic reference Earth models are used while calculating the propagation phase and (iii) both distance-and frequencydependence of overtone interference are considered while applying quality-control criteria. Multimode dispersion across overtone branches may be difficult to disentangle and measure at short epicentral distances where there is not enough separation in arrivals of different modes. Overtone interference in Love waves can extend to teleseismic distances, contaminating fundamental-mode-only measurements (GDM52, Scripps14). Processing choices such as filtering and windowing may not be adequate for avoiding all interference issues (e.g. Ekström et al. 1997); accounting for mode coupling across all branches may be needed for the calculation of synthetic seismograms.

The reference data set presented here represents the current consensus on the available observations of surface wave dispersion between pairs of locations on the Earth. Other types of observations have been reported based on the routine processing of surface wave arrivals on arrays of three-component broad-band seismometers. The amplitudes (e.g. Selby & Woodhouse 2002;Dalton & Ekström 2006), group velocities (e.g. Ritzwoller et al. 2002;Ma & Masters 2014), arrival angles and polarization (e.g. Laske & Masters 1996;Foster et al. 2014a), and ratios between vertical and horizontal components (ZH ratio), are all theoretically sensitive to fine-scale lateral variations (e.g. Larson et al. 1998;Tanimoto & Rivera 2008). Focusing and defocusing of rays due to lateral heterogeneity can influence the amplitudes of surface waves (A F , e.g. Woodhouse & Wong 1986;Park 1987;Romanowicz 1987;Wang & Dahlen 1994), and thereby constrain phase-velocity variations. Amplitude measurements of fundamental-mode Rayleigh waves can also constrain lateral variations of shear wave attenuation in the upper mantle (e.g. Romanowicz 1995;Gung et al. 2003;Dalton et al. 2008;Adenis et al. 2017;Karaoglu & Romanowicz 2018). Such emerging methods of characterizing surface wave arrivals provide sensitivity complementary to that of the propagation phase, especially to constrain small-scale elastic variations that are beyond the scope of the REM3D project. We focused on constructing a reference data set of large and diverse catalogues of propagation phase anomalies currently available from the community.

We adopt the centroid locations based on the Global CMT project and do not solve simultaneously for the earthquake hypocentres in our phase-velocity inversions. In case of inversions that utilize surface waves in isolation, source relocations are typically small (<15 km) but can be substantial in localized regions (e.g. Andes) and influence the patterns of azimuthal anisotropy in Rayleigh waves (Ma & Masters 2015). The average error in centroid locations of the Global CMT catalogue due to lateral heterogeneity and the presence of noise is ∼10 km (Hjorleifsdottir & Ekström 2010), which is unlikely to influence the long-wavelength anisotropic patterns in our study. Such small errors in the centroid locations may be due to the incorporation of longer periods, which are sensitive to different parts of the Earth's structure, and using the full waveform of three wave types (body, surface and mantle) rather than using the first arriving body waves in isolation (e.g. Smith & Ekström 1996). A systematic assessment of the trade-offs between azimuthal anisotropy and centroid locations would require joint source-structure inversions with multiple data types, which is beyond the scope of this study.

Propagation phase measurements from diverse catalogues of multimode surface wave dispersion imply similar large-scale variations in the Earth's mantle. Lateral phase-velocity variations based on summary data sets can explain most of the data variance in the much larger Downloaded from https://academic.oup.com/gji/article/228/3/1808/6408466 by University of California School of Law (Boalt Hall) user on 03 January 2022 (10-1000 times) clean catalogues along original paths. Features of heterogeneity derived from the clean catalogues and summary data set are highly correlated (R > 0.8) in their long-wavelength variations for both fundamental-mode measurements (l max = 15) and overtones (l max = 8). Only the twofold (2ζ ) azimuthal variations in fundamental-mode Rayleigh waves are mapped consistently across catalogues and improve significantly the fits to the reference data set. Inclusion of both major-arc (R2, L2) and higher-orbit arrivals (R3-R5, L3-L5) improves constraints on the even-degree phase-velocity variations. Our inferences on interior structure are based on the simplifying assumption that surface waves travel along the great circle connecting the source and receiver; other theoretical formulations need to be evaluated in future work. As the current consensus data set of multimode surface wave dispersion, this reference seismological data set will provide robust constraints on the 3D reference Earth model (REM3D).

This material is based on work supported by National Science Foundation (NSF) Grant EAR-1345082 and the David and Lucile Packard Foundation. We also thank the Computational Infrastructure for Geodynamics (http://geodynamics.org), which is funded by the NSF Grants EAR-0949446 and EAR-1550901.

. Scatter density plots of raw overtone phase-anomaly measurements (δφ). This is an expanded version of the subplot in Fig. 9. Note the caption in Fig. 5 for reference. Figure S24. Scatter density plots of raw overtone phase-anomaly measurements (δφ). This is an expanded version of the subplot in Fig. 9. Note the caption in Fig. 5 for reference. Figure S25. Scatter density plots of raw overtone phase-anomaly measurements (δφ). This is an expanded version of the subplot in Fig. 9. Note the caption in Fig. 5 for reference.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

Measurements of surface wave dispersion need to be cross-validated against source and station metadata to reference the appropriate seismic waveform and facilitate reproducible research. However, the analysis of large amounts of data from diverse catalogues is computationally inefficient with conventional ASCII and other text-based formats. We have devised reference seismic data formats (RSDFs) in consultation with the community to efficiently process diverse data sets typically used in seismic tomography. While a general definition of a surface wave format is desirable, practical considerations such as diverse processing techniques, computational expertise and utilization of legacy code prevent such an outcome throughout the community. RSDFs format guidelines (Table A1) encourage easy tracking of critical metadata information as header information that will work within existing processing schemes. Headers of RSDF ASCII files store all notes and details relevant to the analysis such as the radial reference Earth model and the associated phase velocity (PVEL) used for calculating the reference phase ( 0 P , eq.7). Other metadata fields correspond to the assumptions used in calculating the various columns in the data such as those contributing to the predicted phase anomaly (δ pred , Moulik & Ekström 2016, Appendix A). The format presented here includes the fields necessary for data reconciliation; other (meta)data specific to a processing scheme can be preserved at the discretion of the analyst.

Table A1. Reference surface wave data format: each ASCII file contains header information followed by columns describing each measurement. All ASCII files are assimilated into a compressed HDF5 container format for distribution and usage in high performance computing (HPC) platforms. The observed uncorrected phase anomaly (delobsphase, δ obs ) w.r.t. the reference phase (refphase, 0 P ) may be compared with the total prediction (delpredphase, δ pred ), which includes contributions from both azimuthal (δ ζ ) and isotropic variations in surface wave phase slowness, crust (delcruphase, δ crust ) and due to Earth's ellipticity of figure (delellphase, δ ellip ), all of which are provided in seconds (Moulik & Ekström 2016 Global reference data sets: surface waves 1849

The archival of millions of surface wave measurements and associated metadata in the RSDF definition requires an efficient container format. Each catalogue is stored in a data container with the metadata as attributes and folders for each wave type, overtone branch and frequency. We chose HDF5 (Hierarchical Data Format version 5;The HDF Group 1997-2015) after extensive testing on various computational platforms. First, HDF5 is widely used as a standard format for data exchange and is operable on various operating systems and platforms. Programs in major languages (e.g. FORTRAN, C and Python) can interface with HDF5 files and there is an active ecosystem with an abundance of libraries and tools. Its usage is steadily increasing in seismology with the related NetCDF 4 definition regularly used during the archival of Earth models and in the processing of waveforms (e.g. Krischer et al. 2016). Second, HDF5 fulfills our requirement of efficient parallel I/O with MPI (message passing interface; MPI Forum 2009) that is needed to process and compare various sets of catalogues. Third, allied features such as the built-in data compression algorithms and data corruption tests in the form of check summing facilitate efficient archival. In contrast to binary formats, HDF5 container formats do not need to account for the endianness of the data and the associated compatibility issues across computational environments. Downloaded from https://academic.oup.com/gji/article/228/3/1808/6408466 by University of California School of Law (Boalt Hall) user on 03 January 2022