Phase velocities of Rayleigh and Love waves in central and northern Europe from automated, broad-band, interstation measurements

Abstract

The increasingly dense coverage of Europe with broad-band seismic stations makes it possible to image its lithospheric structure in great detail, provided that structural information can be extracted effectively from the very large volumes of data. We develop an automated technique for the measurement of interstation phase velocities of (earthquake-excited) fundamental-mode surface waves in very broad period ranges. We then apply the technique to all available broad-band data from permanent and temporary networks across Europe. In a new implementation of the classical two-station method, Rayleigh and Love dispersion curves are determined by cross-correlation of seismograms from a pair of stations. An elaborate filtering and windowing scheme is employed to enhance the target signal and makes possible a significantly broader frequency band of the measurements, compared to previous implementations of the method. The selection of acceptable phase-velocity measurements for each event is performed in the frequency domain, based on a number of fine-tuned quality criteria including a smoothness requirement. Between 5 and 3000 single-event dispersion measurements are averaged per interstation path in order to obtain robust, broad-band dispersion curves with error estimates. In total, around 63,000 Rayleigh- and 27,500 Love-wave dispersion curves between 10 and 350 s have been determined, with standard deviations lower than 2 per cent and standard errors lower than 0.5 per cent. Comparisons of phase-velocity measurements using events at opposite backazimuths and the examination of the variance of the phase-velocity curves are parts of the quality control. With the automated procedure, large data sets can be consistently and repeatedly measured using varying selection parameters. Comparison of average interstation dispersion curves obtained with different degrees of smoothness shows that rough perturbations do not systematically bias the average dispersion measurement. They can, therefore, be treated as random but they do need to be removed in order to reduce random errors of the measurements. Using our large new data set, we construct phase-velocity maps for central and northern Europe. According to checkerboard tests, the lateral resolution in central Europe is ≤150 km. Comparison of regional surface-wave tomography with independent data on sediment thickness in North-German Basin and Polish Trough confirms the high-resolution potential of our phase-velocity measurements. At longer periods, the structure of the lithosphere and asthenosphere around the Trans-European Suture Zone (TESZ) is seen clearly. The region of the Tornquist-Teisseyre-Zone in the southeast is associated with a stronger lateral contrast in lithospheric thickness, across the TESZ compared to the region across the Sorgenfrei-Tornquist-Zone in the northwest.

Time-series analysis, Surface waves and free oscillations, Seismic tomography, Wave propagation; Europe

1 INTRODUCTION

The rapid growth of seismic networks in western and central Europe has produced dense coverage of the entire region with broad-band stations (Fig. 1). The improving coverage makes it possible to map the lithospheric structure of Europe with increasingly high resolution (e.g. Yang et al.2007; Legendre et al.2012; Zhu et al.2012).

Figure 1.

Seismic station coverage in Europe. Triangles indicate permanent and temporary stations that provide data since before 1990 (yellow) and after 1990 (red). Top right: the stations of the German Regional Seismic Network (GRSN) used to illustrate the measurements below.

Phase velocities of the fundamental Rayleigh and Love surface-wave modes are essential observables for the study of the Earth's lithosphere and sublithospheric mantle. In the period range 5–300 they are sensitive to the structure of the crust, mantle lithosphere and asthenosphere (e.g. Aki & Richards 1980; Ben-Menahem & Singh 1981; Dahlen & Tromp 1998; Laske et al.2011). In order to maximize the resolution, we wish to obtain highly accurate measurements in broad period ranges and to utilize as much data as are available. Given the large data volumes, full automation of the measurements is thus highly desirable.

In this paper, we present an automated technique for the measurement of phase velocities of (earthquake-excited) fundamental-mode surface waves in very broad period ranges and an application of the technique to a very large volume of data from across Europe. Phase-velocity maps computed with the new data show detailed lateral structural variations. We test the resolution of the images both with standard resolution tests and by means of comparisons with tectonic boundaries and sediment thickness data.

Phase velocities of Rayleigh waves depend primarily on the velocity of vertically polarized S-waves, but also on density and P-wave velocity. Love-wave phase velocities depend mainly on the velocity of horizontally polarized S waves and on density. The dispersive nature of surface waves (Gutenberg 1924) is caused by the increasing sensitivity of the waves at longer periods to deeper structure within the Earth, with an additional effect from the sphericity of the Earth (Alterman et al.1961). Main characteristics of Rayleigh-wave dispersion curves include a phase-velocity increase from about 2 to 3 km s⁻¹ at short periods below 10 s to about 4 km s⁻¹ at intermediate periods (20–50 s), reflecting the velocity increase across the Moho (Fig. 2, bottom right). Another strong increase, to more than 5 km s⁻¹, occurs at periods longer than 100 s. The strong sensitivity of Rayleigh-wave phase velocities to the Moho depth and their weaker but measurable sensitivity to the lithosphere-asthenosphere boundary depth have been analysed and discussed by, for example, Bartzsch et al. (2011) and Lebedev et al. (2013). In these studies, it was concluded that highly accurate phase-velocity measurements with errors well below 1 per cent are needed in order to resolve the essential details of the S-wave velocity structure of the crust and upper mantle.

Figure 2.

Example of an interstation, Rayleigh-wave, phase-velocity measurement for the station pair BFO (Black Forest Observatory, Germany) and CLL (Collmberg, Germany) (see Fig. 1 for station locations.). The 2-hr long seismograms and the time–frequency representations of their waveforms are shown in the left panels. Top right: the cross-correlation function and its time–frequency representation. Bottom right: the 2π ambiguous phase-velocity curves (blue lines), plotted together with the reference model (black dashed); red line: the manually selected portion of the curve that is accepted as the dispersion measurement.

Love-wave phase velocities are higher than Rayleigh-wave phase velocities for periods lower than about 120 s and show a steadier increase with period (see, e.g. Section 3.2). The latter property is related to the less focussed depth sensitivity of Love compared to Rayleigh waves: for Love waves, the sensitivity at greater depths increases with increasing periods but the sensitivity to a specific depth range is less pronounced compared to Rayleigh waves (e.g. Dahlen & Tromp 1998; Lebedev et al.2013). At very long periods Love-wave phase velocities are normally lower than Rayleigh-wave phase velocities.

Early systematic surface-wave measurements were carried out as single-station measurements of their group velocities by Gutenberg (1924) and Ewing & Press (1950), with group traveltimes from the source to stations estimated as a function of period in the time domain. Time–frequency analysis for single-station group-traveltime measurements was introduced by Landisman et al. (1969). It was improved and used extensively in the following decades (e.g. Levshin et al.1989; Ritzwoller & Levshin 1998; Vdovin et al.1999). Recent examples of single-station measurements are given by Lebedev et al. (2006), based on waveform inversion, and by Foster et al. (2014). Single-station measurements yield phase or group velocities along the entire source-station paths and are affected by errors in the source location and the source mechanism, as pointed out, for example, by Muyzert & Snieder (1996) and Levshin et al. (1999).

In order to minimize errors due to the event location and source mechanism uncertainties, Brilliant & Ewing (1954) measured, in the time domain, Rayleigh-wave phase differences between stations. This eliminated the phase shifts due to the source mechanism and due to the propagation of the fundamental modes from the source to the first of the stations. In order to diminish the influence of the source on phase-velocity measurements, Toksöz & Ben-Menahem (1963) measured phase velocities in the frequency domain using successive passages of surface waves at a single station. McEvilly (1964) measured the phase difference between two stations in the frequency domain for both Love and Rayleigh waves and discovered the so called Love-Rayleigh discrepancy (caused by radial anisotropy, i.e. the difference in the velocities of horizontally and vertically polarized S-waves). With the two-station method, distant sources may be used to study the average local structure between stations in the region under consideration. The influence of the unknown (or uncertain) source parameters is minimized and the problem of the 2π ambiguity of phase-velocity measurements is diminished because of the relatively small interstation distances, compared to epicentral distances in single-station measurements.

Over the last decade, interstation phase-velocity measurements have been widely performed using stations of broad-band arrays and regional networks (e.g. Meier et al.2004; Endrun et al.2004, 2008, 2011; Prindle & Tanimoto 2006; Yao et al.2006; Lebedev et al.2006, 2009; Deschamps et al.2008; Darbyshire & Lebedev 2009; Yoshizawa & Ekström 2010; Beghein et al.2010; Adam & Lebedev 2012; Agius & Lebedev 2013). The interstation dispersion measurements can yield phase-velocity curves in a very broad frequency band. For regional networks (station spacing on the order of 100 km) the period range may span from below 10 to 200–400 s. Lebedev et al. (2006) compared the interstation phase-velocity measurements made using cross-correlation with those computed from single-station measurements. In general, the bandwidth of the dispersion curves is broader for the interstation cross-correlation measurements, in particular at higher frequencies. Interstation phase velocities from single-station measurements are complementary, extending the bandwidth towards lower frequencies. Using local networks with finer station spacing (on the order of 10 km), phase velocities may be obtained for periods as low as 3 s (Endrun et al.2004).

Interstation group or phase velocities may also be determined from cross-correlation of ambient noise or of the coda of surface waves excited by earthquakes (e.g. Campillo & Paul 2003; Shapiro & Campillo 2004; Shapiro et al.2005; Yang et al.2007). An advantage of the ambient noise cross-correlation is the independence from earthquake sources of seismic waves. Inhomogeneous distributions of seismic events or ambient noise sources may introduce a bias, which can be estimated and corrected, to an extent (Paul et al.2005; Roux et al.2005; Yao et al.2006). There is a considerable overlap of the period ranges investigated using the cross-correlation of ambient noise (about 5–50 s, in regional-scale applications) and the cross-correlation of earthquake data (from 5–10 to 400 s, in regional to global applications). Phase velocities calculated from cross-correlations of ambient noise and of earthquake waveforms yield complementary information and may be jointly inverted for 3-D S-wave velocity structure (e.g. Yao et al.2006; Yang et al.2008; Köhler et al.2012).

It has been shown that lateral heterogeneity may lead to considerable perturbations in single-event interstation phase-velocity measurements (Wielandt 1993; Friederich et al.1994; Forsyth & Li 2005; Pedersen et al.2006, 2013; Maupin 2011). In order to overcome this problem the propagation direction in the region under consideration may be measured, as has been suggested already by Press (1956). Several algorithms have been proposed for single-event array analysis, including the f-k analysis based on the assumption of a plane wave (e.g. Capon 1970; Maupin 2011; Sun & Helmberger 2011), the approximation of the incoming wavefield by two interfering plane waves Forsyth & Li (2005) or the determination of local phase velocities from phase as well as amplitude measurements using the Helmholtz equation and allowing for curved wave fronts (Friederich et al.1994). Prindle & Tanimoto (2006) measured the propagation direction from particle motion analysis in order to correct interstation measurements and Foster et al. (2014) extended the single-station approach including a correction of the propagation direction, estimated from densely spaced mini-arrays.

An alternative approach to suppress the errors due to off-path propagation is the rejection of perturbed interstation phase-velocity measurements and the averaging of many single-event measurements. Using sensitivity kernels for two-station measurements, De Vos et al. (2013) showed that averaging of single-event measurements from events on both sides of the station pair decreases the influence of the structure off the interstation path. However, some sensitivity to strong off-path lateral heterogeneity still remains. Because of the strong frequency dependence of the diffraction effects, they normally result in roughness of phase-velocity curves. Consistent rejection of such perturbed (non-smooth) portions of phase-velocity curves and the averaging over many single-event measurements has been shown to yield smooth and reliable interstation phase velocities, provided that (1) suitably elaborate signal-processing techniques are applied (see Section 2.1), (2) noisy and unreliable measurements are consistently rejected and (3) measurements are averaged for many events, from both propagation directions when possible (Meier et al.2004).

The distribution of permanent and temporary broad-band stations in Europe is highly uneven (Fig. 1). Furthermore, the instrumentation, the time of deployment and also the quality of the stations are heterogeneous. Array methods can only be applied to sub-regions and limited time periods. Interstation measurements, in contrast, are well suited to determine phase velocities of fundamental modes for the entire, very large available data volumes (Yang et al.2007). Previous applications of broad-band interstation cross-correlation measurements of phase velocities using earthquake signals have been limited to relatively small regions, due to the substantial burden of the manual data processing. Typical applications in the Baja California region, central Europe or Ireland, for example, used 17, 10, and 24 stations, respectively (Zhang et al.2007; Roux et al.2011; Polat et al.2012). Here, we present an automated processing scheme designed to obtain smooth, path-average fundamental mode dispersion curves using large data sets recorded by permanent and temporary stations. The main elements of the scheme include (1) a solution for the 2π ambiguity, (2) the automated determination of the frequency range of reliable measurements, (3) the identification of smooth segments of the dispersion curve and the rejection of rough perturbations, (4) the detection of quality problems caused by wrong response information or timing problems and (5) the rejection of outliers.

The advantages of automated processing are the consistency in the determination of phase velocities and the option to repeat the processing with different parameter sets, for example, in order to examine the influence of rough phase-velocity perturbations on path-average dispersion curves. In this study, broad-band interstation Rayleigh and Love phase velocities are determined for Europe. Their quality is evaluated statistically and by the analysis of the phase-velocity maps computed using the new measurements. Our new measurements using the heterogeneous European data set offer new insight into the potential of surface wave tomography, in the region and beyond.

2 METHOD

2.1 Cross-correlation measurement of phase velocities

In interstation measurements the phase difference of Rayleigh or Love waves propagating nearly along the great circle path between two considered stations, is measured on the vertical or transverse component, respectively. The Fourier transform of the fundamental mode waveform u₀(t) may be expressed as

\begin{eqnarray} U_0 (\omega ) = | U_0 (\omega ) | \exp (i \phi _0(\omega )) \approx | U_0 (\omega ) | \exp \left(-i \int _0^\Delta k(\omega , s) {\rm d}s\right) \nonumber\\ \end{eqnarray}

(1)

where ϕ₀(ω) is the phase spectrum of the fundamental mode, k(ω, s) is the frequency dependent wavenumber along the path and Δ the epicentral distance. An estimate of the path average wavenumber, k(ω), or the average phase velocity, c(ω), may then be obtained simply from the phase ϕ(ω) of the spectral ratio between the Fourier transforms of the fundamental mode waveforms recorded at the station closer to the epicentre U₀₁(ω) and the more distant station U₀₂(ω) because of

\begin{equation} \frac{U_{01} (\omega )}{U_{02} (\omega )} = \frac{|U_{01} (\omega )|}{|U_{02} (\omega )|} \exp (i \phi (\omega )) \end{equation}

(2)

with

\begin{equation} \phi (\omega ) = \phi _{01} (\omega ) - \phi _{02} (\omega ) \approx k(\omega ) (\Delta _2 - \Delta _1) \end{equation}

(3)

and

\begin{equation} c(\omega ) \approx \frac{\omega ( \Delta _2 - \Delta _1) }{\phi (\omega ) + 2 n \pi } \end{equation}

(4)

where ϕ₀₁(ω) and ϕ₀₂(ω) are the phase spectra of the fundamental modes at stations 1 and 2, respectively, and Δ₁ and Δ₂ are the epicentral distances for station 1 and station 2, respectively. It is important to note here that the use of epicentral distances instead of the interstation distance means that there is no bias in phase velocity if the event is (slightly) off the great-circle between the two stations. This is discussed further in Section 3. Because of the ambiguity in the phase, the correct n has to be determined by comparison with phase velocities for a background model at a reference frequency.

We note that there are two sources of errors in the estimation of the phase velocity according to eq. (4). The first results from deviations of the fundamental mode propagation from a ray theoretical approximation. If the lateral heterogeneity in the vicinity of the stations is not smooth, then the wave front may deviate significantly from a locally plane wave. The second results from errors in the estimation of the phase difference between the fundamental modes. As has already been pointed out by Wielandt et al. (1987) and Prindle & Tanimoto (2006), the spectral ratio method gives poor estimates of the interstation phase difference if it is determined from spectral ratio of the waveforms recorded at the two stations. This is because the spectrum of the waveform is also influenced by the spectra of higher modes U_j(ω), j > 1, and noise N(ω):

\begin{equation} U(\omega ) = U_0 (\omega ) + \sum _j U_j (\omega ) + N(\omega ). \end{equation}

(5)

Therefore, Levshin et al. (1989), Kulesh et al. (2005) and Laske et al. (2011) apply time–frequency analysis, the wavelet transform or a multitaper technique, respectively in order to isolate the contribution from the fundamental mode in the time–frequency domain before calculating the spectral ratio. In the time–frequency domain, u(ω_n, t), higher modes are separated from the fundamental mode. The higher modes may therefore be downweighted in order to obtain a cleaned time–frequency representation u_w(ω_n, t) that contains only the contribution of the fundamental mode and some noise

\begin{equation} u_w(\omega _n, t) \approx u_0 (\omega _n, t) + n(\omega _n, t). \end{equation}

(6)

From the cleaned time–frequency representations the spectral ratio may be determined in order to estimate the phase ϕ_w(ω_n):

\begin{equation} \frac{U_{w1} (\omega _n)}{U_{w2} (\omega _n)} = \frac{|U_{w1} (\omega _n)|}{|U_{w2} (\omega _n)|} \exp (i \phi _w (\omega _n)). \end{equation}

(7)

This is a much better approximation of the phase difference of the fundamental modes than the spectral ratio of the unprocessed waveforms, although it may be perturbed by noise.

Alternatively to eq. (2), the phase difference of the fundamental modes may be determined from the cross-correlation of fundamental mode waveforms ρ₀(t), with its Fourier transform

\begin{equation} {\it F\lbrace \rho _0 (t)\rbrace = |U_{\rm 01} (\omega )|} {|U_{02} (\omega )|} \exp (i \phi (\omega )). \end{equation}

(8)

where ϕ(ω) is the same as in eq. (3).

The cross-correlation function (CCF) calculated from waveforms, ρ(t), is less affected by uncorrelated noise, and the contribution of the large-amplitude fundamental mode to the CCF is enhanced because of the product of the amplitude spectra. Time–frequency analysis may be performed after the cross-correlation in order to obtain ρ(ω_n, t). Again, windowing is applied (Wielandt et al.1987; Laske & Masters 1996; Meier et al.2004) in order to downweight cross-correlations between the fundamental mode and higher modes and between higher modes. This results in ρ_w(ω_n, t). Time windows may be positioned around the maximum of the filtered CCF (Meier et al.2004) in order to isolate the contribution of the fundamental mode. By this procedure the contribution of higher modes and uncorrelated noise is strongly reduced and the filtered and weighted CCF of the seismograms ρ_w(ω_n, t) is approximately equal to the filtered and weighted CCF of the fundamental modes ρ₀(ω_n, t):

\begin{equation} \rho _{w} (\omega _n, t) \approx \rho _{0} (\omega _n, t). \end{equation}

(9)

The phase difference between the fundamental modes may then be determined from the phase of the filtered and weighted CCF, ϕ_CCF(ω_n), in the frequency domain: for each ω_n a Fourier transform is applied to ρ_w(ω_n, t) before measuring the phase at the frequency ω_n. The phase difference between the fundamental modes at the two stations is then determined according to:

\begin{equation} \phi _{\rm CCF}(\omega _n) \approx \phi _{02} (\omega _n) - \phi _{01} (\omega _n). \end{equation}

(10)

In our implementation, we apply Gaussian filters

\begin{equation} F(\omega , \omega _n) = \exp \left(-\alpha _f \left(\frac{\omega }{\omega _n} - 1\right)^2\right) \end{equation}

(11)

with a width of the Gaussian filter that depends linearly on frequency: |$\alpha _f = \gamma _f^2 \omega _n \Delta t$| in order to optimize the time–frequency resolution where the parameter γ_f is chosen usually between 12 and 20 (here 16) and Δt is the sampling interval in the time domain. Windowing in the time domain is easy after the cross-correlation because dispersion is diminished compared to the original seismogram, narrower windows may be applied, and the window is easier to position. We apply Gaussian windows w(t) in the time domain according to:

\begin{equation} w(t) = \exp \left( - \omega _n^2 (t - t_{\text{max}})^2 / (4 \alpha _w) \right) \end{equation}

(12)

where t_max is the time of the maximum amplitude of the CCF after filtering and the width of the Gaussian window is again a linear function of frequency: |$\alpha _w = \gamma _w^2 \omega _n \Delta t$|⁠. Because the dispersion is stronger for longer interstation paths the parameter γ_w varies linearly so that it is 20 for 400 km interstation distance and 50 for 3000 km interstation distance. It has been pointed out by Laske & Masters (1996); Meier et al. (2004) that a careful selection of the frequency range in which the phase-velocity measurement is accepted has to be performed. This becomes obvious by comparing the two examples in Figs 2 and 3. The left panel in Fig. 2 shows the vertical-component time series at the two stations and the corresponding amplitude of the time–frequency representation where the group arrival of the fundamental mode Rayleigh wave clearly emerges as a ridge across the time–frequency map. Cross-correlation of the two waveforms (top right) reduces the dispersion strongly. The phase ϕ(ω_n) is extracted from the complex spectrum of the band-passed and weighted CCF in the frequency domain (bottom right), and we get a bundle of candidate phase-velocity curves (blue lines) due to the inherent 2π ambiguity of the phase measurement. In manual processing, the analyst would have to pick the correct branch, in this example easily identified by the proximity of one branch to a background model, and select a frequency band in which the measurement is acceptable (red segment).

Figure 3.

An example of a noisy interstation measurement of Rayleigh-wave phase velocity for the station pair MOX (Moxa, Germany) and CLL (Collmberg, Germany) (Locations in Fig. 1). Plotting conventions are as in Fig. 2.

The strength of the cross-correlation technique becomes evident when looking at the higher frequencies in this example. Both waveforms are strongly scattered at frequencies above ca. 60 mHz, a typical observation reflecting heterogeneities in crustal structure along the wave-propagation path. However, the CCF shows that these scattered waves are highly correlated between the two stations, implying that scattering must have primarily occurred outside of the inter station path at distant heterogeneity. The scattered wavefields propagate still approximately on the same great circle path as the direct wave front. It is thus possible, despite scattered signals, to determine a consistent phase velocity well above 60 mHz.

In contrast, the vertical-component waveforms in Fig. 3 are less consistent. As the wavefield at both stations is already quite disturbed, with no clear group arrival in the time–frequency map, it is not possible to determine a smooth Rayleigh-wave dispersion curve from the phase of the CCF. Although one branch falls close to our background model, it is rather bumpy and shows rough perturbations around the expected, smooth dispersion curve. This is indicative of an inconsistent phase of the CCF across the frequency range of the measurement. Obviously, this may occur due to strong noise, interference of different modes or when scattered wavefields are not as well correlated as in Fig. 2, leading to strong variations in spectral amplitude across frequencies and larger uncertainties in the phase determination.

As illustrated by these examples, an automated procedure for the selection of phase-velocity curves must be able to identify the correct branch (the correct n) and select smooth segments of the curve that are in an acceptable distance to a background model. The automated procedure allows us then to process and measure large data sets not only consistently but also repeatedly with varying selection constraints, in order to assess the influence of rough perturbations on path-average phase-velocity curves.

2.2 Automated selection of individual phase-velocity curves

Due to the broad depth span of the fundamental mode depth sensitivity kernels (the Fréchet derivatives of phase velocity with respect to shear wave speed as a function of depth) and due to the gradual changes of these kernels with frequency (e.g. Dahlen & Tromp 1998), smooth dispersion curves are expected for any realistic 1-D Earth model. In an automated interstation method, we thus wish to select the smooth parts of an observed dispersion curve while rejecting those parts that exceed a certain level of roughness.

Figure 4.

Illustration of the selection criteria for automated dispersion measurements. (a) Background model criterion. Grey lines show the dispersion curve computed for a background model (solid) and the maximum deviation thresholds (dashed). Thin black lines are 2π ambiguous curves and the multicoloured curve is the selected branch on which the selection criteria are applied. (b) Smoothness (solid, left y-axis) and length (dashed, right y-axis) criterion thresholds are shown in grey. The multicoloured curve is the ‘smoothness curve’ (eq. 14) and dots show the length of the same-coloured curve segments, plotted at the mean frequency of the segment. In the multicoloured curves, blue segments are rejected by the background model criterion, green segments violate the smoothness criterion, the orange segment violates the length criterion. The red segment is accepted.

We introduce a set of criteria to automatically select smooth parts of the phase-velocity dispersion curves which are illustrated in Fig. 4 and formulated as follows:

• Background model criterion. In the first step, we need to select the correct 2π branch on which the measurement should be performed. As can be seen in the examples in Figs 2 and 3, the branches are strongly diverging with decreasing frequency and identification of the desired 2π branch is thus straightforward by comparing the branches at intermediate to lower frequencies to a reference dispersion curve.

   The divergence of 2π branches at lower frequencies decreases with increasing interstation distance because the velocity difference of neighbouring branches is smaller for greater interstation distances. Therefore, a lower reference frequency is needed than for shorter interstation distances. Here we have chosen to vary it linearly between 0.02 Hz at 400 km and 0.0083 Hz at 3000 km interstation distance.

   As phase jumps may occur across the frequency band, we also test neighbouring curves as they may fulfil the criterion at higher frequencies and in case of segments overlapping with respect to frequency the segment close to the reference model is considered. After identification of the desired 2π branch, we apply our first selection criterion, which is defined by the difference of the measured curve from a path-specific background model:

  \begin{equation} \left|\frac{c(\omega _i) - c_o({\omega _i})}{c_o({\omega _i})}\right| {\times 100}\, { per \, cent} < \text{th}_{\Delta C }, \end{equation}

  (13)

  where th_ΔC is the maximum deviation from the background model in percent, c_o(ω_i) is the phase velocity of the background model, and c(ω_i) is the measured phase velocity at every sample ω_i in the frequency domain.

   Segments of the curve that exceed the defined threshold are rejected (blue segment in Fig. 4). Note that we reject also a number of frequency samples before and after a violating segment, to account for the finite resolution of the measurements in the frequency domain.

   Although a global 1-D model would be sufficient for branch identification, we apply path-specific reference models, which we calculate - for each station pair - as path averages through a 3-D model composed of CRUST2.0 (Bassin et al.2000) and PREM (Dziewonski & Anderson 1981). These models take into account the first order structural variations to which surface wave dispersion are strongly sensitive, primarily Moho depth or the presence of sedimentary basins along the interstation path. These path-averaged models proved to be sufficient to define one globally valid parameter set in our routines.

   The choice of the threshold parameter influences many of the individual measurements, and we shall discuss the parameter choice and the sensitivity of the automated measurements to it in detail in Section 3.1. In general, we prefer to keep this parameter as ‘loose’ as possible, so as to account for sufficient variation from our reference models and to avoid a bias of the measurements towards them.

• Smoothness criterion. While criterion (i) is designed to be relatively loose, we introduce a stricter second criterion to enforce a desired degree of smoothness on the accepted measurements. We quantify the roughness of the curve by calculating the first derivative of phase velocity with respect to frequency c^′(ω), and compare it with the corresponding value of the background model, |$c_o^{\prime } ({\omega })$|⁠:

  \begin{equation} S(\omega _i) = \sum _{\omega _j=\omega _i-d(\omega _i) }^{\omega _i+d(\omega _i)} \left| \frac{{c^{\prime }}(\omega _j) - c^{\prime }_o({\omega _j})}{c_o({\omega _j})}\right|<\text{th}_{S }, \end{equation}

  (14)

  where th_S is, again, a constant threshold. A typical value of this threshold is 150 s for both Love and Rayleigh waves.

   We perform a summation of the first derivative deviation from the reference model over a moving window in the frequency domain. The absolute value is taken so that positive as well as negative deviations are treated equally. The frequency range of the summation, 2d(ω_i), is increasing linearly with frequency, to account for the greater first derivative of the phase velocities at lower frequencies. Therefore, a frequency independent threshold can be applied. We note also that this empirical criterion is sensitive to the length (i.e. the number of samples) of the CCF in the time domain, as a longer time series decreases the sampling interval in the frequency domain and the summation in eq. (14) is carried out over an increased number of frequency samples. The length of the cross-correlation is varied linearly with interstation distance, such that it is 1000 s for an interstation distance of 400 km and 2000 s for an interstation distance of 3000 km. That means the criterion becomes more sensitive for larger distances as required because the smoothness of the phase-velocity curve increases with increasing interstation distance. Therefore, the threshold is independent of frequency and interstation distance and is easy to apply. Again, we reject a number of samples before and after a violating segment as can be seen by the green segments in Fig. 4.

• Length criterion. The routine may select any number of smooth segments of the phase-velocity curve with variable length. Very short segments, however, are of little use and are typically determined with less confidence, particularly at higher frequencies where deviations in the dispersion curves from the background model may be greater, due to greater heterogeneities in the crust. In order to avoid too short segments, we apply the length criterion. This is frequency dependent, relaxed at lower frequencies and more strict at higher frequencies (Fig. 4B). The length of the accepted segments has to be greater than the threshold th_Δω:

  \begin{equation} {\rm th}_{\Delta\omega} = \max \left \{a \ast \log(\omega _{m})+b,\; \text{min threshold value} \right \}, \end{equation}

  (15)

  where ω_m is the centre frequency of the segment. Typical values are a = 0.0088 Hz, b = 0.0524 Hz and min threshold value = 1/200 Hz. These values have been determined empirically.

  With increasing interstation distance the ambiguity in the phase velocity due to cycle skipping is increasing. Therefore, an upper frequency limit of the phase-velocity curve is chosen such that the difference between neighbouring phase-velocity curves is larger than a certain threshold, chosen here to be 0.1 km s⁻¹.

  In the example in Fig. 4, two segments have passed the background and smoothness criteria, a short segment starting at 7 mHz (orange) and a long segment starting at 20 mHz (red). While the first segment falls below the threshold of the length criterion, the second one passes the criterion and is accepted as a dispersion curve measurement between 20 and 70 mHz.

2.3 Automatic averaging of interstation phase velocities

After making measurements for all available events for a given station pair, we obtain a bundle of measurements, as in the example shown in Fig. 5 for the station pair BFO-CLL of the German Regional Seismic Network (GRSN). In this example, 20 yr of data, in total 476 events, were processed, resulting in 368 acceptable dispersion curves with our preferred smoothness parameter (see Sec-tion 2.4, Fig. 5C). Obviously, the measurements are highly mutually consistent, although with relatively greater variability at frequencies below 60 mHz. We note that the measurements for the two different propagation directions (black and red curves) are highly consistent: the red curves cover the black curves nearly fully. To obtain a single, robust interstation curve from this bundle of measurements we need to average the single-event curves and assess the measurement quality statistically.

Figure 5.

Automatic averaging of individual dispersion measurements for the station path BFO-CLL with a ‘strict’ (A, B), ‘conservative’ (C, D) and ‘loose’ (E, F) set of measurement parameters. (A, C, E) All selected individual measurements, where black and red curves are measured in opposite propagation directions. (B, D, F) Average dispersion curve and standard deviation, after outlier rejection. (G) Map of station locations (triangles) and processed events (circles) for this station pair. (H) The difference of the average curve from the background model and its standard deviation, for the ‘loose’ (red), ‘conservative’ (blue) and ‘strict’ (green) parameters.

In that process we apply the following conditions and selection steps: (1) outliers are rejected (here, the 10 per cent outermost values are rejected), (2) at each frequency a minimum number of measurements (here, 5) is required, (3) the mean phase velocity and standard deviation are calculated separately for the two directions and if the phase-velocity difference for the two directions exceeds a threshold (here, the maximum standard deviation of the two directions) the measurement is rejected at the corresponding frequency, (4) the standard deviation of all measurements should be smaller than a certain threshold (here, 3 per cent), and (5) finally the length criterion as explained in Section 2.2 is also applied to the resulting average curve.

After the application of this secondary selection procedure, the remaining curves are highly mutually consistent and yield a final path-average phase-velocity measurement with low standard deviation (e.g. Fig. 5C).

2.4 Choice of the smoothing parameter

Naturally, the choice of parameters and thresholds is an essential step in fine tuning the method for application to a large data set. Out of the three criteria, the smoothness criterion turns out to be the most sensitive to the measurement quality. If we relax the smoothness criterion, the bandwidth will increase, the smoothness of the curve will decrease and the standard deviation may also increase. Therefore, there is a trade-off between maximum bandwidth, least standard deviation and smoothness. The first criterion, the limit on the difference with the background model, has only a minor influence compared to the smoothness constraint. The effect of the length criterion is obvious.

We illustrate the effects of varying smoothness parameters by comparing a ‘loose’ (250s), ‘conservative’ (150s) and ‘strict’ (50s) parameter setting. The effect becomes clear in Fig. 4(b) where the green segments of the curve would be accepted or rejected in the ‘loose’ or ‘conservative’ parameter sets. Clearly, the ‘loose’ setting finds rougher curves acceptable. With a ‘strict’ choice, only a narrowband segment between 25 and 50 mHz would be accepted.

Fig. 5 shows a comparison of all measurements for the station pair BFO-CLL located in Germany using these three parameter settings.

Acceptance rates drop significantly in the ‘strict’ setting, while they are not very different for the ‘loose’ and ‘conservative’ choices. The bundle of accepted curves (Figs 5C and e) is broader and rougher in the ‘loose’ setting. After outlier rejection and averaging (Figs 5B, d and f), the curve bundles are smoother, with the standard deviations increasing from the ‘strict’ to the ‘loose’ settings. At the same time, the ‘strict’ setting results in a narrower bandwidth of the average curve, compared to the other two settings, and this appears to be the most significant overall difference in the average curves (Fig. 5h). The narrowing of the frequency range of the measurements is undesirable, particularly at higher frequencies. We thus exclude the ‘strict’ setting from further consideration and from the more detailed parameter analysis of the effect of ‘loose’ and ‘conservative’ setting, discussed in Section 3.1. Although the value of the smoothing parameter may not be very intuitive in itself, it can be chosen appropriately by performing only a few tests.

3 APPLICATION TO CENTRAL AND NORTHERN EUROPE

We applied the method to a large data set, recorded by 1073 seismic stations in and around Europe (Figs 1 and 7). We used the EHB (Engdahl et al.1998) when available and NEIC catalogues and searched for suitable events from January 1990 to October 2013.

For each station pair with an interstation distance between |$1{\deg }$| and |$30{\deg }$|⁠, we conduct a search in the catalogue for all events that have (1) an epicentral distance between |$5{\deg }$| and |$120{\deg }$|⁠, (2) a minimum magnitude of 4–6 (increasing linearly with epicentral distance) and (3) a maximum depth of 100 km. Furthermore, the backazimuth of the selected events must lie within |$7{\deg }$| from the azimuth of the great circle connecting the two stations.

The choice of the maximum backazimuthal deviation for the events is sometimes thought to bias the apparent phase velocity of a plane wave between the two stations by cos⁻¹(ε), where ε is the deviation of the backazimuth to the interstation azimuth (Pedersen et al.2006). This happens, however, only if the interstation distance is used in eq. (4). We avoid this in our implementation by using the epicentral distances instead. A bias may also occur due to wave propagation off the expected great-circle path, caused by lateral heterogeneity between the source and the stations (Weidle et al.2010). Finally, although the source-station geometry is accounted for correctly in eq. (4), the deviations from the interstation great circle should still be reasonably small, because for larger deviations the measured phase difference may be less representative of the structure between the stations, and the measurement error in the phase velocity will increase (see Section 3.2).

A typical choice of a |$5{\deg }$|–|$10{\deg }$| limit is commonly applied and allows for robust dispersion measurements (e.g. Endrun et al. (2008), Adam & Lebedev (2012)). To verify that this assumption does not introduce a systematic bias in dispersion measurements of real, non-planar surface waves, we separate measurements from events with great-circle deviations of up to |$3{\deg }$| from the remaining events with deviations |$3{\deg }< {\varepsilon }\; {\le }\; 7{\deg }$| (Fig. 6).

Figure 6.

(A) Comparison of the automated phase-velocity measurements for subsets of events with different backazimuthal deviation from the interstation great circle (station pair BFO-CLL, as in Fig. 5). Events with backazimuths in the ranges |$0{\deg }-3{\deg }$| and |$3{\deg }-7{\deg }$| are shown in green and red, respectively. Note that the latter has more than twice as many events. (B) All measurements plotted together. (C, D) Measurements for the two opposite propagation directions, for |$0{\deg }-3{\deg }$| and |$3{\deg }-7{\deg }$|⁠, respectively. The total number differs slightly from that in Fig. 5(c) because individual segments of phase-velocity curves are counted here (can be >1 for any event), as compared to the number of events used for the measurements in Fig. 5. (E) Average curves and their standard deviations. (F) Differences of the curves in panel (E) from the background model.

A choice of |${\varepsilon }\; {\le }\; 3{\deg }$| reduces the number of events by, roughly, a factor of three, with the remaining individual dispersion measurements overall smoother than for the events in the range |$3{\deg }<{\varepsilon } \; {\le }\; 7{\deg }$| (Figs 6C and d). However, after the application of the averaging routine, the resulting average phase-velocity curves are barely distinguishable (Fig. 6e), with slightly larger standard deviations when ε is increased. Note also that the sign of the difference in phase velocity between the subsets is not constant and appears to vary randomly (Fig. 6f), showing that there is no systematic bias to faster velocities when increasing ε to |$7{\deg }$|⁠.

We retrieve waveform data through the European Integrated Data Archive (EIDA) infrastructure. An automated routine based on ObsPy (Beyreuther et al.2010) is used to request approximately 4.6 million waveforms through the Arclink interface of WebDC (http://www.webdc.eu) for a total of 364.939 station pairs.

Of all requests, we retrieved around 1.37 million waveforms and performed more than 12 million interstation measurements with acceptable measurements for 164 231 station pairs for Rayleigh waves and 133 736 station pairs for Love waves. As a result, we have obtained 1.6 million acceptable individual (single-event) and 63 000 interstation average phase-velocity curves for Rayleigh waves and ca. 1.3 million individual and 27 500 interstation average curves for Love waves.

The number of processed events per station pair (Fig. 7B) is greater for station pairs with an NE to SE azimuth, owing to the high event rates in East Asia. For paths with a more northerly azimuth, there are significantly fewer events, with moderate event rates in eastern Mediterranean and the northeastern Pacific. Another factor that is seen here is that most of the older seismic stations are located in Central Europe, covering a larger portion of the total time span of the data set. When applying our automated dispersion measurement procedure to the data set, on average, 1/4 of the measurements per interstation path are accepted (Fig. 7C). Typically, this amounts to around 200–400 measurements for station pairs with older stations in continental Europe, and around 20–30 for more recently or temporarily installed stations (Fig. 7d).

Figure 7.

Overview of processed data and retrieved measurements.

3.1 Quality control

The advantages of the automated processing are its consistency and the possibility to assess the effect of different parameter settings on the final average phase-velocity curves of the entire data set. Here, we examine the differences between the results given by two parameter sets, one with a ‘loose’ and another with a tighter, ‘conservative’ setting (smoothness parameters 250 and 150 s, respectively).

This comparison illustrates the influence of rough perturbations in the individual phase-velocity curves on the average phase-velocity curve. The influence of random noise is expected to cancel out during averaging. The influence of higher modes, however, may bias the phase velocity towards higher velocities. Also, deviations of locally plane waves from the expected propagation direction may cause a positive bias of the average phase velocity. As discussed previously, this does not immediately apply to earthquakes that are slightly off the interstation great-circle because in eq. (4) the epicentral distances are used.

Large deviations of surface-wave rays from the predicted great-circle direction—around |$20\deg$|—are observed (e.g. Weidle et al.2010), but these phenomena affect only a specific frequency range, for example around 20–40 s when the wave front is deflected by strong variations in the crustal structure or the Moho depth. Such effects are reflected by distinct ‘bumps’ in the individual phase-velocity measurements that are usually rejected by our selection criteria. By comparing wave fronts propagating in both directions between stations of a given station pair, such strong wave-front deviations are easily detectable, because the waves propagate through different tectonic settings between the sources and the stations.

Furthermore, the approximation of the wave front by a plane wave may be insufficient (e.g. Wielandt 1993; Friederich et al.2000; Cotte et al.2000; Forsyth & Li 2005; Pedersen et al.2006; Bodin & Maupin 2008). A curved and scattered wave front may be convex or concave, leading to frequency dependent perturbation in the phase velocities that may bias the phase velocities towards larger or smaller values. It is possible but very unlikely that a sharp lateral heterogeneity will result in smooth perturbations in the phase-velocity curve. These different perturbations are thus expected to cause rough perturbations in the phase velocity, and the proposed automated processing is well suited to test if they bias the results towards higher velocities.

Therefore, we expect that a ‘loose’ setting might lead to ‘faster’ measurements, as more deflected and scattered waves pass our measurement acceptance criteria. In a statistical sense, when we compare distributions of velocity differences ‘loose’ − ‘conservative’, the distributions should thus tend to positive values, both in the mean and the skewness. Fig. 8 shows this comparison and the statistical parameters are summarized in Tables 1 and 2. To assess potential effects at different frequency ranges, we calculate the distributions in three frequency bands. Most importantly, from Fig. 8 and Tables 1 and 2, it becomes obvious that the average value as well as the most likely value do not change significantly between the two selections.

Figure 8.

Difference (‘loose’ − ‘conservative’) in phase-velocity measurements for the entire data set obtained with a ‘loose’ (250 s) and a ‘conservative’ (150 s) roughness threshold used in data selection.

Table 1.

Statistical parameters of the distribution shown in Fig. 8(a) (Rayleigh waves).

Period	Mean	Std	Var	Skewness	Kurtosis
100–200 s	0.048	0.49	0.24	1.86	53.47
25–100 s	0.044	0.28	0.08	2.24	46.63
10–25 s	0.049	0.52	0.27	−0.08	4.86

Table 2.

Statistical parameters of distribution shown in Fig. 8(b) (Love waves).

Period	Mean	Std	Var	Skewness	Kurtosis
100–200 s	0.008	1.10	1.21	2.72	38.26
25–100 s	0.091	0.63	0.40	5.34	98.98
10–25 s	−0.007	0.58	0.34	−0.06	9.36

For Rayleigh waves the mean difference is slightly but insignificantly positive with a maximum velocity deviation of 0.05 per cent. Only for the intermediate period Love waves the mean is slightly higher (0.091 per cent). Standard deviations are overall on the order of 0.5 per cent, slightly elevated for long-period Love waves where number of data and measurement quality is poorest. The skewness values tend to slightly positive values, implying more positive outliers in the velocity differences. The kurtosis of all distributions is positive and >3. A kurtosis of 3 is expected for a normal distribution. The large values of the kurtosis reflect the overall peak-like character of the distributions and imply rather few larger outliers beyond the standard deviation. The difference in kurtosis between 10–25 s and 25–100 s Love waves—where the distributions in Fig. 8 are rather similar—is explained by larger outliers in the 25–100 s period range. The distribution for 100–200 s Love waves is significantly different as only approximately half the number of data are available as compared to the lower period bands (ca. 3500 versus ca. 7000). The distribution of the difference between the ‘loose’ and ‘conservative’ selection implies thus in conclusion that the variability we measure is due to data noise and rather random perturbations in the phase velocities caused by complicated wavefields.

When we compare a potential bias in the phase-velocity measurements as a function of interstation distance (Fig. 9) we note that the bias is, as seen before, very small particularly at interstation distances larger than approx. 300 km. At short distances, a bias may be present but on average not larger than around 0.5 per cent for both Rayleigh and Love waves. We may conclude that the bias due to rough perturbations is smaller than expected. But it is necessary to reject rough perturbations as they can introduce significant positive or negative errors in the average phase-velocity curves. The consistent selection and averaging of individual phase velocities is an efficient tool to obtain reliable phase-velocity measurements. We note that numerical modelling is needed to quantify a possible bias of the average phase-velocity curve by smooth perturbations in the individual curves that are similar for both propagation directions. Such perturbations are possible but not likely as they require strong but smooth lateral heterogeneity that varies only slightly with depth. In conclusion, we thus prefer the ‘conservative’ setting (smoothness parameter 150 s) which results in an overall improved standard deviation.

Figure 9.

Difference (‘loose’ − ‘conservative’) in average interstation phase velocities between ‘loose’ and ‘conservative’ parameter sets as function of distance.

The data quality may not only be affected by noise and complexities in the wavefield, but also by technical issues such as incorrect response information, timing problems, or wrong polarities. The automated evaluation of the average deviation from the background model and the comparison of average phase velocities for both directions, allows us to identify stations with data quality problems, which may be indicative of technical issues (Weidle et al.2013).

This is illustrated in Fig. 10 where we show the average standard deviation of all measurements per station and the frequency-averaged difference of velocity measurements in two directions. Large values in either of these maps are mere indicators of potential data quality problems, for example, a difference in the velocities measured in two propagation directions might be related to an incorrect station response information (see Weidle et al.2013, for details).

Figure 10.

Average standard deviation (a) and difference between two propagation directions (b) per station, averaged over all frequencies.

3.2 Phase velocities

A 2-D histogram of our entire phase-velocity data set (Fig. 11) shows that most measurements are in the 10–60 mHz frequency range, with the number of measurements at higher frequencies decreasing gradually, probably due to stronger scattering from shallow crustal heterogeneities. We note, however, that our method produces phase-velocity measurements in uniquely broad period ranges. Many thousands of phase-velocity curves were obtained at periods of 10 s and shorter. The low frequency bound for the majority of the measurements is around 6–8 mHz for Rayleigh and 8–10 mHz for Love waves (although, again, thousands of measurements were obtained at frequencies lower than these).

Figure 11.

2-D histograms of all automatically measured phase-velocity curves of the entire data set. Shown are all individual curves after outlier rejection in the averaging step. Note that the unit of the scale annotation is thousands of measurements.

The total amount of Love wave measurements is only about 20 per cent of the number of Rayleigh wave measurements which is mostly due to the generally larger noise amplitudes on the horizontal than on the vertical components. At periods greater than 60–80s, the reduction of Love wave measurements can also be caused by interference of higher order Love modes with the fundamental one (Schaeffer & Lebedev 2015), or contamination of the transverse component by Rayleigh waves in the presence of lateral heterogeneity and anisotropy.

In Fig. 12, we show the 2-D histograms of the standard deviation and standard error of all the measurements in the data set as a function of interstation distance. For the standard deviation they show the expected decrease of the standard deviation with increasing interstation distance. This is readily explained because according to eqn. (4) a small error in the phase measurement ε ≪ 2π will result in a phase-velocity error ε_c that amounts to |$\vert \varepsilon _c \vert = \frac{\omega ( \Delta _1 - \Delta _2) }{\zeta (\omega )} \cdot (\frac{\zeta (\omega )}{\varepsilon }+1)^{-1}$|⁠, where ζ(ω) is the unwrapped (absolute) phase and thus much larger than ε. As the unwrapped phase increases with increasing interstation distance, i.e. more 2π cycles are added to the phase, the second term and hence the phase-velocity error, decreases.

Figure 12.

Standard deviations and standard errors as a function of interstation distance for all measurements in the entire data set.

In general, we find standard deviations of <2 per cent for interstation distances of and below a few hundred kilometres, decreasing to <1 per cent for paths longer than about 1000 km, both for Rayleigh and Love waves. Given the large number of measurements, the standard error (the standard deviation of the mean) is much lower, indicating an overall uncertainty in the estimated mean phase velocities of <0.3 per cent and <0.5 per cent for Rayleigh and Love waves, respectively.

In order to test the consistency of the phase-velocity measurements, phase-velocity maps and synthetic checkerboard tests were calculated for up to 54 022 station pairs for Rayleigh and 23 852 station pairs for Love waves (as, e.g. in Darbyshire & Lebedev (2009)). In Fig. 13, recovered checkerboards are shown for Rayleigh and Love waves at 12, 30 and 60 s period. In the central and southern parts of our study area, the anomalies are well recovered in terms of their amplitude and shape, as compared to the north, southeast, and southwest. Slight lateral smearing is visible towards the east. These checkerboard tests indicate that the lateral resolution reaches 100 km for periods <∼30s Rayleigh waves in regions with the densest path coverage (see Fig. 7) and about 150 km elsewhere.

Figure 13.

Isotropic checkerboard tests with anomalies of approximately 100 km (left) and 150 km (right) size and 100 m s⁻¹ input amplitude. Solid and dashed outlines mark positive and negative input anomalies, respectively.

The resulting tomographic maps (Fig. 14) show, as can be expected, that at short periods (12 s), both Rayleigh and Love waves are strongly sensitive to the thick sedimentary cover, as in the North German Basin and in the Polish Trough. Comparison with the thickness of the sediments from the model EUCRUST07 (Tesauro et al.2008) shows that our regional surface wave tomography is consistent with compilations of the Deep Seismic Sounding results and confirms the high resolution for the 3-D crustal structure. This agreement also holds when we compare our phase-velocity maps with the model of the Central European Basin System by Maystrenko & Scheck-Wenderoth (2013), probably the most detailed sedimentary-cover model of the area to date.

Figure 14.

Selected isotropic phase-velocity maps. Some major tectonic lines are given for orientation (TTZ, Tornquist-Teisseyre Zone; STZ, Sorgenfrei-Tornquist-Zone; AF, Alpine Front). (a,c,e) Rayleigh and (b,d,f) Love wave phase-velocity maps are shown for (a,b) 12s, (c,d) 30s and (e,f) 60s period. The 12s maps in panels (a,b) are overlain by sedimentary thickness from EUCrust07 (Tesauro et al.2008) (black).

At longer periods, the properties of the lithosphere-asthenosphere system in the Trans-European Suture Zone (TESZ) and in Central Europe are clearly imaged. Northeast of the TESZ high velocities indicates the thick cratonic mantle lithosphere of the East European Craton and the Baltic Shield. It is interesting to note the lateral differences in the lithosphere–asthenosphere system between the Sorgenfrei-Tornquist-Zone (STZ) in the northern part of the TESZ and the Tornquist-Teysseire-Zone (TTZ). The region of the TTZ is associated with a stronger lateral decrease in the lithospheric thickness from the East European Platform towards the southwest than the region of the STZ. Furthermore, the shallow asthenosphere in Central Europe in the region of the Cenozoic volcanic centres like Eifel or Vogelsberg and in the Pannonian Basin is clearly visible. It is also remarkable that the Eurasian slab in the western Alps is imaged clearly, thanks to the high station density in this region. The Love-wave phase velocities at 60 s are mainly sensitive to the mantle lithosphere, as indicated by high velocities in central to northern Europe, and to the shallow asthenosphere, as in the Pannonian Basin where Love-wave phase velocities are low. Love waves also indicate a shallow asthenosphere in the Eifel region. Overall, these results confirm the high quality and rich information content of the large, new set of Rayleigh and Love phase-velocity measurements. A more detailed discussion of the tomography and its geodynamic implications would be beyond the scope of this paper and will be presented in a separate contribution.

4 CONCLUSIONS

An automated processing scheme for the determination of broad-band, Rayleigh- and Love-wave phase velocities is suggested, based on the cross-correlation of earthquake data. Robust dispersion curves are computed as averages over numerous curves and curve segments determined from data from different earthquakes and the same station pair. In our algorithm for the automated determination of interstation, fundamental mode phase velocities, segments of the phase-velocity curves are selected based on (i) their difference from a path-dependant background model, (ii) their degree of smoothness and (iii) their bandwidth in the frequency domain. The thresholds have to be fine-tuned to the data set under consideration and the comparison of phase-velocity measurements for both propagation directions is essential in order to ensure that there are no systematic deviations between the two, which would be indicative of errors. A by-product of the procedure is that stations with erroneous response information or timing problems may be detected.

With the automated processing, large data sets can be processed and the influence of rough perturbations on the phase-velocity curves can be quantified by comparing results from runs with different selection parameters. Rough perturbations (rough curve segments) can have a strong influence on the average dispersion curves and should be rejected. Interestingly, our tests suggest that the deviations can be both positive and negative, with the overall bias due to rough perturbations relatively small and detectable only for distances smaller than about 300 km. Averaging of many individual (single event) phase-velocity measurements combined with the rejection of outliers is necessary to obtain reliable fundamental mode dispersion curves.

The automated procedure has been applied to available data from 1073 permanent stations in central and northern Europe, including more than 1.3 million waveforms. The application of the procedure resulted in a set of high-quality interstation fundamental mode dispersion curves in the period range between about 10 and 350 s for both Rayleigh and Love waves. For most of the dispersion curves, the value of standard deviation was below 2 per cent, and the value of standard error below 0.5 per cent. The low values of the standard deviation and standard error show that the proposed criteria are well suited to select acceptable individual phase velocities measurements.

We inverted our measurements for isotropic phase-velocity maps. According to checker board tests, the lateral resolution in Central Europe varies between about 100 and 150 km and is the highest for intermediate period Rayleigh waves. This represents an increase in the lateral resolution compared to previous tomographic studies of the lithosphere-asthenosphere system in Central Europe.

The new automated procedure makes the processing of large data sets feasible and enables consistent determination of interstation phase-velocity curves. The method is applicable both to permanent stations and networks and to temporary networks, where the smaller amount of data may require loosening of the selection parameters and the available events may not provide sufficient observations in both directions along each path. The obtained data set is complementary to dispersion measurements from ambient noise as the considered frequency range overlaps considerably at short periods but is extended towards longer periods.

This work was funded by German Science Foundation grant ME 1320/4-1. SL acknowledges funding from Science Foundation Ireland grant 09/RFP/GEO2550 and 13/CDA/2192. Waveform data were obtained through webdc (http://www.webdc.eu/). We are grateful to network operators that provide data to the EIDA archives. Furthermore, the usage of data of a number of temporary experiments like TOR, SVEKOLAPKO, EGELADOS, and PASSEQ is acknowledged. This work was a part of the PASSEQ 2006 − 2008 project (Wilde-Piórko et al.2008), station equipment for the PASSEQ experiment were provided by GIPP (Geophysical Instrument Pool Potsdam). Generic Mapping Tools (Wessel & Smith 1998) were used to make some of the figures.

PASSEQ working Group: M. Wilde-Piórko, W. H. Geissler, J. Plomerová, M. Grad, V. Babus̆ka, E. Brückl, J. C̆yz̆ienė, W. Czuba, R. England and E. Gaczyński and R. Gazdova and S. Gregersen, A. Guterch, W. Hanka, E. Hegedűs, B. Heuer, P. Jedlic̆ka, J. Lazauskiene, G. Randy Keller, R. Kind, K. Klinge, P. Kolinsky, K. Komminaho and E. Kozlovskaya, F. Krüger, T. Larsen, M. Majdański, J. Málek, G. Motuza, O. Novotný, R. Pietrasiak, Th. Plenefisch, B. Råuz̆ek, S. Sliaupa, P. Środa, M. Świeczak, T. Tiira, P. Voss, P. Wiejacz.