# Traveltime dispersion in an isotropic elastic mantle: Strong lower mantle signal in differential-frequency residuals

Bernhard S. A. Schuberth, Christophe Zaroli, Guust Nolet

► To cite this version:

Bernhard S. A. Schuberth, Christophe Zaroli, Guust Nolet. Traveltime dispersion in an isotropic elastic mantle: Strong lower mantle signal in differential-frequency residuals. Geophysical Journal International, 2015, 203 (3), pp.2099-2118. 10.1093/gji/ggv389. hal-01256878

HAL Id: hal-01256878

https://hal.science/hal-01256878

Submitted on 9 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Traveltime dispersion in an isotropic elastic mantle: strong lower-mantle signal in differential-frequency residuals

Bernhard S. A. Schuberth,<fnref n="1" /> Christophe Zaroli<fnref n="2" /> and Guust Nolet<fnref n="3" />

Accepted 2015 September 14. Received 2015 August 25; in original form 2015 March 10

## SUMMARY

We study wavefield effects of direct *P*- and *S*-waves in elastic and isotropic 3-D seismic structures derived from the temperature field of a high-resolution mantle circulation model. More specifically, we quantify the dispersion of traveltime residuals caused by diffraction in structures with dynamically constrained length scales and magnitudes of the lateral variations in seismic velocities and density. 3-D global wave propagation is simulated using a spectral element method. Intrinsic attenuation (i.e. dissipation of seismic energy) is deliberately neglected, so that any variation of traveltimes with frequency can be attributed to structural effects. Traveltime residuals are measured at 15, 22.5, 34 and 51 s dominant periods by cross-correlation of 3-D and 1-D synthetic waveforms. Additional simulations are performed for a model in which 3-D structure is removed in the upper 800 km to isolate the dispersion signal of the lower mantle. We find that the structural length scales inherent to a vigorously convecting mantle give rise to significant diffraction-induced body-wave traveltime dispersion. For both *P*- and *S*-waves, the difference between long-period and short-period residuals for a given source–receiver pair can reach up to several seconds for the period bands considered here. In general, these ‘differential-frequency’ residuals tend to increase in magnitude with increasing short-period delay. Furthermore, the long-period signal typically is smaller in magnitude than the short-period one; that is, wave-front healing is efficient independent of the sign of the residuals. Unlike the single-frequency residuals, the differential-frequency residuals are surprisingly similar between the ‘lower-mantle’ and the ‘whole-mantle’ model for corresponding source–receiver pairs. The similarity is more pronounced in case of *S*-waves and varies between different combinations of period bands. The traveltime delay acquired in the upper mantle seems to cancel in these differential signals depending on the associated wavelengths and the length scales of structure at shallow depth. Differential-frequency residuals may thus prove useful to precondition tomographic inversions for the lower-mantle structure such as to reduce the influence of the upper mantle for certain length scales. Overall, standard deviations of the diffraction-induced traveltime dispersion between the longest (51 s) and the shortest (15 s) period considered here are 0.6 and 1.0 s for *P*- and *S*-waves, respectively. For comparison, the corresponding standard deviations of the 15 s residuals are 1.0 s and 2.8 s. In the lower-mantle model, standard deviations are 0.3 and 0.6 s, respectively, which gives an average lower-mantle contribution to the total dispersion of 50 per cent for *P*-waves and 60 per cent for *S*-waves.

**Key words:** Mantle processes; Body waves; Computational seismology; Wave scattering and diffraction; Wave propagation; Dynamics of lithosphere and mantle.

## 1 INTRODUCTION

Detailed knowledge of the seismic structure of Earth’s mantle is of fundamental importance for understanding its dynamics and in particular to quantify the buoyancy forces that drive convection. Both magnitude and inherent length scales of seismic heterogeneity bear key information on the underlying lateral variations in temperature and chemical composition. One example of length scales inherent

<footnote n="1">¹Department of Earth and Environmental Sciences, Geophysics Section, Ludwig-Maximilians-Universität München, Theresienstr. 41, D-80333 München, Germany. E-mail: bernhard.schuberth@geophysik.uni-muenchen.de</footnote>

<footnote n="2">²Institut de Physique du Globe de Strasbourg (UMR 7516 CNRS, Université de Strasbourg/EOST), 5 rue René Descartes, F-67084 Strasbourg Cedex, France</footnote>

<footnote n="3">³Université de Nice Sophia-Antipolis, Centre National de la Recherche Scientifique (UMR 6526), Observatoire de la Côte d'Azur, Géoazur, Les Lucioles 1, Sophia Antipolis, 250 Rue Albert Einstein, F-06560 Valbonne, France</footnote>
to the mantle is the thickness of subducted oceanic lithosphere, which is on the order of 100 km. This thickness reflects the upper thermal boundary layer of the convective system and is governed by the temperature gradient across the lithosphere as well as material parameters such as thermal diffusivity, coefficient of thermal expansion, and viscosity. The same applies, in principle, to the lower thermal boundary layer of the mantle just above the core–mantle boundary (CMB). Thus, inherent length scales in the mantle are ultimately linked to the vigour of convection as expressed through the Rayleigh number (e.g. Kennett & Bunge 2008). The spectrum of seismic mantle structures may then further be altered by the potential presence of chemical heterogeneities, such as for example the differences between crustal and lithospheric compositions that enter the mantle during subduction (e.g. Allègre & Turcotte 1986; Agranier et al. 2005). The distribution and length scales of chemical heterogeneities and their influence on the dynamics of the mantle, however, are still a matter of debate (Christensen & Hofmann 1994; Davaille 1999; Kellogg et al. 1999; Tackley 2000; Hansen & Yuen 2000; Montague & Kellogg 2000; Davaille et al. 2002; Stegman et al. 2002; Jellinek & Manga 2004; Nakagawa & Tackley 2004; McNamara & Zhong 2005; Tan & Gurnis 2007; Bull et al. 2009; Schuberth et al. 2009a,b; Simmons et al. 2010; Davies et al. 2012; Mosca et al. 2012; Schuberth et al. 2012; Li & McNamara 2013; Li et al. 2014).

Despite great progress in the field of global seismic tomography in the last decades, it is yet not possible to image the 3-D structure of the mantle at high enough resolution to cover the full range of dynamically inherent length scales (e.g. Ritsema et al. 2011; Trampert et al. 2013; French & Romanowicz 2014). Tomographic studies on the global scale suffer from uneven data coverage and non-uniqueness of the solution. One is always forced to regularize the inversions, which may result in significantly ‘damped’ and low-amplitude images of the true 3-D mantle structure (e.g. Ritsema et al. 2007). Thus, it is a non-trivial task to accurately constrain the amplitude of seismic velocity variations (e.g. de Wit et al. 2012).

Recently, we have introduced an approach for studying the deep mantle that is complementary to tomographic inversions and which allows for a quantitative assessment of geodynamic models directly against seismic observations (Schuberth et al. 2012). It is based solely on forward modelling and thus largely avoids the problems of limited resolution and non-uniqueness. In this approach, fully synthetic seismic waveforms are computed using a spectral element method for 3-D seismic structures derived from mantle circulation models (MCMs). An important advantage of this approach is that such seismic models feature realistic spectra in case that the simulations of mantle flow are performed at earth-like convective vigour (for a detailed statistical analysis of these models, in particular their power spectra, the reader is referred to Schuberth et al. 2009b).

In a series of earlier papers, we have shown that purely thermal models of mantle flow are compatible with estimates of true polar wander, the geoid, and most important, the magnitude of seismic heterogeneity as imaged by global tomography (Schaber et al. 2009; Schuberth et al. 2009a,b; Davies et al. 2012). In Davies et al. (2012), we have also explicitly demonstrated that thermochemical MCMs exhibit very large variations in temperature in the lowermost mantle and thus predict shear wave velocity variations that are on average two to three times too large compared to those observed in tomographic studies. In the comparisons between MCMs and tomography of Schuberth et al. (2009a) and Davies et al. (2012), the effects of limited resolution were taken into account, but both studies investigated shear wave velocity models only. In Schuberth et al. (2012), we therefore tested MCMs directly against seismic

observations using the new forward modelling approach described above. Simulating 3-D wave propagation in the MCMs revealed that—for both P- and S-waves— the standard deviations of long-period traveltime variations in purely thermal models agree well with those of the observed traveltime residuals. Most important, the synthetic data reproduced both their different trends with ray turning depth as well as their magnitude.

The finding of a predominance of thermal variations in driving mantle flow seems at odds with the potential presence of chemical heterogeneities that are expected, for example, from the subduction of oceanic crust together with the depleted part of the lithosphere mentioned earlier. We note, however, that our studies do not rule out the possibility of chemical heterogeneity, in particular on short length scales. The synthetic seismic wavefields generated in Schuberth et al. (2012) were band-pass filtered at 15 s dominant period before measuring body-wave traveltimes by cross-correlation. Owing to the large width of the Fresnel zone at the ray turning point for waves at this period, the measurements are effectively influenced by structures with length scales of a few hundred kilometres and larger. Any structure smaller than this could potentially ‘hide’ due to wavefield effects. In other words, this wave propagation scenario largely falls into the diffraction regime as defined, for example, in Wu & Aki (1988). Structural length scales in the geodynamic model of ≥100 km are similar to the wavelengths of P- and S-waves with 15 s shortest period (also ≥100 km), and propagation distances are in the teleseismic range of roughly 3000–10 000 km (see Nolet 2008, pp. 58 ff., for details on the four main wave propagation regimes).

Wielandt (1987) was among the first to systematically investigate diffraction of seismic waves, and he proposed that its effects may play an important role in the interpretation of teleseismic delay times. He also pointed out that tomographic studies using onset times may be biased towards faster models due to ‘wave-front healing’ effects, which manifest themselves in different manner for fast or slow seismic anomalies. Along those lines, Nolet (1991), who computed Fresnel zones in Earth’s mantle, stated that ‘The lower mantle may very well seem to possess only small velocity variations because the waves lack the resolution to resolve this. Further research is needed to substantiate this. Such research should be aimed at establishing the frequency dependence of delay times’. The goal of the present study is to do exactly that: to systematically investigate and quantify the diffraction-induced dispersion of teleseismic arrival times in a MCM; that is, in seismic structures with dynamically consistent length scales and realistic magnitudes of the velocity variations.

Several studies have addressed the question of the frequency dependence of seismic traveltimes, scattering and wavefield phenomena in heterogeneous media before. So far, however, all analyses have been limited to media with random distributions of seismic heterogeneity (e.g. Müller et al. 1992; Nolet & Moser 1993; Roth et al. 1993; Shapiro et al. 1996; Igel & Gudmundsson 1997; Baig et al. 2003; Baig & Dahlen 2004; Hong et al. 2005) or single anomalies with simple geometry placed within a homogeneous background model (e.g. Hung et al. 2001; Spetzler & Snieder 2001; Malcolm & Trampert 2011). A number of studies have (mainly qualitatively) looked at wavefield effects of a single slab or plume in an otherwise unperturbed mantle (e.g. Vidale 1987; Cormier 1989; Igel & Ita 1997; Tilmann et al. 1998; Igel et al. 2002). Recently, the influence of 3-D structure in specific regions within Earth’s mantle on traveltime residuals has been investigated based on wave propagation simulations in tomographic models (e.g. To & Romanowicz 2009; Cottaar & Romanowicz 2012). Studies of the latter type, however,
are limited in that tomographic models do not feature all relevant
length scales as mentioned before.

For the purpose of improving the interpretation of seismic data
in terms of Earth structure, theoretical work has also been devoted
to capturing the finite-frequency nature of seismic waves in form of
sensitivity (Fréchet) kernels based on the linear first-order Born ap-
proximation (Marquering et al. 1999; Dahlen et al. 2000; Dahlen &
Baig 2002). Since then, the performance of these so-called banana-
doughnut kernels has been validated against full numerical solutions
of the wave equation (e.g. Hung et al. 2000, 2001; Malcolm & Tram-
pert 2011; Mercerat & Nolet 2012; Mercerat et al. 2014). The work
of Nolet & Dahlen (2000), Hung et al. (2001) and Malcolm & Tram-
pert (2011) has also shown that wavefield effects, such as wave-front
healing, depend not only on the sign, size and magnitude of a seis-
mic heterogeneity and the wavelength of the wave, but also on the
exact location between the source and the receiver. For example, the
traveltime residual imposed by a certain anomaly will be different
if it lies on the ray path (small or zero sensitivity) compared to
an off-ray location where sensitivity has a maximum in the first
Fresnel zone. The residual will also change when moving the
anomaly closer to the source or the receiver, as the overlap with the
volume of non-zero sensitivity (practically speaking) will change.

The situation is much more complicated in case of realistic struc-
tures, which do not consist of a single anomaly within a homoge-
neous half-space only. They rather are combinations of anomalies
of varying size, magnitude, sign and shape. The complex nature of
the wavefield in Earth's interior has, for example, been shown by
Zaroli et al. (2010), who observed traveltime dispersion curves for
body waves that both increase as well as decrease as a function of
increasing period for both early as well as late arrivals (i.e. for waves
that have mainly passed through fast or slow anomalies). It is there-
fore important to statistically analyse the frequency-dependence of
traveltime residuals in a realistic scenario. Furthermore, following
the development of finite-frequency theory, cross-correlation trav-
eltime residuals are nowadays increasingly measured in multiple
frequency bands on regional and global scales (e.g. Sigloch & No-
let 2006; Obrebski et al. 2010; Zaroli et al. 2010; Hosseini & Sigloch
2015). A good understanding of such multifrequency data sets is
essential to fully exploit the contained information on deep Earth
structure. With the large computational resources at hand nowadays,
it is straightforward to solve the full wave equation directly for com-
plex global models of mantle structure, and it is currently possible
to cover a large part of the teleseismic frequency spectrum in such
simulations. Following this approach here, we wish to character-
ize body-wave traveltime dispersion in Earth's mantle from a pure
forward modelling perspective and to illuminate some previously
unrecognized aspects of multifrequency cross-correlation data sets.

The key questions that we want to answer are:

(1) Do mantle structures resulting from thermal convection at
earth-like Rayleigh-number lead to diffraction-induced dispersion
of cross-correlation body wave traveltime residuals?

(2) Is the dispersion significant enough to be relevant for global tomographic inversions?

(3) How large is the contribution of the lower mantle to the seismic dispersion?

To this end, we generate synthetic seismograms for our synthetic
Earth model and measure cross-correlation traveltime residuals for
direct *P*- and *S*-waves in four different frequency bands. This is
done not only for the original MCM, but also for a modified version
in which heterogeneity is removed in the upper mantle. The purpose
of this exercise is to isolate and investigate separately the wavefield

effects associated with structure in the lower mantle. The synthetic
wavefields are computed for purely elastic and isotropic structures
(i.e. no dissipation of seismic energy due to intrinsic attenuation),
and thus, any variation of traveltimes with frequency can be taken
as a sign of structural (i.e. diffraction-induced) dispersion. In other
words, we do not face the problem to separate structural and intrinsic
dissipative dispersion as is the case when dealing with real data.
Note that in fact one could turn the problem around and try to provide
some additional constraints on the potential extrinsic contribution
of mantle structures to the observed ‘apparent’ attenuation. This,
however, is not the purpose of the present study.

We start with a brief description of the models and methods used to generate the synthetic traveltime measurements (Section 2). In Section 3.1, we then take a look at the variation of the geographic distribution of the residuals with frequency to first get a qualitative impression of the dispersion related to elastic 3-D mantle structure. Following up on this, we investigate the difference in traveltime residuals between different period bands for a given source receiver pair. In Section 3.2, we present statistical distributions of the long-period–short-period difference as a function of the associated short-period signal and quantify the seismic dispersion in both the original geodynamic model as well as its modified ‘lower-mantle’ counterpart. In Section 4, we then discuss the dispersion in the sense of ‘differential-frequency’ residuals and elaborate on their similarity between the whole-mantle and the lower-mantle model.

## 2 METHODOLOGY

As noted in Section 1, we aim at analysing the frequency-
dependence of synthetic traveltime residuals of direct body wave
phases. To this end, we follow the approach of Schuberth et al.
(2012): Synthetic 3-D wavefields and seismograms are obtained
from simulations of global seismic wave propagation through geo-
dynamic models of mantle heterogeneity. We use model S09-M2-
Q of Schuberth et al. (2009b) in which variations in seismic ve-
locities and density are predicted for a pyrolite composition by
post-processing the temperature field of an isochemical 3-D MCM
with thermodynamic models of mantle mineralogy (e.g. Stixrude &
Lithgow-Bertelloni 2005, 2011; Piazzoni et al. 2007). Mineralog-
ical models provide elastic parameters and density for the range
of mantle pressures and temperatures. The underlying equilibrium
phase assemblages are computed by Gibbs free energy minimiza-
tion in the CFMAS (CaO–FeO–MgO–Al₂O₃–SiO₂) system. As in
Schuberth et al. (2012), we use the model of Stixrude & Lithgow-
Bertelloni (2005, 2007) here, with seismic velocities corrected for
the effect of anelasticity. Mantle general circulation was modelled
with the parallel finite element code TERRA (Bunge & Baumgard-
ner 1995; Bunge & Richards 1996; Bunge et al. 1996, 1997) that
was implemented on a dedicated compute cluster designed for large-
scale geophysical modelling (Oeser et al. 2006). TERRA solves for
the momentum and energy balance at infinite Prandtl number (i.e.
no inertial forces) in a spherical shell. Compressibility effects are
incorporated in form of the anelastic liquid approximation.

The philosophy behind model S09-M2-Q was to generate mantle
structures based on a minimum set of assumptions on the driving
mechanisms of flow while ensuring earth-like convective vigour
through sufficiently high numerical resolution (~25 km grid spac-
ing globally). The assumptions are: (1) a large-scale flow structure
related to the history of plate motions (Ricard et al. 1993; Lithgow-
Bertelloni & Richards 1998), (2) a radial three-layer viscosity profile
that agrees with post-glacial rebound and geoid observations (Hager
1984; Paulson et al. 2007) as well as with the inferred rate of true polar wander (Schaber et al. 2009), and (3) isochemical whole-mantle flow with strong plume flux (Bunge 2005). A high core heat flow of ~12 TW (about 35 per cent of the outflow at the top of the mantle) is obtained by keeping the CMB temperature constant at 4200 K. The three layers of the viscosity profile are identified as lithosphere, upper mantle and lower mantle, separated at 100 and 650 km depth with assigned viscosities of 10²³, 10²¹ and 10²³ Pa s, respectively. The history of plate motions, which is assimilated into the model as time-dependent surface velocity boundary conditions, controls the location of past and recent subduction, and thus the regions of cold and hence seismically fast material at depth. A detailed statistical analysis of the seismic heterogeneity in model S09-M2-Q demonstrated that it features an earth-like spectrum (Schuberth et al. 2009b). In particular, when the effects of limited tomographic resolution are specifically taken into account, its spectral characteristics compare well to those of tomographic models for the length scales covered by the latter (Schuberth et al. 2009a).

For solving the 3-D wave equation in spherical geometry on the full globe we use the software SPECFEM3D_GLOBE (Komatitsch & Tromp 2002a,b). The size of the numerical grid is chosen such as to ensure sufficient numerical accuracy down to a shortest period of 10 s, and we use the 3-D distribution of absolute seismic velocities and density predicted by model S09-M2-Q projected onto this grid (for details on the implementation of the geodynamically derived seismic model in SPECFEM3D_GLOBE, the reader is referred to Schuberth et al. 2012). Apart from 3-D heterogeneity in the mantle, complexity in the wave propagation simulations was reduced to a minimum by switching off all other options available in SPECFEM3D_GLOBE to compute realistic seismograms; that is, we excluded the effects of anisotropy, topography/bathymetry, ellipticity, etc., in addition to excluding the effects of intrinsic attenuation. For the crust and the core, 1-D seismic velocities and densities of the ‘continental’ version of PREM (i.e. without the ocean layer; Dziewonski & Anderson 1981) are used.

In order to measure traveltime residuals by cross-correlation, reference waveforms are needed. These are also computed with SPECFEM3D_GLOBE using the average 1-D structure of S09-M2-Q for the mantle. The 1-D model is implemented on the exact same numerical grid as the 3-D version of the model to assure that any measured signal is the result of mantle heterogeneity and not

related to numerical artefacts. In particular, using the same mesh ensures that both the 1-D and 3-D signals are affected in the same way by numerical dispersion, which thus should not affect the cross-correlation measurements. For both, the 1-D and 3-D simulations, we directly use the numerical Green’s function in the frequency range resolved accurately by the spectral element simulation without additional convolution with a source time function (i.e. we assume the source is a Heaviside function in time and seismograms are low-pass filtered at 0.1 Hz before any further processing).

For part of the current study, we re-use the set of synthetic seismograms from Schuberth et al. (2012). In addition, we performed a suite of simulations for a modified version of S09-M2-Q in which 3-D heterogeneity is present only below 800 km depth. This second model is called S09-M2-Q-LM and will also be referred to as ‘lower-mantle model’ in the following. The second set of simulations allows us to specifically identify wavefield effects arising in the lower mantle without the strong overprint from structures in the lithosphere and upper mantle.

For both 3-D models, a total of 17 earthquakes was simulated (see Table 1; strictly speaking, it is not correct to use the term ‘earthquakes’ here, as our synthetic model is not the Earth, but for ease of reading, we will use it interchangeably with ‘seismic events’). The events are taken from the Global CMT catalogue (www.globalcmt.org, last accessed 6 October 2015; Ekström et al. 2012) and were selected such as to provide an almost even sampling of our mantle model with seismic waves. The seismic wavefield of each earthquake was ‘recorded’ with a network of 42250 equidistant seismic stations (~100 km apart, on land and at sea), in order to obtain a dense and homogeneous data coverage even with a rather small number of events.

## 2.1 Automated measurement of traveltime residuals

With the event-receiver configuration described above, we obtain around 3 million seismograms for each model. For these 6 million traces, we want to measure traveltime residuals by cross-correlation (see e.g. Nolet 2008, for a definition of such residuals). In synthetic studies, as is the case here, the ‘observed waveform’ is given by the wavefield simulation using the 3-D model and the reference waveform by the one with the 1-D average structure. Measurements

<table><caption>**Table 1.** Earthquakes simulated in this study (taken from the Global CMT catalogue at www.globalcmt.org).</caption><thead><tr><th>No.</th><th>Event name</th><th>Region</th><th>Date</th><th>Latitude (°)</th><th>Longitude (°)</th><th>Depth (km)</th><th>Moment magnitude (M<sub>w</sub>)</th></tr></thead><tbody><tr><td>1</td><td>052382A</td><td>Gilbert Islands Region</td><td>23/05/1982</td><td>-3.35</td><td>177.40</td><td>11</td><td>5.7</td></tr><tr><td>2</td><td>120789A</td><td>Southern Iran</td><td>07/12/1989</td><td>25.94</td><td>59.00</td><td>15</td><td>5.9</td></tr><tr><td>3</td><td>042190C</td><td>Near East Coast of Eastern USSR</td><td>21/04/1990</td><td>47.46</td><td>138.96</td><td>520</td><td>5.5</td></tr><tr><td>4</td><td>041399B</td><td>Fiji Islands Region</td><td>13/04/1999</td><td>-21.42</td><td>-176.46</td><td>173</td><td>6.8</td></tr><tr><td>5</td><td>0062100A</td><td>Iceland</td><td>21/06/2000</td><td>63.98</td><td>-20.76</td><td>15</td><td>6.4</td></tr><tr><td>6</td><td>082100A</td><td>South Atlantic Ocean</td><td>21/08/2000</td><td>-53.02</td><td>-45.97</td><td>15</td><td>6.1</td></tr><tr><td>7</td><td>062003D</td><td>Western Brazil</td><td>20/06/2003</td><td>-7.61</td><td>-71.72</td><td>556</td><td>7.0</td></tr><tr><td>8</td><td>11604D</td><td>Central Mid-Atlantic Ridge</td><td>16/01/2004</td><td>7.64</td><td>-37.70</td><td>15</td><td>6.2</td></tr><tr><td>9</td><td>2070804C</td><td>Southern East Pacific Rise</td><td>08/07/2004</td><td>-25.06</td><td>-115.96</td><td>12</td><td>5.9</td></tr><tr><td>10</td><td>030605A</td><td>North of Severnaya Zemly</td><td>06/03/2005</td><td>84.94</td><td>99.14</td><td>12</td><td>6.3</td></tr><tr><td>11</td><td>200506070534A</td><td>Pacific-Antarctic Ridge</td><td>07/06/2005</td><td>-62.43</td><td>-161.50</td><td>12</td><td>5.8</td></tr><tr><td>12</td><td>200508150753A</td><td>North of Ascension Island</td><td>15/08/2005</td><td>-1.68</td><td>-13.05</td><td>12</td><td>5.3</td></tr><tr><td>13</td><td>200610110600A</td><td>Northern East Pacific Rise</td><td>11/10/2006</td><td>8.40</td><td>-103.17</td><td>12</td><td>5.7</td></tr><tr><td>14</td><td>200704070520A</td><td>South-west Indian Ridge</td><td>07/04/2007</td><td>-39.79</td><td>46.18</td><td>19</td><td>5.7</td></tr><tr><td>15</td><td>200811281342A</td><td>Off Coast of Northern California</td><td>28/11/2008</td><td>40.35</td><td>-126.98</td><td>23</td><td>5.9</td></tr><tr><td>16</td><td>200909071612A</td><td>South of Java, Indonesia</td><td>07/09/2009</td><td>-10.20</td><td>110.63</td><td>36</td><td>6.2</td></tr><tr><td>17</td><td>200909122006A</td><td>Near Coast of Venezuela</td><td>12/09/2009</td><td>10.70</td><td>-67.92</td><td>12</td><td>6.4</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Example of multifrequency cross-correlation traveltime residuals in the geodynamic model S09-M2-Q for a shallow event at the central mid-Atlantic ridge (event 8 in Table 1). Measurements were performed in four different frequency bands using overlapping Gaussian filters with dominant periods of 15, 22.5, 34 and 51 s (grey lines). Top panel: *P*-wave residuals; bottom panel: *SH*-wave residuals. In both cases, the magnitude of the negative delay time decreases with increasing period. Dispersion of the seismic waves results solely from the heterogeneous mantle structure, as the effect of intrinsic dissipative attenuation was not included in the wave propagation simulations.</figcaption></figure>

were done on the vertical and transverse component for the direct *P*- and *S*-waves, respectively. A minimum epicentral distance range of 30° was chosen to ensure a clear separation of the direct phases from later arrivals and to avoid problems related to upper-mantle triplications. As we concentrate on direct *P*- and *S*-waves and wish to minimize interference with *PcP* and *ScS*, the maximum distance is taken to be 90°.

Automated processing of the synthetic seismic data is done using a newly generated *Python* implementation of the algorithm of Zaroli *et al.* (2010). This software performs an automated window selection and subsequent measurement of frequency-dependent traveltimes based on cross-correlation of the trace of interest with the reference seismogram. The underlying algorithm is an extension of *FLEXWIN* (Maggi *et al.* 2009), tailored more specifically to the selection of windows for a set of specific target phases, rather than separating seismograms into portions of coherent seismic energy and portions of noise. After window selection, cross-correlation traveltime residuals were measured in four frequency bands using overlapping Gaussian filters with 15, 22.5, 34 and 51 s dominant

<table><caption>**Table 2.** Standard deviations of the *P*- and *S*-wave traveltime residuals in the four frequency bands investigated in this study for (column 3) the whole-mantle and (column 4) the lower-mantle model. The SMAD (defined in Section 2.1) is employed as measure of scale, which is a robust (i.e. unbiased) estimator for the standard deviation (Kleiner & Graedel 1980). Column 5 gives the ratio between the corresponding SMAD values of the two models.</caption><thead><tr><th>Phase</th><th>Period band<br>(s)</th><th>SMAD(δt)<sub>WM</sub><br>(s)</th><th>SMAD(δt)<sub>LM</sub><br>(s)</th><th>SMAD(δt)<sub>LM</sub><br>SMAD(δt)<sub>WM</sub></th></tr></thead><tbody><tr><td><i>P</i></td><td>15.0</td><td>1.04</td><td>0.45</td><td>0.43</td></tr><tr><td><i>P</i></td><td>22.5</td><td>1.04</td><td>0.45</td><td>0.43</td></tr><tr><td><i>P</i></td><td>34.0</td><td>1.19</td><td>0.45</td><td>0.38</td></tr><tr><td><i>P</i></td><td>51.0</td><td>1.33</td><td>0.45</td><td>0.33</td></tr><tr><td><i>S</i></td><td>15.0</td><td>2.82</td><td>1.33</td><td>0.47</td></tr><tr><td><i>S</i></td><td>22.5</td><td>2.82</td><td>1.33</td><td>0.47</td></tr><tr><td><i>S</i></td><td>34.0</td><td>2.97</td><td>1.33</td><td>0.45</td></tr><tr><td><i>S</i></td><td>51.0</td><td>2.97</td><td>1.33</td><td>0.45</td></tr></tbody></table>

period. Note that the length of the time window used for the cross-correlation measurements increases with increasing period. Further details on the measurement procedure can be found in Zaroli *et al.* (2010). In the paper, we will refer to the measurements at 15 s as ‘short-period’ and to those at 51 s as ‘long-period’ residuals. Fig. 1 shows an example of a multifrequency measurement for model S09-M2-Q. The short-period delay time at 15 s period is negative for both *P*- and *S*-waves, and the magnitude of this negative residual decreases with increasing period (possibly related to more severe wave-front healing at the longer periods).

Only those measurements were retained for which the cross-correlation coefficient (CC) is greater than 0.9. Given the fact that we are dealing with synthetic noise-free waveforms here, the number of measurements removed by the criterion is rather small. Only about 1 percent of the *P*-wave and 6 to 8 percent (depending on the frequency band) of the *S*-wave measurements fall below this threshold. For each model, we obtain a total of around 3 million traveltime measurements (1.5 million for *P*- and *S*-waves each). Table 2 lists the standard deviations of traveltime residuals in the four frequency bands. The scaled median absolute deviation (SMAD) is employed as measure of scale, which is a robust (i.e. unbiased) estimator for the standard deviation (Kleiner & Graedel 1980), and both terms will be used interchangeably in the following. For a set of univariate data $X = (x_1, x_2, \dots, x_n)$, the SMAD is defined as:

$$ \text{SMAD}(X) = s \cdot \text{MAD}(X), \qquad (1) $$

with scaling factor $s = 1.4826$ and the median absolute deviation

$$ \text{MAD}(X) = \text{median}_i(|x_i - \text{median}_j(x_j)|), \qquad (2) $$

which, in other words, is the median (index *i*) of absolute deviations of the data from their median (index *j*). The scaling factor *s* actually is a normalizing constant $1/0.6745$ that strictly is valid only in case of normal distributions and $n \to \infty$.

In the Supporting Information (Appendix SOM2), we provide a possible explanation for the observed increase in standard deviation of the traveltime residuals with increasing period.

## 3 RESULTS

Before we will present statistical analyses of the multifrequency cross-correlation traveltime residuals, we will take a closer look at their geographic distribution for one event and investigate the variation of patterns as a function of frequency.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Maps of traveltime residuals in model S09-M2-Q for an earthquake in Iceland (event 5 in Table 1) plotted at the location of their respective receiver. Left part: P-wave residuals in the four period bands investigated in this study (indicated above each map). Right part: same for S-waves. As we concentrate on direct P- and S-waves, epicentral distances are restricted to the range from 30° (outer edge of white annulus around each CMT solution) to 90° (rim of each map). Note the different colour-scales for P- and S-waves and the dominance of positive delay times, which is a consequence of the event being located in a region of predominantly slow seismic velocities.</figcaption></figure>

**3.1 Geographic patterns of traveltime residuals and their variation with frequency**

Fig. 2 shows maps of P- and S-wave traveltime residuals for each of the four frequency bands. The residuals are plotted at their corresponding station locations for the event in Iceland (event 5 in Table 1). P-wave residuals are about a factor of 2 smaller than the S-wave residuals as expected from the ratio of compressional to shear wave velocity variations in our model. For each period band, all measurements with CC ≥ 0.9 are plotted independent of whether the residuals at the other frequencies are also above this threshold. Given that only a small amount of data is removed by the CC criterion, this choice does not bear on the following discussion (i.e. only very minor effects show up in the maps due to, for example, interpolation).

The pattern of traveltime residuals clearly shows a strong imprint of the near surface structure in model S09-M2-Q in large parts of the maps. This results from the combination of strong heterogeneity in the lithosphere and the increasing concentration of sensitivity of the waves in a smaller volume close to the surface (Schuberth *et al.* 2012). Owing to the fact that the location of structures in our geodynamic model (in particular in the upper mantle) is controlled by a plate motion history model, we are able to use real-world geographical features in the following discussion of traveltime anomalies.

The largest negative residuals are clearly related to regions of subduction (e.g. all around the Pacific Ocean at the Aleutian Arc, Japan and East Asia as well as under North America, where remnants of the ancient Farallon plate are visible as three different slab segments; further strong fast arrivals are visible from the Mediterranean to the Himalayas, which are related to the subduction of the ancient Tethys Ocean). Overall, however, the maps are dominated

by positive residuals as a consequence of the predominantly slow structures around the hypocentre. The dominance of near source structure on the mean residuals for one event is a general and expected feature (Supporting Information Fig. S1 shows a second example with mainly negative residuals for an event located in a region of fast seismic velocities).

Comparing the maps for different frequency bands in Fig. 2, the strongest differences can be seen for P-wave residuals between 15 and 51 s period, as expected. A streak of large positive and negative residuals appears at 34 and 51 s period that runs in radial direction from the event through eastern Siberia and into the Sea of Okhotsk. This linear feature coincides with the direction of a nodal plane in the earthquake source mechanism. As discussed in Schuberth *et al.* (2012), we attribute such features in the patterns of traveltime residuals to diffraction of seismic energy around heterogeneities in the model. Large residuals occur close to the nodal planes when energy leaving the source in one quadrant is diffracted around 3-D structure ending up in the neighbouring quadrant. This results in an arrival of a wave with reversed polarity compared to a wave travelling in the corresponding 1-D structure, which gives rise to large delay times when measured by cross-correlation. Such effects of polarity reversal have also been observed in studies investigating the effects of upper-mantle slab structures (e.g. Silver & Chan 1986; Igel *et al.* 2002).

Close to the nodal planes, the amplitudes of the direct waves are rather small. For real data, measurements along such azimuths would thus likely be discarded in most cases due to a low signal-to-noise ratio. This is not the case here, as noise is absent in our synthetic data. Interestingly, strong signals from polarity reversals show up in our data even with a minimum CC of 0.9. We checked
that the measurements excluded by this threshold in case of the 15 and 22 s period bands do not fall within the geographic region considered here (i.e. along the azimuths of the nodal planes).

A similar, but much weaker feature exists in the maps of S-wave traveltime residuals at 34 and 51 s period. There, a faint linear band of slightly positive values can be recognized that runs radially away from the event through Venezuela. This again is interpreted as diffraction of energy around 3-D structure together with the different radiation pattern for S-wave energy. Further tests would be necessary to verify our interpretation, which, however, is not the scope of this study.

More important, evidence for wave-front healing can also be observed in Fig. 2. The negative signals related to the Farallon slab under North America, for example, experience severe reduction in magnitude with increasing period in case of P-waves. The same is true for the Tethys slab. As expected from their smaller wavelengths at the same frequencies, wavefield effects in the various period bands are smaller for S- than for P-waves (e.g. less reduction of the slab-related signals with increasing period). Nevertheless, a change in shape of the residual pattern can be recognized, for example, for the negative residuals in the Pacific Ocean (which change in amplitude only to a minor extent).

In contrast to the slab-related signals, the positive residuals in Fig. 2 increase in magnitude with increasing period in case of the P-waves. The trend in these signals seems to indicate that the near-source structure (slow seismic velocities surrounding the event in Iceland) has a stronger influence on long-period than on short-period waves. Fig. 3(a) shows the frequency-dependence of the mean P-wave residual for all events. In most cases, the magnitude of the event-mean increases with increasing period, and most pronounced so for those with a large short-period mean. The relation between near-source structure and mean residual is further highlighted in Fig. 3(b). It shows that the mean residual at 51 s period behaves like the P-wave velocity perturbation (inverted in sign) close to each respective event. It is clear that the mean residual is also influenced to some extent by the near-receiver structure and the fact that, depending on the location of the event, more fast or slow structures may be sampled. This may explain why also for seismic events in regions with little heterogeneity around the source, the mean residual is sometimes non-zero (e.g. event 1 in Fig. 3). However, the change of the mean with frequency generally seems to be quite small for these events. Event-means for both P- and S-waves are given Table 3. The S-wave mean residuals show a similar, but weaker and slightly less systematic frequency-dependence compared to the P-waves.

To test whether the frequency dependence of the event-mean is indeed governed by shallow structure, we additionally plot the mean residuals for model S09-M2-Q-LM in Fig. 3(a). In this lower-mantle model, mean residuals are much smaller than in case of the whole-mantle model. The fact that they are not exactly zero indicates a small influence of preferential sampling of fast or slow regions in the lower mantle for certain events. More important, they do not show a strong dependence on frequency (not only in absolute terms, but also relative to the magnitude of the short-period mean), which confirms our above interpretation.

### 3.1.1 Lower-mantle contribution to traveltime residuals

As discussed above, a large part of the measured traveltime residual at each frequency is acquired close to the source and the receiver. The lower-mantle signal is overprinted and partly masked by the effects of structures in the upper mantle and lithosphere and thus somewhat hidden in the measurements. Fig. 4 shows that in case

<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a) Mean P-wave traveltime residuals for each earthquake individually at periods of 15, 22, 34 and 51 s (blue, red, green and magenta colours, respectively). Solid lines/symbols and dashed lines denote the mean residuals in the whole-mantle and the lower-mantle model, respectively. Events are labelled according to Table 1 and sorted along the x-axis according to the mean P-wave residual at 15 s. Note the general trend of increasing magnitude of mean residuals with increasing period. (b) Comparison of (blue line) the mean residual at 51 s period in the whole-mantle model with (green line) the P-wave velocity perturbation (inverted in sign) averaged in a small region with 25 km radius around each respective event.</figcaption></figure>

of heterogeneity being absent in the upper mantle, the traveltime residuals tend to be smaller on average, as expected. The occurrence of positive and negative residuals is more balanced than in case of the whole-mantle model (i.e. the mean residual is close to zero as discussed above; see also Supporting Information Figs S1 and S2 for an event in a region of fast seismic velocities). Interestingly, the geographic pattern of residuals is now strongly related to structure close to the ray turning point. This is illustrated in Fig. 7 and discussed in more detail in Section 4.1. There, we also show that it seems possible (at least in case of our geodynamic model) to reduce the influence of shallow structure on the residuals in the whole-mantle model by combining measurements in different frequency bands and thus to recover the contained information on lower-mantle structure.

Most important, and also in contrast to the whole-mantle model, not only P-wave but also S-wave traveltime residuals show discernible diffraction-induced dispersion in the lower-mantle model. The magnitude of the residuals decreases with increasing period, while the geographic patterns change in shape and tend to appear less sharp at long periods than at short periods (e.g. the large region
<table><caption>**Table 3.** Mean traveltime residuals of *P*- and *S*-waves at different periods (specified in the top row) for each earthquake individually. Events are labelled according to Table 1 and sorted according to the mean *P*-wave residual at 15 s.</caption><thead><tr><th rowspan="2">No.</th><th rowspan="2">Region</th><th colspan="4">⟨δt<sub>P</sub>(T)⟩ (s)</th><th colspan="4">⟨δt<sub>S</sub>(T)⟩ (s)</th></tr><tr><th>15 s</th><th>22.5 s</th><th>34 s</th><th>51 s</th><th>15 s</th><th>22.5 s</th><th>34 s</th><th>51 s</th></tr></thead><tbody><tr><td>2</td><td>Southern Iran</td><td>-1.42</td><td>-1.59</td><td>-2.10</td><td>-2.75</td><td>-4.10</td><td>-4.12</td><td>-4.27</td><td>-4.62</td></tr><tr><td>16</td><td>South of Java</td><td>-0.55</td><td>-0.53</td><td>-0.67</td><td>-0.98</td><td>-1.27</td><td>-1.34</td><td>-1.36</td><td>-1.52</td></tr><tr><td>17</td><td>Near Coast of Venezuela</td><td>-0.23</td><td>-0.34</td><td>-0.53</td><td>-0.76</td><td>-0.35</td><td>-0.42</td><td>-0.51</td><td>-0.62</td></tr><tr><td>7</td><td>Western Brazil</td><td>-0.06</td><td>0.13</td><td>0.25</td><td>0.12</td><td>0.20</td><td>0.21</td><td>0.16</td><td>0.17</td></tr><tr><td>10</td><td>North of Severnaya Zemly</td><td>0.17</td><td>0.18</td><td>0.31</td><td>0.50</td><td>0.24</td><td>0.28</td><td>0.23</td><td>0.18</td></tr><tr><td>3</td><td>Near East Coast of Eastern USSR</td><td>0.25</td><td>0.28</td><td>0.21</td><td>-0.01</td><td>1.10</td><td>1.08</td><td>1.01</td><td>0.97</td></tr><tr><td>1</td><td>Gilbert Islands Region</td><td>0.26</td><td>0.28</td><td>0.27</td><td>0.30</td><td>0.73</td><td>0.74</td><td>0.76</td><td>0.80</td></tr><tr><td>4</td><td>Fiji Islands Region</td><td>0.33</td><td>0.09</td><td>-0.18</td><td>-0.36</td><td>0.25</td><td>0.25</td><td>0.13</td><td>0.01</td></tr><tr><td>14</td><td>Southwest Indian Ridge</td><td>0.45</td><td>0.47</td><td>0.54</td><td>0.70</td><td>0.84</td><td>0.78</td><td>0.85</td><td>0.98</td></tr><tr><td>5</td><td>Iceland</td><td>0.48</td><td>0.61</td><td>1.01</td><td>1.21</td><td>1.03</td><td>1.18</td><td>1.36</td><td>1.37</td></tr><tr><td>6</td><td>South Atlantic Ocean</td><td>0.61</td><td>0.59</td><td>0.55</td><td>0.63</td><td>1.51</td><td>1.52</td><td>1.54</td><td>1.56</td></tr><tr><td>12</td><td>North of Ascension Island</td><td>0.62</td><td>0.76</td><td>0.93</td><td>0.99</td><td>1.08</td><td>1.13</td><td>1.29</td><td>1.38</td></tr><tr><td>15</td><td>Off Coast of Northern California</td><td>0.64</td><td>0.75</td><td>1.01</td><td>1.13</td><td>1.66</td><td>1.65</td><td>1.70</td><td>1.93</td></tr><tr><td>8</td><td>Central Mid-Atlantic Ridge</td><td>0.80</td><td>0.87</td><td>1.05</td><td>1.10</td><td>1.86</td><td>1.80</td><td>1.86</td><td>2.01</td></tr><tr><td>11</td><td>Pacific-Antarctic Ridge</td><td>1.01</td><td>0.99</td><td>1.18</td><td>1.32</td><td>2.69</td><td>2.73</td><td>2.78</td><td>3.04</td></tr><tr><td>9</td><td>Southern East Pacific Rise</td><td>1.04</td><td>1.50</td><td>1.94</td><td>1.94</td><td>2.63</td><td>2.93</td><td>3.19</td><td>3.28</td></tr><tr><td>13</td><td>Northern East Pacific Rise</td><td>1.09</td><td>1.39</td><td>1.57</td><td>1.63</td><td>2.18</td><td>2.41</td><td>2.74</td><td>2.88</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Same as Fig. 2 for the lower-mantle model S09-M2-Q-LM (3-D heterogeneity below 800 km depth only). Note the generally smaller magnitudes of traveltime residuals and the reduced dominance of positive delays compared to the whole-mantle model.</figcaption></figure>

of positive residuals in the Pacific Ocean). Further wavefield effects can be seen from changes in the short-scale patterns, which tend to get smoother and fuzzier as the period of the waves increases.

## 3.2 Quantifying diffraction-induced seismic dispersion in an isotropic elastic mantle

Now, we move to the quantification of diffraction-induced dispersion in our MCM. The following statistical analysis is based on the full data set (i.e. encompassing all 17 earthquakes).

We note again here that only those measurements were retained for which the CC was greater than 0.9.

Fig. 5 shows distributions of $\delta t_{\text{lp}} - \delta t_{\text{sp}}$ (i.e. the difference between the long-period and short-period traveltime residuals measured for a given source-receiver pair, denoted LP-SP difference in the following) plotted as a function of the short-period residual $\delta t_{\text{sp}}$. By plotting the dispersion in this way, the dependence of wavefield effects on the sign and also magnitude of the heterogeneity can be highlighted and possible non-linearities in the relation between the differential signal and the short-period residual may become
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Diffraction-induced traveltime dispersion in model S09-M2-Q. The model features perfectly elastic and isotropic material parameters and structural length scales that are controlled by the dynamics of mantle flow. Statistical representations of the dispersion are shown for (a,c) *P*-waves and (b,d) *S*-waves in (a,b) the whole-mantle and (c,d) the lower-mantle model. Plotted are distributions of the difference between traveltime residuals at the longest and shortest period analysed in this study (51 and 15 s, respectively) as a function of the associated short-period residual (red: median, black: mean, cyan box: 25–75 percentiles, whiskers: 5 and 95 percentiles). The long-period–short-period differences increase systematically in magnitude towards large positive and negative short-period residuals, which indicates wave-front healing. In the lower-mantle model, dispersion relative to the short-period residual seems to be stronger than in the whole-mantle model.</figcaption></figure>

apparent. Positive and negative short-period residuals indicate that the wave has travelled predominantly through slow and fast seismic anomalies, respectively, especially in case of large magnitude measurements. In this sense, Fig. 5 is similar in nature to fig. 17 of Hung *et al.* (2001). In case of simple wave-front healing, one would expect that residuals at long period will generally be smaller in magnitude than at short period. Thus, the LP-SP difference should be positive for fast arrivals (negative short-period residuals) and negative for slow arrivals (positive short-period residuals). Furthermore, the differences should be larger for large magnitude residuals than for small ones, which in the simplest case would lead to a negative linear relation between the LP-SP differences and their associated short-period residual.

One question that we want to answer here is whether the synthetic traveltimes from our geodynamic model meet these expectations

of a linear negative trend or whether they show strong departures from this simple scenario. To check whether wavefield effects are different between *P*- and *S*-waves, the dispersion graphs in Fig. 5 are all plotted on the same *x*- and *y*-scales to ease direct comparison. In fact, *P*-waves are expected to experience stronger wave-front healing than *S*-waves at the same frequency based on their larger Fresnel zone.

Looking at the *S*-wave dispersion in the whole-mantle model first (Fig. 5b), one can see a clear trend of the distributions from positive differences at large negative short-period residuals to negative differences at large positive residuals. This means that the synthetic *S*-wave traveltimes roughly show the expected behaviour. However, the trend of the synthetic data is not exactly linear. The median of the LP-SP differences varies non-monotonically and stays close to zero for short-period residuals between −5 and +5 s. It reaches values
of ~1 s in magnitude for short-period residuals between -10 and -5 s and between +5 and +10 s. For short-period residuals beyond ±10 s, median values reach up to around ±1.5 s. Dispersion is found to be largely symmetric between negative and positive residuals. A slight asymmetry, as is observed here, would be in accordance with the 'Wielandt effect' mentioned earlier; that is, wavefield effects and their impact on traveltime residuals are different for positive and negative velocity perturbations of same shape and magnitude (Wielandt 1987; Hung et al. 2001; Malcolm & Trampert 2011; Mercerat & Nolet 2013). Note that Hung et al. (2001), for example, find stronger wave-front healing for fast anomalies (in particular for strong ones) in case of large propagation distances behind the object. However, the small asymmetry may in our case also be the consequence of the skewed model statistics and different wavefield effects for positive and negative anomalies due to their different dominant length scales and characteristic shapes in the geodynamic model (i.e. elongated slabs versus narrow columnar plumes).

Turning to the results for *P*-waves in the whole-mantle model (Fig. 5a), a similar trend in the dispersion can be recognized (positive LP-SP differences at large negative short-period residuals and vice versa). For short-period residuals between -3 and +3 s, the median of the LP-SP differences is close to zero, also similar to the *S*-wave behaviour. However, deviations from a strictly linear trend are much larger for *P*-waves, as is evident from the negative differences at short-period residuals of around -5 s. Most important to note is that not only the non-linearity but also the dispersion is stronger for *P*- than for *S*-waves. We recall here that in our geodynamically derived seismic model, relative variations in $v_p$ are about a factor of two smaller than their corresponding $v_s$ counterparts. Still, the dispersion of *P*-wave residuals shows large median values of up to ±2 s at short-period residuals of -9 and +8 s. In case of *S*-waves, median values reach only up to about ±1.5 s; and those values only occur at larger short-period residuals of -13 and +14 s there. If the magnitude of seismic anomalies were the same for *P*- and *S*-waves, the difference in dispersion would be even larger. Our result are thus in full agreement with the expectation of stronger wavefield effects for *P*- than for *S*-waves at the same frequency.

In the lower-mantle model (Figs 5c and d), dispersion relative to the short-period residual, as plotted here, seems to be stronger than in the whole-mantle model. Both *P*- and *S*-wave dispersion also show a more linear trend and a stronger asymmetry (with stronger dispersion for negative than for positive short-period residuals). In case of *P*-waves, the median of the LP-SP differences is around

2.5 s at short-period residuals of -5 s and around -1 s at short-period residuals of +4 s. For *S*-waves, the median reaches up to 5.5 s at residuals of -13 s and lies somewhere between -2 and -2.5 s for positive residuals between +10 s and +13 s. Note that the total number of measurements in each distribution diminishes towards larger short-period residuals. Therefore, the distributions beyond ±10 s represent the less frequent extreme cases, which do not contribute much to the overall variance of the dispersion for the entire data set discussed further below.

In Table 4, mean and standard deviation (again given by the SMAD introduced in Section 2.1) of the LP-SP differences computed over the entire data set are listed for all combinations of period bands. Note that the mean values of *P*- and *S*-wave traveltime dispersion are quite small (ranging from 0.024 s to 0.113 s in the whole-mantle model) so that the standard deviation is equivalent to the quadratic mean.

We shortly note here that our findings of small yet non-zero mean values of dispersion might bear on studies of attenuation in the mantle (recall that our simulations are based on perfectly elastic material properties). It is generally accepted that dispersion due to anelasticity (i.e. intrinsic dissipative attenuation) is the main reason for the discrepancy between seismic 1-D models based on long-period normal mode observations and high-frequency body wave arrival times (e.g. Gilbert & Dziewonski 1975; Dziewonski & Anderson 1981). However, the problem of disentangling intrinsic and extrinsic (i.e. scattering and focusing/de-focusing) contributions to the apparent attenuation is a long-standing issue that has not yet been resolved satisfactorily (e.g. Romanowicz & Mitchell 2007; Ricard et al. 2014). A variety of studies have focused on scattering effects in random media, and they consistently found an apparent 'velocity shift' in the traveltimes (e.g. Müller et al. 1992; Nolet & Moser 1993; Roth et al. 1993; Shapiro et al. 1996). In other words, the waves arrive on average earlier than expected from the average propagation velocity in the structure, as they try to avoid slow regions and travel preferentially in the fast regions (in accordance with Fermat's principle). It was proposed that part of the aforementioned discrepancy might be related to such elastic wavefield effects (Nolet & Moser 1993) or to preferential sampling of fast regions owing to the distribution of seismicity and receiver locations (Davies & Bunge 2001). The small mean values of our LP-SP differences could be taken as evidence for little influence of diffraction-induced dispersion on the apparent attenuation. However, though intrinsic attenuation will generally lower the velocities in the less than two octaves (15-51 s),

<table><caption>**Table 4.** Mean and standard deviation of *P*- and *S*-wave traveltime dispersion for all combinations of frequency bands investigated in this study. $\Delta\delta t = \delta t_p1 - \delta t_p2$ is the differential-frequency (DF) traveltime residual with $p1$ and $p2$ the respective period bands, and $\langle \cdot \rangle$ denotes the arithmetic mean. Columns 5 and 6 show the standard deviation (SMAD) for the whole-mantle and the lower-mantle model, respectively. Column 7 gives the ratio between the corresponding SMAD values of the two models.</caption><thead><tr><th>Phase</th><th>DF band<br>(s)</th><th>⟨Δδt⟩<sub>WM</sub><br>(s)</th><th>⟨Δδt⟩<sub>LM</sub><br>(s)</th><th>SMAD(Δδt)<sub>WM</sub><br>(s)</th><th>SMAD(Δδt)<sub>LM</sub><br>(s)</th><th>SMAD(Δδt)<sub>LM</sub><br>SMAD(Δδt)<sub>WM</sub></th></tr></thead><tbody><tr><td><i>P</i></td><td>15.0-22.5</td><td>-0.057</td><td>-0.011</td><td>0.30</td><td>0.15</td><td>0.50</td></tr><tr><td><i>P</i></td><td>15.0-34.0</td><td>-0.113</td><td>-0.010</td><td>0.44</td><td>0.15</td><td>0.33</td></tr><tr><td><i>P</i></td><td>22.5-34.0</td><td>-0.056</td><td>0.001</td><td>0.30</td><td>0.15</td><td>0.50</td></tr><tr><td><i>P</i></td><td>15.0-51.0</td><td>-0.077</td><td>-0.006</td><td>0.59</td><td>0.30</td><td>0.50</td></tr><tr><td><i>P</i></td><td>22.5-51.0</td><td>-0.019</td><td>0.004</td><td>0.59</td><td>0.15</td><td>0.25</td></tr><tr><td><i>P</i></td><td>34.0-51.0</td><td>0.037</td><td>0.004</td><td>0.44</td><td>0.15</td><td>0.33</td></tr><tr><td><i>S</i></td><td>15.0-22.5</td><td>-0.036</td><td>0.001</td><td>0.44</td><td>0.30</td><td>0.67</td></tr><tr><td><i>S</i></td><td>15.0-34.0</td><td>-0.073</td><td>-0.003</td><td>0.74</td><td>0.44</td><td>0.60</td></tr><tr><td><i>S</i></td><td>22.5-34.0</td><td>-0.037</td><td>-0.003</td><td>0.44</td><td>0.30</td><td>0.67</td></tr><tr><td><i>S</i></td><td>15.0-51.0</td><td>-0.098</td><td>-0.038</td><td>1.04</td><td>0.59</td><td>0.57</td></tr><tr><td><i>S</i></td><td>22.5-51.0</td><td>-0.062</td><td>-0.039</td><td>1.04</td><td>0.59</td><td>0.57</td></tr><tr><td><i>S</i></td><td>34.0-51.0</td><td>-0.024</td><td>-0.036</td><td>0.74</td><td>0.44</td><td>0.60</td></tr></tbody></table>
its effect is still small compared to the average diffraction-induced dispersion observed in our model of ~0.1 s: for a wave that spends 40 s crossing the asthenosphere twice, a lowering of Q from 300 to 150 causes a differential attenuation dispersion of only 0.02 s. However, further research is needed to evaluate the influence of the choice of earthquakes used in the simulations, which we postpone to a future publication.

To finally fully quantify the diffraction-induced dispersion in our geodynamic model, we turn to the standard deviation of LP-SP differences in Table 4. For *P*-wave residuals in the whole-mantle model, values range between 0.30 and 0.60 s, while they are slightly larger in case of the *S*-wave residuals for which they range from 0.45 to 1.04 s. In the lower-mantle model, the standard deviations range between 0.15 and 0.30 s for *P*-wave residuals, and between 0.30 s and 0.60 s in case of *S*-wave residuals. Dispersion in an average sense thus is smaller in the lower-mantle model compared to the whole-mantle model. Most interesting, however, *S*-wave dispersion in the lower mantle already accounts for roughly two-thirds of the entire *S*-wave dispersion observed in the whole-mantle model at the frequencies considered here. In case of the *P*-wave residuals, the signal from the lower mantle amounts to around one-half of the whole-mantle signal.

# 4 DISCUSSION

The main result of this study is that mantle structures with length scales that are controlled by the dynamics of flow lead to significant diffraction-induced traveltime dispersion for *P*- and *S*-waves. Interestingly, the dispersion signal relative to the associated short-period residual seems to be stronger in the lower-mantle model than in the whole-mantle model. However, a closer look at the *S*-wave dispersion in Fig. 5(b) reveals that the variation with short-period residual (i.e. the slope of medians) in the whole-mantle model for short-period residuals between -15 and -5 s is similar to that for residuals in the lower-mantle model between -10 and 0 s. The same is true for positive residuals between +5 and +15 s (WM model) and between 0 s and +10 s (LM model), respectively. Between -5 and +5 s, all medians are very close to zero in the whole-mantle model. The same general trend is also discernible from the plots of *P*-wave dispersion, but arguably even more complexity exists there.

The non-linear dependence of the dispersion on the short-period residual observed in case of the whole-mantle model mainly results from two factors. The first, but probably minor, contribution is related to the statistical distribution of seismic velocity variations in the underlying model with which the wavefield interacts. In model S09-M2-Q, the variations in seismic velocities show considerable skewness that changes with depth (positive skew in the upper part of the lower-mantle and negative skew close to the CMB; for details see Schuberth et al. 2009b). As we measured traveltime residuals for a large range of source-receiver distances (i.e. ray-turning depths), it is difficult to tell upfront the exact effect of model skewness on their global statistics. It would be interesting to explore how dispersion in our model depends on the ray-turning depth, but here we opt to stick with the global statistics for now.

A second factor that likely contributes much more to the complexity of the distributions in Figs 5(a) and (b) is that in the whole-mantle model, residuals are offset from zero by varying amounts; that is, according to the mean residual for each individual earthquake (see discussion in Section 3.1). Averaging over a large set of traveltime residuals that belong to a number of events, and which in turn span a broad range of mean values, seems to induce a tendency of lev-

elling out the median LP-SP differences for moderate short-period residuals.

The similarity of the slopes at large short-period residuals between the whole-mantle and lower-mantle model in Fig. 5 mentioned above is intriguing. It could possibly mean that the dispersion acquired in the lower mantle is still imprinted also in the whole-mantle residuals, but that the corresponding LP-SP differences are shifted to larger positive or negative short-period residuals depending on the event-mean; that is, due to the upper-mantle contribution.

It turns out that indeed, for individual earthquakes, there is a shift of the LP-SP differences of the whole-mantle model mainly along the horizontal axis (as a consequence of the non-zero short-period mean residual), but also to a minor extent along the vertical axis (as a consequence of the typically quite small but non-zero mean of the LP-SP differences). This is shown in Fig. 6 for the two events with the largest negative and positive mean short-period residual (events 2 and 9, respectively). Most important, the median of the distributions of the LP-SP differences show a similarly large slope as those of the whole data set in the lower-mantle model. However, they are shifted towards negative short-period residuals and negative LP-SP differences in case of event 2, and towards positive short-period residuals and positive LP-SP differences in case of event 9. This could indicate that the difference between cross-correlation traveltime residuals measured at different periods for the same seismic station and a given event may be influenced by the upper-mantle structure only to a minor extent. Turned around, the diffraction-induced dispersion seems to be dominated by structure in the lower mantle. If so, the LP-SP differences (i.e. combinations of residuals measured in different frequency bands for a given source-receiver pair) could represent a new datum for tomographic inversions focusing on lower-mantle structure. We will call this datum 'differential-frequency traveltime residual' and demonstrate its similarity between whole-mantle and lower-mantle model in the following. Note that out of the four frequency bands used for our finite-frequency traveltime measurements, six differential-frequency residuals can be constructed.

## 4.1 Imprint of the lower-mantle signal in differential-frequency traveltime residuals

The fact that the differential-frequency residuals indeed do carry substantial information on lower-mantle structure is illustrated in Fig. 7. There, several maps of the seismic heterogeneity in the model, and of *S*-wave traveltime residuals for single-frequency as well as differential-frequency residuals are plotted for the Iceland earthquake (event 5). The leftmost map in the upper row shows $v_s$ variations in model S09-M2-Q at 50 km depth beneath the seismic stations. The subducting slabs all along the west coasts of the Americas and along the Aleutian arc can clearly be recognized as well as the mid-Atlantic ridge. The African and Eurasian continents both show fast velocities at this depth. As we are interested in the traveltime signal related to the lower-mantle structure, we additionally show the $v_s$ perturbation at the ray turning point, but plotted at the location of the respective receiver. Note that these values of seismic heterogeneity thus actually correspond to a different geographic location at depth, half-way between the earthquake and the seismic station, and that the depth increases with epicentral distance from top to bottom of the lower mantle. For the discussion later, we highlight some particular features in the seismic structures at the ray turning point in the following that reappear as characteristic patterns in the maps of traveltime residuals:
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Same as Figs. 5(a) and (b) for subsets of the data from the whole-mantle model belonging to two individual earthquakes (a,b: event 2; c,d: event 9; labels according to Table 1). A shift of the distributions along both axes can be observed in these single-event cases. Events 2 and 9 are the earthquakes with largest negative and positive mean residuals (at 15 s), respectively.</figcaption></figure>

(1) The first feature (F1) is a band of fast material on the left of the map running from the Caribbean through Colombia into the Pacific Ocean.

(2) The second feature (F2) is another band of fast material stretching from Florida through Mexico into the Pacific Ocean with direction roughly towards Hawaii.

(3) The third feature (F3) is the ‘blob’ of strong positive shear wave velocity variations that is separated from Feature F2 by a thin stripe of negative anomalies.

(4) The fourth feature (F4) consists of two bands of seismically slow and fast material starting off in the middle of the North Atlantic, running through Brazil and into the South Atlantic.

(5) The fifth feature (F5) consists of two bands of slow and fast material running from close to the Sea of Okhotsk through all of China.

The above features are also indicated in the bottom left of Fig. 7 by their labels F1–F5. Features F1 and F2 correspond to material that entered the mantle during early stages of the Farallon plate

subduction. The broadening of the bands close to the edges of the maps (i.e. at 90° epicentral distance) results from the horizontal spreading of downwelling material when it hits the CMB.

The other maps in Fig. 7 show traveltime residuals. In the second column from left, single-frequency residuals at 15 s period are plotted for the whole-mantle model at the top and the lower-mantle model at the bottom. As discussed in Section 3.1, the residual map for the whole-mantle model shows strong resemblance to the seismic structure of model S09-M2-Q at shallow depth (i.e. at 50 km depth), except for the effect that all residuals are shifted to positive values due to the event being located in a predominantly slow region. The residual map for the lower-mantle model, in contrast, shows little similarity with the near-surface structure, but much closer resemblance to the heterogeneity at the ray turning point. For example, the bands of fast material corresponding to Features F1 and F2 show up as negative residuals (for Feature F1 the signal is somewhat faint in the Caribbean, but stronger near Ecuador and Peru, and Feature F2 is clearly discernible between Florida and Hawaii running across Mexico). Feature F3 appears as a small
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Comparison of single-frequency and differential-frequency S-wave traveltime residuals in the whole-mantle and the lower-mantle model and their relation to the seismic heterogeneity in model S09-M2-Q. Column 1 (from left): shear wave velocity variations at (top) 50 km depth and (bottom) at the ray-turning point of the S-wave plotted at the location of the respective receiver. Column 2: single-frequency S-wave residuals at 15 s period for event 5 in (top) the whole-mantle and (bottom) the lower-mantle model. Columns 3 and 4: same as column 2 for differential-frequency residuals with combination of period bands 15–22 and 15–51 s, respectively. Single-frequency residuals in the whole-mantle model are clearly related to heterogeneity near the surface, while in the lower-mantle model they bear a strong resemblance to the heterogeneity at the ray-turning point. Note the strong similarity of the differential-frequency traveltime residuals between the whole-mantle and the lower-mantle model.</figcaption></figure>

region of negative residuals (in fact the strongest in the Pacific), which is comparatively thin compared to its corresponding counterpart in the seismic structures. Features F4 and F5 are also discernible in the traveltime map of the lower-mantle model, but mainly only present as positive residuals that are related to the slow part of the features. Interestingly, the signal from Feature F4 is, to a lesser extent, also discernible in the residual map for the whole-mantle model.

In addition to the maps of single-frequency traveltime residuals, Fig. 7 shows maps of differential-frequency residuals for two combinations of periods (15–22.5 and 15–51 s in the third and fourth column from left, respectively). Not surprising, the residuals overall have smaller magnitudes than in the single-frequency maps (note the different colour scales). As expected based on the larger differences between the size of their Fresnel zones, larger values occur in case of the 15–51 s differential-frequency residuals compared to the 15–22.5 s residuals. A second observation is that the geographic pattern of the differential-frequency residuals shows more short-scale variations than in case of the single-frequency measurements. This is related to the fact that these differential traveltimes are influenced by the seismic heterogeneity in a more complex way. Owing to the oscillating nature of single-frequency sensitivity kernels, the difference of two such kernels will oscillate on even shorter length scales (two such differential-frequency sensitivity kernels for the case of a direct S-wave are shown in Fig. 8). The differential-frequency sensitivity kernels can easily be computed—based on the

Born approximation—as a linear combination (i.e. difference) of two classical single-frequency kernels corresponding to the periods involved (i.e. similar to the kernels for differential traveltimes involving different seismic phases as given, for example, in Dahlen et al. (2000) and Hung et al. (2000)):

$$ \delta t_i^{\text{SF}}(T) = \int_{\oplus} K_i(T, \mathbf{r})m(\mathbf{r})d^3\mathbf{r} \qquad (3) $$

$$ \begin{aligned} \delta t_i^{\text{DF}}(T_1, T_2) &= \int_{\oplus} K_i(T_1, \mathbf{r})m(\mathbf{r})d^3\mathbf{r} \\ &\quad - \int_{\oplus} K_i(T_2, \mathbf{r})m(\mathbf{r})d^3\mathbf{r} \\ &= \int_{\oplus} (K_i(T_1, \mathbf{r}) - K_i(T_2, \mathbf{r}))m(\mathbf{r})d^3\mathbf{r} \\ &= \int_{\oplus} \hat{K}_i(T_1, T_2, \mathbf{r})m(\mathbf{r})d^3\mathbf{r} \end{aligned} \qquad (4) $$

with

$$ \hat{K}_i(T_1, T_2, \mathbf{r}) = K_i(T_1, \mathbf{r}) - K_i(T_2, \mathbf{r}). \qquad (5) $$

Here, $\delta t$ denotes the cross-correlation traveltime residual and $T$ the period of the wave; SF and DF refer to ‘single-frequency’ and ‘differential-frequency’, respectively; $K_i$ and $\hat{K}_i$ are the single-frequency and differential-frequency traveltime kernels associated
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Illustration of the short-scale oscillations of differential-frequency sensitivity kernels. Top panels show vertical slices through 3-D S-wave sensitivity kernels for a receiver at 70° epicentral distance and for measurements at (a) 15, (b) 22 and (c) 51 s dominant periods. Bottom panels show resulting differential-frequency kernels (i.e. the difference between two single-frequency kernels) for the combination of period bands (d) 15–22 and (e) 15–51 s.</figcaption></figure>

with receiver *i* for a given earthquake, and *m* denotes the medium properties; that is, the relative perturbation in seismic velocities (e.g. $dlnv_s$); and *r* is the radius vector. The volume integral formally is taken over the entire volume of the Earth $\oplus$, but in tomographic practice typically is restricted to the region in which the amplitude of the finite-frequency sensitivity kernels is significant (e.g. Zaroli *et al.* 2013). The computation of the differential-frequency kernels in equation (5) exploits the linearity of equation (3), which has been verified for velocity perturbations of up to $\pm 10$ per cent by Mercerat & Nolet (2013). The differential-frequency residuals that we construct here are similar in nature to the ‘optimal observables’ recently introduced by Bernauer *et al.* (2014). The construction of the latter involves an optimization step with the goal to enhance sensitivity with respect to certain model parameters such as density. Our differential-frequency measurements on the other hand turn out to highlight a certain region of the mantle without any optimization involved; that is, they are trivial to construct.

Coming back to Fig. 7, the most striking observation is the strong similarity between the differential-frequency residuals of the whole-mantle and the lower-mantle model. Looking at the 15–22.5 s residuals in column 3 first, one can recognize that all five features that we distinguished in the heterogeneity maps show up in the residual maps of both models (for Feature F3 and F5 to a lesser extent). Most important, the almost exact same overall pattern is visible in the top and bottom maps of column 3. The comparison between the single-frequency and differential-frequency residuals for the lower-mantle model reveals the influence of the kernel oscillating on shorter-scales for the latter. For example, Feature F4 is clearly discernible in the bottom map of column 3, however not as a pattern with just positive residuals, but now as two adjacent stripes of both positive and negative residuals.

The strong similarity between differential-frequency residuals of the lower-mantle and the whole-mantle model can also be seen for other combination of periods, for example 15–51 s (column 4

from left in Fig. 7; see also Supporting Information Fig. S3 for the differential-frequency residual maps for the other combinations of periods). For the 15–51 s case, the features identified in the structures at the ray turning point are also imprinted in both maps, but with stronger amplitudes. The only exception is Feature F5, which appears in the bottom map for the lower-mantle model as adjacent positive and negative stripes (stronger than at 15–22.5 s). In the map for the whole-mantle model, however, this feature cannot be seen, but seems to be overprinted by a long band of negative residuals that runs straight away from the event through Russia. Speculating, this band could be related to the S-wave radiation pattern of the earthquake, as it lies roughly at 45° to the nodal planes of the focal mechanism.

Fig. 9 shows the same comparison between the differential-frequency measurements in the whole-mantle and the lower-mantle model for the P-wave residuals. There, one can also see some similarity between the two models in case of the 15–22.5 s maps; that is, the closest lying period bands. However, the differential-frequency residuals in the whole-mantle model seem to be dominated by negative values overall, which degrades the similarity to some extent. The same can be seen in the 15–51 s case, where the tendency towards negative values is even more severe. In addition, the characteristic band of extreme values shows up again that is related to the interaction of the wavefield with structure near the source as discussed in Section 3.1. Still, some imprint of Features F1 and F4 is visible in the map for the whole-mantle model, and signals from the upper mantle are largely cancelled. In the 15–51 s map, the negative residuals related to the subduction of the Pacific plate and the shallow lying remnants of the Farallon slab, which are the dominant features in the 15 s single-frequency map for the whole-mantle model, are only faintly visible. In case of the 15–22.5 s differential-frequency residuals, this upper-mantle signal is almost completely cancelled as well as the signal related to the Tethys subduction.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Same as Fig. 7 (columns 2–4) for *P*-wave traveltime residuals.</figcaption></figure>

4.1.1 *Quantifying the similarity between residuals in the whole-mantle and in the lower-mantle model*

It is clear that a simple visual comparison of maps is unsatisfactory to judge on their similarity between the whole-mantle and the lower-mantle model. To this end, we analysed the similarity between the maps in more detail based on 1) the root-mean-square (RMS) difference of either single-frequency or differential-frequency residuals for all stations, and 2) a ‘direct station-by-station’ (DSS) similarity. We define the DSS similarity as the number of stations that have an absolute difference between the lower-mantle residual and the whole-mantle residual of less than 0.2 s (two times the precision of our cross-correlation traveltime residuals), relative to the total number of stations. The DSS thus neglects minor differences between the differential-frequency residuals and is a scalar number that can vary between 0 (no similarity) and 1 (perfect similarity). In Appendix SOM1 in the Supporting Information, we provide further details including an analysis based on an additional measure of graphical similarity between maps of seismic heterogeneity and maps of traveltime residuals.

The results of the similarity analysis are given in Fig. 10, where the mean and standard deviation of the DSS similarity and RMS differences are shown for the set of 17 earthquakes. In case of single-frequency measurements, the similarity between both models is quite small, with mean values close to 0.15 for all period bands in case of *P*-waves and below 0.1 in case of *S*-waves. As expected from the visual inspection of the maps discussed above, the similarity is substantially larger in case of the differential-frequency measurements with mean values between 0.4 and 0.8 in case of *P*-waves and between 0.3 and 0.7 in case of *S*-waves. The smallest value for both types of waves is obtained for the largest difference in dominant period (15–51 s) and the largest similarity for the smallest difference in dominant period (15–22.5 s). Interestingly, the increase in similarity between the single-frequency and the differential-frequency residuals is larger for the *S*-waves compared to that of the *P*-waves.

The case for differential-frequency traveltimes gets even stronger when looking at the RMS difference of residuals between the two models (Figs 10c and d). Mean RMS differences for the *P*-wave single-frequency residuals all fall within the range from around 0.8 to 1.5 s and those of the *S*-waves within the range from 2.0 to 2.5 s. In contrast, RMS differences are much smaller in case of the differential-frequency residuals, the largest ones being ~0.8 s for *P*-waves and ~1.1 s for *S*-waves at the period combination 15–51 s. The smallest ones are around 0.25 and ~0.3 s, respectively, for the 15–22.5 s period combination.

The above analysis confirms the visual impression from Fig. 7 that there is significant similarity between the differential-frequency measurements of the whole-mantle model and the lower-mantle model. This in turn indicates that differential-frequency residuals (i.e. the diffraction-induced dispersion) indeed are less influenced by the upper-mantle structure of our model than the single-frequency residuals. In case of *S*-waves, they even are dominated by seismic heterogeneity in the lower mantle, as can be seen from the ratio of the standard deviations of dispersion between the lower-mantle and the whole-mantle model in Table 4.

4.2 Potential benefit of differential-frequency
measurements for tomography

We speculate that differential-frequency traveltime residuals may be beneficial when used as data for tomographic inversions that try to illuminate the deep Earth structure. It seems that by taking the difference between cross-correlation measurements in different frequency bands, the signal from shallow mantle structure is (at least partially) removed. In fact, any constant offset coherent between the two measurements, such as for example errors related to event mislocation or origin time that affect real data, will be cancelled. The degree to which the signal of the upper mantle and lithosphere is eliminated depends on the choice of dominant periods
<figure><img src="image_12.png" /><figcaption>**Figure 10.** (a,b) Direct station-by-station similarity and (c,d) root-mean-square difference between traveltime residuals in the whole-mantle and the lower-mantle model. Upper and lower parts of each panel show results for single-frequency and differential-frequency residuals, respectively. Red squares and error bars show the mean and standard deviation, respectively, for the set of 17 earthquakes simulated in this study.</figcaption></figure>

(controlling the length-scale of oscillations of the sensitivity kernels) in combination with the characteristic length scales of structure in the lithosphere and upper mantle. For example, if the structural length scales at shallow depth are larger than the width of the Fresnel zones of both waves there, the resulting contribution to the traveltime will be similar for each single frequency and thus the difference will be close to zero. In essence, for large-scale structures in the upper mantle the differential-frequency measurements are similar in nature to other classical (i.e. ray-theoretical) differential traveltime measurements, such as PP-P/SS-S or PcP-P/ScS-S. The advantage of the new datum, however, is that the waves corresponding to both single-frequency residuals have travelled along exactly the same ray-theoretical path (just having Fresnel zones of different width), which is not true for the classical differential traveltimes. It remains to be tested whether the sensitivity of differential-frequency traveltime residuals can be represented with sufficient accuracy using ray-theoretical kernels (excluding the near- and mid-field terms and neglecting the P-wave cross-dependence)

or whether one has to resort to full waveform based adjoint kernels. Note that for the former case this would involve a very fine sampling of the differential-frequency sensitivity kernels in order to correctly capture their short-scale oscillations.

In this study, we have provided some first evidence that differential-frequency measurements may bear advantages for tomographic studies. Our findings are based on synthetic seismic structures, which are derived from a geodynamic model of mantle flow. In earlier work, we have shown that the seismic heterogeneity in our model is very well compatible with body-wave data and tomographic models in a statistical sense (Schuberth et al. 2009a,b, 2012). In particular, there is very good agreement in the spectral characteristics of the synthetic mantle structures with those imaged by tomography, which means that the characteristic length scales are earth-like. However, lithosphere structure may be too simplistic in our model, as it is represented only as a pure thermal boundary layer. It thus remains to be seen whether the addition of more realistic structures in the upper 100–150 km (including the ocean-continent dichotomy, etc.), for example from models like 3SMAC (Nataf & Ricard 1996) or Litho1.0 (Pasyanos et al. 2014), will lead to the same conclusion.

It is also clear that this study just provides a starting point and that a formal inversion is needed to quantify the potential gain of using differential-frequency traveltime residuals in a realistic tomographic scenario. With a similar intention, Zaroli et al. (2014) recently demonstrated that errors in tomographic inversions due to uncertainties in earthquake parameters (e.g. hypocentre location and origin time) can be reduced by using suitable receiver pairs. Based on the results of this study, one question that we wish to address in future will be whether the differential-frequency residuals can also help in reducing the mislocation-related bias. Another question is whether these data would help in treating the signal of the crust. Crustal corrections typically pose considerable problems in tomographic studies, as the structure in the shallowest part of the Earth is not perfectly known (e.g. Ferreira et al. 2010). This is true not only for ray-theoretical body-wave as well as full waveform surface wave inversions (e.g. Lekić et al. 2010), but in particular also for tomographic studies based on finite-frequency theory (e.g. Yang & Shen 2006; Ritsema et al. 2009). There, special care needs to be taken to mimic the waveform distortion due to crustal reverberations when computing the reference waveform for the cross-correlation measurements (Zaroli et al. 2010). This processing step could potentially introduce systematic errors and would thus preferably be avoided.

## 5 CONCLUSIONS

In this study, we have investigated the strength of diffraction-induced traveltime dispersion of seismic body waves in mantle models with isotropic and perfectly elastic seismic heterogeneity. The 3-D distribution of seismic velocities and densities is based on the temperature field of an isochemical MCM. To specifically quantify the dispersion related to the lower mantle, we investigated a model in which heterogeneity was removed in the upper 800 km, in addition to the original model in which heterogeneity exists in the entire mantle. 3-D global wave propagation in the two models was simulated using a spectral-element method deliberately excluding effects of intrinsic attenuation as well as 3-D structure in the crust. Using the waveforms of direct P- and S-waves in the synthetic seismograms, we have measured traveltime residuals in four frequency bands using a cross-correlation technique. The frequencies used
here are at the long period end of the usable spectrum for body waves, in particular for P-waves. The restriction to periods above 10 s is related to the intense computational requirements of wave propagation simulations on a global scale. With the ever growing capabilities of high-performance computing infrastructures, it will be possible in near future to increase the frequencies in the simulations and to cover a larger range of the teleseismic body-wave spectrum. Nevertheless, the synthetic multifrequency measurements obtained here allow us—for the first time—to characterize in a relevant part of the frequency range the wavefield effects to be expected in a mantle with earth-like structural length scales and realistic magnitudes of seismic heterogeneity.

The most important result is that our synthetic data do indeed show significant diffraction-induced dispersion between the longest and shortest period, in both the whole-mantle and the lower-mantle model. The following general trend can be observed: Differences between long-period and short-period residuals tend to increase with increasing short-period signal. The fact that these differences tend to be positive for negative short-period residuals and negative for positive ones indicates ubiquitous wave-front healing. At first sight, the dispersion signal relative to the short-period residual appears to be larger in the lower-mantle model compared to the whole-mantle model. The smaller relative dispersion in the whole-mantle model, however, is due to non-zero mean short-period residuals for individual earthquakes related to the near-surface effects. For data from a single event, the variation of the diffraction-induced dispersion with the short-period residual is similar between both models if the event-mean is removed.

The dispersion, in the way it is defined here as the difference between traveltime residuals measured in different frequency bands for a given station-event pair, can also be regarded as a new seismic datum: a differential-frequency time residual. We have quantified the similarity of the source-receiver specific differential-frequency residuals between the whole-mantle and the lower-mantle model. A second important outcome of our study is that the dispersion is not only similar in a statistical sense, but also for the individual differential-frequency residuals. This leads us to conclude that this datum is to some extent insensitive to upper-mantle structure. The similarity of the dispersion between our two models is more pronounced in case of S-waves and varies among the different combinations of frequency bands (i.e. different combinations of wavelengths/Fresnel zones). The extent to which differential-frequency residuals are insensitive to upper-mantle structure thus seems to depend on the relation between the wavelengths involved and the characteristic length scales of heterogeneities at shallow depth. Differential-frequency traveltime residuals may therefore potentially prove useful to learn more about the characteristic length scales of $v_p$ and $v_s$ variations in the lithosphere and upper mantle. Using differential-frequency data may also help to improve tomographic models of the lower mantle, as they might result in elimination of systematic errors and precondition the inversions such as to reduce influence of upper mantle for certain length scales. Formal inversions on synthetic test cases will be necessary to elucidate whether these data can indeed be beneficial for tomographic studies. In addition, future studies need to investigate the degree of contamination of the differential-frequency residuals by the seismic signature of the crust and its impact on our conclusions as well as the effects of a potential frequency-dependence of the source-time function.

Overall, the dispersion of traveltime residuals introduced by virtue of diffraction at 3-D heterogeneities (which here are derived from the temperature field of a mantle flow model) is quite

substantial and is therefore of relevance to seismic studies. This notion is based on our findings of rather large standard deviations of the LP-SP differences compared to the average magnitude of the single-frequency residuals. For the longest and shortest periods considered here (i.e. 51 and 15 s, respectively), standard deviations of dispersion are 0.6 s for P-waves and 1.0 s for S-waves. This represents a large fraction of the total short-period residual. The standard deviations of single-frequency residuals at 15 s are 1.0 and 2.8 s for P- and S-waves, respectively (note that these values agree with the inferred mantle contribution to the total variance of observed traveltime residuals; Bolton & Masters 2001; Schuberth et al. 2012). Relative to the short-period residuals, P-waves thus show stronger diffraction-induced dispersion than S-waves, which is expected from their larger Fresnel zone at the same frequency.

The strength of dispersion found here is significant not only with respect to the magnitude of the residuals themselves, but also with respect to the dispersion observed in real data. Zaroli et al. (2010) have measured cross-correlation traveltime residuals in the global seismic network for S, SS and ScS waves at the same frequencies as used in the present study. They also found values of dispersion for the S-waves of around 1–2 s, similar to those found here. We note again that we did not include crustal effects in our study, which might explain the slightly smaller dispersion in our case. If the information contained in this dispersion signal is fully exploited through multifrequency inversions, Earth’s deep interior can likely be imaged at unprecedented resolution and accuracy. In particular, it will hopefully be possible to better constrain the gradients of velocity anomalies as well as their magnitude. To this end, it will be necessary to reduce (at best eliminate) uncertainties in the data related to event mislocation, inexact origin time and effects from crustal corrections using imperfect models. All these effects actually introduce errors that are of same order as the average traveltime residual itself (e.g. Bolton & Masters 2001). Zaroli et al. (2014) have recently explored new ways to reduce the influence of event mislocation on tomographic images, but further efforts in this direction are necessary.

## ACKNOWLEDGEMENTS

We want to thank Diego Mercerat and Lorenzo Colli for fruitful and stimulating discussions and two anonymous reviewers for their careful comments that helped to improve the manuscript. BSAS received support from a Marie Curie Intra European Fellowship within the 7th European Community Framework Programme [FP7/2007-2013] under grant agreement no. 235861. GN and CZ received support from the European Research Council (ERC Advanced grant 226837). CZ also thanks the SEISGLOB ANR 2011 Blanc SIMI 5-6-016-01. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, www.lrz.de). We also thank the DEISA Consortium (www.deisa.eu), co-funded through the EU FP6 project RI-031513 and the FP7 project RI-222919, for support within the DEISA Extreme Computing Initiative. The maps appearing in some of the figures were generated with *The Generic Mapping Tools* (GMT) (Wessel et al. 2013).

## REFERENCES

Agranier, A., Blicherttoft, J., Graham, D., Debaille, V., Schiano, P. & Albarede, F., 2005. The spectra of isotopic heterogeneities along the mid-Atlantic Ridge, *Earth planet. Sci. Lett.*, **238**(1–2), 96–109.
Allègre, C.J. & Turcotte, D.L., 1986. Implications of a two-component marble-cake mantle, *Nature*, **323**(6084), 123-127.

Baig, A.M. & Dahlen, F.A., 2004. Traveltime biases in random media and the S-wave discrepancy, *Geophys. J. Int.*, **158**(3), 922-938.

Baig, A.M., Dahlen, F.A. & Hung, S.-H., 2003. Traveltimes of waves in three-dimensional random media, *Geophys. J. Int.*, **153**(2), 467-482.

Bernauer, M., Fichtner, A. & Igel, H., 2014. Optimal observables for multiparameter seismic tomography, *Geophys. J. Int.*, **198**(2), 1241-1254.

Bolton, H. & Masters, G., 2001. Travel times of P and S from the global digital seismic networks: implications for the relative variation of P and S velocity in the mantle, *J. Geophys. Res.*, **106**(B7), 13 527-13 540.

Bull, A., McNamara, A. & Ritsema, J., 2009. Synthetic tomography of plume clusters and thermochemical piles, *Earth planet. Sci. Lett.*, **278**(3-4), 152-162.

Bunge, H.-P., 2005. Low plume excess temperature and high core heat flux inferred from non-adiabatic geotherms in internally heated mantle circulation models, *Phys. Earth planet. Inter.*, **153**(1-3), 3-10.

Bunge, H.-P. & Baumgardner, J., 1995. Mantle convection modeling on parallel virtual machines, *Comput. Phys.*, **9**, 207-215.

Bunge, H.-P. & Richards, M., 1996. The origin of large-scale structure in mantle convection: Effects of plate motions and viscosity stratification, *Geophys. Res. Lett.*, **23**, 2987-2990.

Bunge, H.-P., Richards, M. & Baumgardner, J., 1996. The effect of depth-dependent viscosity on the planform of mantle convection, *Nature*, **379**, 436-438.

Bunge, H.-P., Richards, M. & Baumgardner, J., 1997. A sensitivity study of 3D-spherical mantle convection at 10exp8 Rayleigh number: effects of depth-dependent viscosity, heating mode and an endothermic phase change, *J. geophys. Res.*, **102**, 11 991-12 007.

Christensen, U.R. & Hofmann, A.W., 1994. Segregation of subducted oceanic crust in the convecting mantle, *J. geophys. Res.*, **99**(B10), 19 867-19 884.

Cormier, V.F., 1989. Slab diffraction of S waves, *J. geophys. Res.*, **94**(B3), 3006-3024.

Cottaar, S. & Romanowicz, B., 2012. An unusually large ULVZ at the base of the mantle near Hawaii, *Earth planet. Sci. Lett.*, **355-356**(0), 213-222.

Dahlen, F.A. & Baig, A.M., 2002. Fréchet kernels for body-wave amplitudes, *Geophys. J. Int.*, **150**(2), 440-466.

Dahlen, F.A., Hung, S.H. & Nolet, G., 2000. Fréchet kernels for finite-frequency traveltimes-I. Theory, *Geophys. J. Int.*, **141**(1), 157-174.

Davaille, A., 1999. Simultaneous generation of hotspots and superswells by convection in a heterogenous planetary mantle, *Nature*, **402**(6763), 756-760.

Davaille, A., Girard, F. & Le Bars, M., 2002. How to anchor hotspots in a convecting mantle?, *Earth planet. Sci. Lett.*, **203**(2), 621-634.

Davies, D.R., Goes, S., Davies, J., Schuberth, B., Bunge, H.-P. & Ritsema, J., 2012. Reconciling dynamic and seismic models of Earth's lower mantle: The dominant role of thermal heterogeneity, *Earth Planet. Sci. Lett.*, **353-354**(0), 253-269.

Davies, H.J. & Bunge, H.-P., 2001. Seismically "fast" geodynamic mantle models, *Geophys. Res. Lett.*, **28**(1), 73-76.

de Wit, R. W.L., Trampert, J. & van der Hilst, R.D., 2012. Toward quantifying uncertainty in travel time tomography using the null-space shuttle, *J. geophys. Res.*, **117**(B3), B03301, doi:10.1029/2011JB008754.

Dziewonski, A.M. & Anderson, D.L., 1981. Preliminary reference Earth model, *Phys. Earth planet. Inter.*, **25**, 297-356.

Ekström, G., Nettles, M. & Dziewonski, A.M., 2012. The global CMT project 2004-2010: centroid-moment tensors for 13,017 earthquakes, *Phys. Earth planet. Inter.*, **200-201**, 1-9.

Ferreira, A. M.G., Woodhouse, J.H., Visser, K. & Trampert, J., 2010. On the robustness of global radially anisotropic surface wave tomography, *J. geophys. Res.*, **115**(B4), B04313, doi:10.1029/2009JB006716.

French, S.W. & Romanowicz, B.A., 2014. Whole-mantle radially anisotropic shear velocity structure from spectral-element waveform tomography, *Geophys. J. Int.*, **199**(3), 1303-1327.

Gilbert, F. & Dziewonski, A., 1975. An application of normal mode theory to the retrieval of structural parameters and source mechanisms for seismic spectra, *Phil. Trans. R. Soc. Lond.*, **A**, **278**, 187-269.

Hager, B.H., 1984. Subducted slabs and the geoid: Constraints on mantle rheology and flow, *J. geophys. Res.*, **89**, 6003-6015.

Hansen, U. & Yuen, D.A., 2000. Extended-Boussinesq thermal-chemical convection with moving heat sources and variable viscosity, *Earth planet. Sci. Lett.*, **176**(3-4), 401-411.

Hong, T.-K., Wu, R.-S. & Kennett, B., 2005. Stochastic features of scattering, *Phys. Earth planet. Inter.*, **148**(2-4), 131-148.

Hosseini, K. & Sigloch, K., 2015. Multifrequency measurements of core-diffracted P-waves (Pdiff) for global waveform tomography, *Geophys. J. Int.*, **203**(1), 506-521.

Hung, S.-H., Dahlen, F. & Nolet, G., 2000. Fréchet kernels for finite-frequency traveltimes-II. Examples, *Geophys. J. Int.*, **141**(1), 175-203.

Hung, S.-H., Dahlen, F. & Nolet, G., 2001. Wavefront healing: a banana-doughnut perspective, *Geophys. J. Int.*, **146**(2), 289-312.

Igel, H. & Gudmundsson, O., 1997. Frequency-dependent effects on travel times and waveforms of long-period S and SS waves, *Phys. Earth planet. Inter.*, **104**(1-3), 229-246.

Igel, H. & Ita, J., 1997. The effects of subduction zones on teleseismic SH waves: A numerical study, in *Upper Mantle Heterogeneities from Active and Passive Seismology*, pp. 333-341, ed. Fuchs, K., Kluwer Academic Publishers.

Igel, H., Nissen-Meyer, T. & Jahnke, G., 2002. Wave propagation in 3D spherical sections: effects of subduction zones, *Phys. Earth planet. Inter.*, **132**(1-3), 219-234.

Jellinek, A.M. & Manga, M., 2004. Links between long-lived hot spots, mantle plumes, D", and plate tectonics, *Rev. Geophys.*, **42**, RG3002, doi:10.1029/2003RG000144.

Kellogg, L.H., Hager, B.H. & van der Hilst, R.D., 1999. Compositional stratification in the deep mantle, *Science*, **283**(5409), 1881-1884.

Kennett, B.L. & Bunge, H.-P., 2008. *Geophysical Continua*, Cambridge Univ. Press.

Kleiner, B. & Graedel, T.E., 1980. Exploratory data analysis in the geophysical sciences, *Rev. Geophys.*, **18**(3), 699-717.

Komatsitsch, D. & Tromp, J., 2002a. Spectral-element simulations of global seismic wave propagation-I. Validation, *Geophys. J. Int.*, **149**, 390-412.

Komatsitsch, D. & Tromp, J., 2002b. Spectral-element simulations of global seismic wave propagation-II. Three-dimensional models, oceans, rotation and self-gravititation, *Geophys. J. Int.*, **150**(1), 303-318.

Lekić, V., Panning, M. & Romanowicz, B., 2010. A simple method for improving crustal corrections in waveform tomography, *Geophys. J. Int.*, **182**(1), 265-278.

Li, M. & McNamara, A.K., 2013. The difficulty for subducted oceanic crust to accumulate at the Earth's core-mantle boundary, *J. geophys. Res.*, **118**(4), 1807-1816.

Li, Y., Deschamps, F. & Tackley, P.J., 2014. The stability and structure of primordial reservoirs in the lower mantle: insights from models of thermochemical convection in three-dimensional spherical geometry, *Geophys. J. Int.*, **199**(2), 914-930.

Lithgow-Bertelloni, C. & Richards, M.A., 1998. The dynamics of Cenozoic and Mesozoic plate motions, *Rev. Geophys.*, **36**(1), 27-78.

Maggi, A., Tape, C., Chen, M., Chao, D. & Tromp, J., 2009. An automated time-window selection algorithm for seismic tomography, *Geophys. J. Int.*, **178**(1), 257-281.

Malcolm, A.E. & Trampert, J., 2011. Tomographic errors from wave front healing: more than just a fast bias, *Geophys. J. Int.*, **185**(1), 385-402.

Marquering, H., Dahlen, F.A. & Nolet, G., 1999. Three-dimensional sensitivity kernels for finite-frequency traveltimes: the banana-doughnut paradox, *Geophys. J. Int.*, **137**, 805-815.

McNamara, A. & Zhong, S., 2005. Thermochemical structures beneath Africa and the Pacific Ocean, *Nature*, **437**, 1136-1139.

Mercerat, E. & Nolet, G., 2012. Comparison of ray-based and adjoint-based sensitivity kernels for body-wave seismic tomography, *Geophys. Res. Lett.*, **39**, L12301, doi:10.1029/2012GL052002.

Mercerat, E. & Nolet, G., 2013. On the linearity of cross-correlation delay times in finite-frequency tomography, *Geophys. J. Int.*, **192**, 681-687.

Mercerat, E., Nolet, G. & Zaroli, C., 2014. Cross-borehole tomography with correlation delay times, *Geophys. J. Int.*, **79**, R1-R12.
Montague, N.L. & Kellogg, L.H., 2000. Numerical models of a dense layer at the base of the mantle and implications for the geodynamics of D", *J. geophys. Res.*, **105**(B5), 11 101–11 114.

Mosca, I., Cobden, L., Deuss, A., Ritsema, J. & Trampert, J., 2012. Seismic and mineralogical structures of the lower mantle from probabilistic tomography, *J. geophys. Res.*, **117**(B6), doi:10.1029/2011JB008851.

Müller, G., Roth, M. & Korn, M., 1992. Seismic-wave traveltimes in random media, *Geophys. J. Int.*, **110**(1), 29–41.

Nakagawa, T. & Tackley, P.J., 2004. Thermo-chemical structure in the mantle arising from a three-component convective system and implications for geochemistry, *Phys. Earth planet. Inter.*, **146**(1–2), 125–138.

Nataf, H.-C. & Richard, Y., 1996. 3SMAC: an a priori tomographic model of the upper mantle based on geophysical modeling, *Phys. Earth planet. Inter.*, **95**(1–2), 101–122.

Nolet, G., 1991. Imaging the deep Earth: technical possibilities and theoretical limitations, in *Proc. XXIIIth Assembly ESC, Barcelona 1990*, pp. 107–115, ed. A. Roca.

Nolet, G., 2008. *A Breviary of Seismic Tomography*, Cambridge Univ. Press.

Nolet, G. & Dahlen, F.A., 2000. Wave front healing and the evolution of seismic delay times, *J. geophys. Res.*, **105**(B8), 19 043–19 054.

Nolet, G. & Moser, T.-J., 1993. Teleseismic delay times in a 3-D Earth and a new look at the S discrepancy, *Geophys. J. Int.*, **114**(1), 185–195.

Obrebski, M., Allen, R.M., Xue, M. & Hung, S.-H., 2010. Slab-plume interaction beneath the Pacific Northwest, *Geophys. Res. Lett.*, **37**(14), L14305, doi:10.1029/2010GL043489.

Oeser, J., Bunge, H.P. & Mohr, M., 2006. Cluster design in the earth sciences–Tethys, in *High Performance Computing And Communications, Proceedings*, vol. 4208, pp. 31–40, eds Gerndt, M. & Kranzlmüller, D., Springer.

Pasyanos, M.E., Masters, T.G., Laske, G. & Ma, Z., 2014. Litho1.0: an updated crust and lithospheric model of the Earth, *J. geophys. Res.*, **119**(3), 2153–2173.

Paulson, A., Zhong, S.J. & Wahr, J., 2007. Limitations on the inversion for mantle viscosity from postglacial rebound, *Geophys. J. Int.*, **168**(3), 1195–1209.

Piazzoni, A.S., Steinle-Neumann, G., Bunge, H.-P. & Dolejš, D., 2007. A mineralogical model for density and elasticity of the Earth's mantle, *Geochem. Geophys. Geosyst.*, **8**(11), Q11010, doi:10.1029/2007GC001697.

Ricard, Y., Spada, G. & Sabadini, R., 1993. Polar wandering of a dynamic Earth, *Geophys. J. Int.*, **113**(2), 284–298.

Ricard, Y., Durand, S., Montagner, J.-P. & Chambat, F., 2014. Is there seismic attenuation in the mantle?, *Earth planet. Sci. Lett.*, **388**, 257–264.

Ritsema, J., McNamara, A.K. & Bull, A.L., 2007. Tomographic filtering of geodynamic models: Implications for model interpretation and large-scale mantle structure, *J. geophys. Res.*, **112**(B1), B01303, doi:10.1029/2006JB004566.

Ritsema, J., van Heijst, H.J., Woodhouse, J.H. & Deuss, A., 2009. Long-period body wave traveltimes through the crust: implication for crustal corrections and seismic tomography, *Geophys. J. Int.*, **179**(2), 1255–1261.

Ritsema, J., Deuss, A., van Heijst, H.J. & Woodhouse, J.H., 2011. S40RTS: a degree-40 shear-velocity model for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements, *Geophys. J. Int.*, **184**(3), 1223–1236.

Romanowicz, B. & Mitchell, B., 2007. Deep Earth structure – Q in the Earth from crust to core, in *Treatise of Geophysics*, Vol. 1, pp. 731–774, ed. Schubert, G., Elsevier.

Roth, M., Müller, G. & Snieder, R., 1993. Velocity shift in random media, *Geophys. J. Int.*, **115**(2), 552–563.

Schaber, K., Bunge, H.-P., Schuberth, B.S.A., Malservisi, R. & Horbach, A., 2009. Stability of the rotation axis in high-resolution mantle circulation models: weak polar wander despite strong core heating, *Geochem. Geophys. Geosyst.*, **10**, Q11W04, doi:10.1029/2009GC002541.

Schuberth, B.S.A., Bunge, H.-P. & Ritsema, J., 2009a. Tomographic filtering of high-resolution mantle circulation models: can seismic heterogeneity

be explained by temperature alone?, *Geochem. Geophys. Geosyst.*, **10**(5), Q05W03, doi:10.1029/2009GC002401.

Schuberth, B.S.A., Bunge, H.-P., Steinle-Neumann, G., Moder, C. & Oeser, J., 2009b. Thermal versus elastic heterogeneity in high-resolution mantle circulation models with pyrolite composition: high plume excess temperatures in the lowermost mantle, *Geochem. Geophys. Geosyst.*, **10**(1), Q01W01, doi:10.1029/2008GC002235.

Schuberth, B.S.A., Zaroli, C. & Nolet, G., 2012. Synthetic seismograms for a synthetic Earth: long-period P- and S-wave traveltime variations can be explained by temperature alone, *Geophys. J. Int.*, **188**(3), 1393–1412.

Shapiro, S.A., Schwarz, R. & Gold, N., 1996. The effect of random isotropic inhomogeneities on the phase velocity of seismic waves, *Geophys. J. Int.*, **127**(3), 783–794.

Sigloch, K. & Nolet, G., 2006. Measuring finite-frequency body-wave amplitudes and traveltimes, *Geophys. J. Int.*, **167**(1), 271–287.

Silver, P.G. & Chan, W.W., 1986. Observations of body wave multipathing from broadband seismograms: evidence for lower mantle slab penetration beneath the Sea of Okhotsk, *J. geophys. Res.*, **91**(B14), 13 787–13 802.

Simmons, N.A., Forte, A.M., Boschi, L. & Grand, S.P., 2010. GyPSuM: A joint tomographic model of mantle density and seismic wave speeds, *J. geophys. Res.*, **115**(B12), B12310, doi:10.1029/2010JB007631.

Spetzler, J. & Snieder, R., 2001. The effect of small-scale heterogeneity on the arrival time of waves, *Geophys. J. Int.*, **145**(3), 786–796.

Stegman, D.R., Richards, M.A. & Baumgardner, J.R., 2002. Effects of depth-dependent viscosity and plate motions on maintaining a relatively uniform mid-ocean ridge basalt reservoir in whole mantle flow, *J. geophys. Res.*, **107**(B6), 2116, doi:10.1029/2001JB000192.

Stixrude, L. & Lithgow-Bertelloni, C., 2005. Thermodynamics of mantle minerals I. Physical properties, *Geophys. J. Int.*, **162**(2), 610–632.

Stixrude, L. & Lithgow-Bertelloni, C., 2007. Influence of phase transformations on lateral heterogeneity and dynamics in Earth's mantle, *Earth planet. Sci. Lett.*, **263**, 45–55.

Stixrude, L. & Lithgow-Bertelloni, C., 2011. Thermodynamics of mantle minerals – II. Phase equilibria, *Geophys. J. Int.*, **184**(3), 1180–1213.

Tackley, P.J., 2000. Mantle convection and plate tectonics: Toward an integrated physical and chemical theory, *Science*, **288**(5473), 2002–2007.

Tan, E. & Gurnis, M., 2007. Compressible thermochemical convection and application to lower mantle structures, *J. geophys. Res.*, **112**(B6), B06304, doi:10.1029/2006JB004505.

Tilmann, F.J., McKenzie, D. & Priestley, K.F., 1998. P and S wave scattering from mantle plumes, *J. geophys. Res.*, **103**(B9), 21 145–21 163.

To, A. & Romanowicz, B., 2009. Finite frequency effects on global S diffracted traveltimes, *Geophys. J. Int.*, **179**(3), 1645–1657.

Trampert, J., Fichtner, A. & Ritsema, J., 2013. Resolution tests revisited: the power of random numbers, *Geophys. J. Int.*, **192**(2), 676–680.

Vidale, J., 1987. Waveform effects of a high-velocity, subducted slab, *Geo-phys. Res. Lett.*, **14**(5), 542–545.

Wessel, P., Smith, W. H.F., Scharroo, R., Luis, J. & Wobbe, F., 2013. Generic mapping tools: Improved version released, *EOS, Trans. Am. geophys. Un.*, **94**(45), 409–410.

Wielandt, E., 1987. On the validity of the ray approximation for interpreting delay times, in *Seismic Tomography*, Vol. 5, pp. 85–98, ed. Nolet, G., Springer.

Wu, R.-S. & Aki, K., 1988. Introduction: seismic wave scattering in three-dimensionally heterogeneous Earth, *Pure appl. Geophys.*, **128**(1–2), 1–6.

Yang, T. & Shen, Y., 2006. Frequency-dependent crustal correction for finite-frequency seismic tomography, *Bull. seism. Soc. Am.*, **96**(6), 2441–2448.

Zaroli, C., Debayle, E. & Sambridge, M., 2010. Frequency-dependent effects on global S-wave traveltimes: wavefront-healing, scattering and attenuation, *Geophys. J. Int.*, **182**(2), 1025–1042.

Zaroli, C., Sambridge, M., Lévêque, J.-J., Debayle, E. & Nolet, G., 2013. An objective rationale for the choice of regularisation parameter with application to global multiple-frequency S-wave tomography, *Solid Earth*, **4**(2), 357–371.

Zaroli, C., Lévêque, J.-J., Schuberth, B.S.A., Duputel, Z. & Nolet, G., 2014. Global S-wave tomography using receiver pairs: an alternative to get rid of earthquake mislocation, *Geophys. J. Int.*, **199**(2), 1043–1057.
## SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this paper:

**Figure S1.** Same as Fig. 2 for an earthquake in Southern Iran (event 2 in Table 1). Traveltime residuals are mostly negative in consequence of the event being located within a broad region of fast seismic velocities related to the convergence of Africa/Arabia with Eurasia and the closure of the Tethys Ocean.

**Figure S2.** Same as Fig. S1 for the lower-mantle model S09-M2-Q-LM (3-D heterogeneity below 800 km depth only). The magnitude of traveltime residuals is lower compared to Fig. S1 and positive and negative residuals are more balanced.

**Figure S3.** Same as Fig. 7 (columns 3 and 4) for the period combinations 15–34, 22.5–34, 22.5–51 and 34–51 s.

**Figure S4.** Comprehensive comparison between S-wave traveltime residuals in the whole-mantle and the lower-mantle model and between traveltime residuals in each model and the seismic heterogeneity. Upper and lower parts of each panel show results for single-frequency and differential-frequency residuals, respectively. Coloured symbols and error bars show the mean and standard deviation, respectively, for the set of 17 earthquakes simulated in this study. (a) Root-mean-square difference, (b) direct station-by-station similarity, (c) correlation and (d) graphical pixel-by-pixel similarity between traveltime residuals in the whole-mantle and the lower-mantle model. (e) Correlation and (f) graphical pixel-by-pixel similarity between traveltime residuals and S-wave velocity heterogeneity. In panels (e) and (f), symbols depict comparisons between (red) residuals in the whole-mantle model and heterogeneity at 50 km depth, (green) same for residuals in the lower-mantle model, (blue) residuals in the whole-mantle model and heterogeneity at the ray turning point and (magenta) same for residuals in the lower-mantle model.

**Figure S5.** Frequency dependence of the standard deviation of traveltime residuals. The scaled median average deviation (SMAD, defined in Section 2.1 of the main text) is employed as measure of scale, which is a robust (i.e. unbiased) estimator for the standard deviation (Kleiner & Graedel 1980). (a,d) Standard deviation of the residuals plotted as a function of ray-turning depth for (blue) P-waves and (red) S-waves. Solid, dashed, dash-dotted and dotted lines correspond to measurements at 15, 22.5, 34 and 51 s dominant period, respectively, in (top) the whole-mantle model and (bottom) the lower-mantle model. (b,c,d,f) Relative change of the SMAD at 22.5, 34 and 51 s period with respect to that at 15 s period (dashed, dash-dotted and dotted lines, respectively).

**Figure S6.** Top maps: Same as Fig. 2 for an earthquake at the central mid-Atlantic ridge (event 8 in Table 1) using modified versions of the whole-mantle model S09-M2-Q: (left) 3-D heterogeneity restricted to $v_s$ (i.e. 1-D profiles for $v_p$ and density are used) and only the P-wave residuals are shown; (right) 3-D heterogeneity restricted to $v_p$ and only the S-wave residuals are shown. The effects of the P-wave cross-dependence can be seen in case of the $v_s$-only model. In contrast, S-wave residuals show very little dependence on variations in compressional velocity, as expected. Bottom: Same as Fig. S5a for event 8 in (left) the $v_s$-only and (right) the $v_p$-only model. (http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggv389/-/DC1).

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.