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

HAL
open science

# Synthetic seismograms for a synthetic Earth: long-period P- and S-wave traveltime variations can be explained by temperature alone

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

► To cite this version:

Bernhard S. A. Schuberth, Christophe Zaroli, Guust Nolet. Synthetic seismograms for a synthetic Earth: long-period P- and S-wave traveltime variations can be explained by temperature alone. Geophysical Journal International, 2012, 188, pp.1393-1412. hal-00704470

HAL Id: hal-00704470

https://hal.science/hal-00704470

Submitted on 7 Mar 2022

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_2.png" /></figure>

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

# Synthetic seismograms for a synthetic Earth: long-period *P*- and *S*-wave traveltime variations can be explained by temperature alone

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

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, 06560 Valbonne, France. E-mail: mail@bernhard-schuberth.de

Accepted 2011 December 10. Received 2011 November 18; in original form 2011 September 6

## SUMMARY

Current interpretations of seismic observations typically argue for significant chemical heterogeneity being present in the two large low shear velocity provinces under Africa and the Pacific. Recently, however, it has been suggested that large lateral temperature variations in the lowermost mantle resulting from a strong thermal gradient across D'' may provide an alternative explanation. In case of a high heat flux from the core into the mantle, the magnitude of shear wave velocity variations in tomographic models can be reconciled with isochemical whole mantle flow and a pyrolite composition. So far, the hypothesis of strong core heating has been tested in a consistent manner only against tomographic S-wave velocity models, but not against P-wave velocity models. Here, we explore a new approach to assess geodynamic models and test the assumption of isochemical whole mantle flow with strong core heating directly against the statistics of observed traveltime variations of both P and S waves. Using a spectral element method, we simulate 3-D global wave propagation for periods down to 10 s in synthetic 3-D elastic structures derived from a geodynamic model. Seismic heterogeneity is predicted by converting the temperature field of a high-resolution mantle circulation model (MCM) into seismic velocities using thermodynamic models of mantle mineralogy. Being based on forward modelling only, this approach avoids the problems of limited resolution and non-uniqueness inherent in tomographic inversions while taking all possible finite-frequency effects into account. Capturing the correct physics of wave propagation allows for a consistent test of the assumption of high core heat flow against seismic data.

The statistics of long-period body wave traveltime observations show a markedly different behaviour for *P* and *S* waves: the standard deviation of *P*-wave delay times stays almost constant with turning depth, whereas that of the *S*-wave delay times increases strongly throughout the mantle. Surprisingly, synthetic traveltime variations computed for the isochemical MCM reproduce these different trends. This is not expected from a ray-theoretical point of view and highlights the importance of finite-frequency effects. Most importantly, the large lateral temperature variations in the lower mantle related to strong core heating are able to explain most of the standard deviation of observed *P*- and *S*-wave delay times. This is a strong indication that seismic heterogeneity in the lower mantle is likely dominated by thermal variations on the length scales relevant for long-period body waves.

**Key words:** Heat flow; Composition of the mantle; Body waves; Computational seismology; Wave scattering and diffraction; Wave propagation.

## 1 INTRODUCTION

Long-standing questions in the study of Earth's deep interior are about the origin of seismic heterogeneity and the nature of flow in

*Now at: Institute de Physique du Globe de Strasbourg (UMR 7516 CNRS, Université de Strasbourg/EOST), 5 rue René Descartes, 67084 Strasbourg Cedex, France.

the mantle. Understanding the dynamic behaviour is important as the flow drives plate tectonics and controls the way the Earth looses its heat. To improve conceptual models of mantle flow, the major challenges in seismology today are to efficiently mine the wealth of information contained in seismic waveforms and to constrain the relative contributions of thermal anomalies and compositional variations to the observed seismic heterogeneity. One specific question is the nature of the two large low shear velocity zones under Africa
and the Pacific consistently mapped out in tomographic models of the mantle (e.g. Dziewonski et al. 1977; Su et al. 1994; Li & Romanowicz 1996; Grand et al. 1997; van der Hilst et al. 1997; Su & Dziewonski 1997; Kennett et al. 1998; Masters et al. 2000; Ritsema et al. 2004; Montelli et al. 2006; Panning & Romanowicz 2006; Simmons et al. 2007; Houser et al. 2008; Kustowski et al. 2008; Li et al. 2008; Simmons et al. 2010). It has been proposed that these ‘superplumes’ are chemically distinct from the rest of the mantle, as they have sharp vertical boundaries and show anticorrelated variations in shear ($v_s$) and bulk sound velocity ($v_\phi$) as well as evidence for an increased density (van der Hilst & Karason 1999; Masters et al. 2000; Ishii & Tromp 2001; Wen et al. 2001; Ni et al. 2002; Ritsema & van Heijst 2002; Trampert et al. 2004; Wang & Wen 2004; To et al. 2005).

However, a large number of studies from geodynamics, seismology and mineral physics point towards a high core-mantle boundary (CMB) temperature and a correspondingly large temperature drop across D'', with associated high heat flow (e.g. Glatzmaier & Roberts 1995; Kuang & Bloxham 1997; Boehler 2000; Bunge et al. 2001; Steinle-Neumann et al. 2001; Alfè et al. 2002; Buffett 2002; Gubbins et al. 2004; Nimmo et al. 2004; Bunge 2005; Mittelstaedt & Tackley 2006; Nolet et al. 2006; Zhong 2006; Alfè et al. 2007; van der Hilst et al. 2007; Lay 2008; Steinberger & Holme 2008). As shown by Schuberth et al. (2009a,b), a large thermal gradient across D'', on the order of 1000–1500 K, may provide an alternative explanation for the seismic heterogeneity in the lowermost mantle. The strong lateral variations in shear wave velocity seen in tomographic models can be reconciled with isochemical whole mantle flow and a pyrolite composition in case of strong core heating. Most importantly, the resulting high plume excess temperatures in the lowermost mantle are capable of explaining the magnitude of the shear wave velocity reduction of as much as –4 per cent within the two large low-velocity provinces. Furthermore, also the observed strong gradients in shear wave velocity (i.e. sharp sides of the African superplume) can be reconciled with isochemical whole mantle circulation in the presence of a high core heat flux.

To warrant entirely consistent comparisons to tomographic models, Schuberth et al. (2009a) modified the mantle circulation models (MCMs) of Schuberth et al. (2009b) to reflect the limited resolution and long-wavelength nature of tomographic S-wave velocity models. They accounted for the effects of uneven data coverage and damping in tomographic inversions by ‘tomographically filtering’ the geodynamically predicted heterogeneity using the resolution operator of model S20RTS (Ritsema et al. 2004, 2007). Doing so allowed them to specifically consider the magnitudes of lateral temperature variations and corresponding seismic heterogeneity in the comparisons. This way, a quantitative assessment of the geodynamic models was possible based on a variety of statistical measures, such as the root-mean-square (rms) profiles, histograms and spectral power of heterogeneity.

In addition to matching the strength of heterogeneity in tomographic shear wave velocity models, Schaber et al. (2009) have shown that a large amount of core heating is also compatible with small rates of true polar wander as inferred from palaeomagnetic data. This strengthens the notion that seismic heterogeneity is likely dominated by thermal variations (e.g. Davies & Davies 2009; Dalton et al. 2009; Simmons et al. 2009), which may place limits on the possible role of chemical heterogeneity in the lower mantle.

As mentioned earlier, the interpretation of seismic structure in terms of chemical heterogeneity is to a large part based on differ-

ences in compressional ($v_p$) and shear wave velocity distributions imaged in the lower mantle. This is not only reflected in the anticorrelation of $\partial \ln v_s$ and $\partial \ln v_\phi$ inferred by tomographic inversion (Masters et al. 2000). In the long-period data set of Bolton & Masters (2001), which was used in this inversion (among other data sets), there are also obvious differences directly visible between the statistics of P- and S-wave traveltime variations.

It is clear that any geodynamical hypothesis and deduced models should stand the test not only against one specific type of data. The conclusions of Schuberth et al. (2009a,b) are mainly based on comparisons of MCMs to tomographic models of shear wave velocity. Only in fig. 9(d) of Schuberth et al. (2009b), they showed rms profiles of predicted P-wave heterogeneity in addition to their comparisons of shear wave velocity models, which indicated a similarly good agreement to tomography. However, repeating the exercise of tomographically filtering the velocity variations predicted from the MCMs in case of the P-wave heterogeneity was not possible at the time. Tomographic inversions for P-wave velocity structure typically involve a much larger number of free parameters compared to the long-wavelength S-wave models, so that the construction of the resolution operator is computationally very demanding. Alternatively, one could produce a synthetic set of traveltime variations by tracing rays through the geodynamically predicted seismic structures and invert this data set (e.g. Bunge & Davies 2001; Davies & Bunge 2001). This approach is equivalent to multiplying the model parameters with the resolution operator, but is computationally less convenient in case of long wavelength S-wave models.

In principle, it is thus possible to also mimic the effects of limited tomographic resolution for the predicted P-wave velocity variations. However, there are a number of issues that render such an analysis inconsistent (described below): (1) effects related to the reparametrization of the model in case of the resolution operator approach, (2) inconsistencies in interpreting both P- and S-wave data, which typically have different frequency content, in the context of ray theory and (3) the limited resolving power of tomographic inversions can be accounted for, but the non-uniqueness problem remains.

(1) In case of using the resolution operator, difficulties arise from the need to project geodynamic models onto the parametrization of the corresponding tomographic model. This reparametrization results in an unphysical loss of power of the seismic velocity perturbations, because of the fact that modern geodynamic models typically are generated on numerical grids that have a much larger number of points and higher resolution than is possible for tomographic models. Schuberth et al. (2009a) showed that it is in principle possible to correct for the effect of the reparametrization when looking at statistical measures such as rms profiles or histograms of velocity perturbations. However, the degree to which this effect should be corrected for depends on the amount of aliasing of short-scale structure into long-wavelength tomographic models, which is not a priori known (e.g. Mégnin et al. 1997; Boschi & Dziewonski 1999; Dahlen 2004). Furthermore, in case of long-period seismic data, wave front healing may play an important role by obscuring the information on short-scale structure in the lower mantle (Wielandt 1987; Marquering et al. 1999; Dahlen et al. 2000; Hung et al. 2000).

(2) Using ray tracing through the geodynamic model and subsequent inversion of predicted traveltimes is appropriate in case of short-period body wave data. For P waves, typically, the ISC (International Seismological Centre) data are used, which consist of a large number of first onset picks made at 1 s period, in which
case ray theory is valid (Nolet 2008). In case of S waves, however, picking the onset is much more difficult, as other phases may interfere with their arrival and they typically contain lower frequencies. Therefore, traveltime measurements are now often done using cross-correlation techniques, for which case ray theory is not appropriate (Nolet 2008).

(3) Any tomographic inversion suffers from the fact that an infinite number of models will fit the data equally well. This problem of non-uniqueness of the solution is inherent to tomographic inversions and is a consequence of the imperfect nature of seismic observations (i.e. errors in the data and inhomogeneous coverage). Tomographers have made great progress in finding new ways to efficiently explore the so-called null space (Deal & Nolet 1996; Trampert et al. 2004), but such approaches still rely on accurate prior information. To improve the resolution of tomographic models and to improve robustness of inversions, high expectations are currently placed on finite-frequency theory (also called banana-doughnut theory; e.g. Dahlen et al. 2000; Montelli et al. 2006; Sigloch et al. 2008) and numerical simulations of full 3-D wave propagation through complex 3-D structures. Inversions using the 'adjoint' method have lately been pioneered with the hope to exploit a larger part of the information contained in waveforms (Tromp et al. 2005; Fichtner et al. 2009; Tape et al. 2009; Fichtner et al. 2010).

In light of the above considerations, it seems desirable to eliminate the need to involve any inversion in the assessment of geodynamic models. In this study, we thus explore a new approach to test geodynamic hypotheses, which is based on forward modelling only: we compute synthetic seismograms using simulations of 3-D global wave propagation based on a spectral element method (SEM) for synthetic elastic structures predicted from an MCM.

Using state-of-the-art techniques to solve the wave equation in 3-D heterogeneous media, we capture the correct physics of wave propagation. Furthermore, we are not restricted to a single traveltime datum for each seismic phase, as in case of ray tracing, but we are able to measure traveltime variations at different frequencies. This way, we will be able to study frequency-dependent waveform effects, such as wave front healing and focusing/defocusing, due to 3-D heterogeneity in a physically consistent manner. A number of recent numerical studies of 3-D seismic wave propagation have, for example, pointed to an artificial bias in the interpretation of seismic observations due to an apparent anisotropy in isotropic media (e.g. Komatitsch et al. 2010; Lin & Ritzwoller 2011). These studies highlight the importance of taking 3-D waveform effects correctly into account. Moreover, seismic data sets are rapidly growing not only due to an increasing number of seismic stations, but also due to the fact that traveltime measurements are now starting to be done at multiple frequencies (e.g. Sigloch et al. 2008; Zaroli et al. 2010). It will thus become increasingly important to understand these multifrequency data sets from a forward modelling perspective complementary to using them for tomographic inversions.

The fact that this approach relies on forward modelling only has the advantage of avoiding all problems inherent to tomographic inversions, namely the limited resolving power and the non-uniqueness of the solutions. Only a set of parameters that enter the governing equations of mantle dynamics and information on mantle mineralogy are needed as input, resulting in synthetic seismograms that are predicted independently of seismic observations. Thus, the potential danger of circular reasoning is averted that may bias interpretations when using tomographic models for the simulation of wave propagation.

The general purpose of this study is to test whether the generation of synthetic wavefields for geodynamically derived seismic velocity models is suited to provide complementary information in a quantitative interpretation of seismic data. As a first application, we will use our approach to test the assumption of isochemical whole mantle flow with strong core heating directly against the statistics of seismic data.

Our approach goes somewhat along the lines of the 'COSY' project (Igel et al. 2000), which aimed at validating numerical waveform modelling techniques using an MCM. In contrast, our goal is to use an MCM in simulations of wave propagation to test geodynamic hypotheses directly against seismic observations. Another difference is that at the time of the COSY project, computing resources did only allow to use axisymmetric methods for 2-D sections, whereas nowadays 3-D simulations in a full sphere can be done on a routine basis on modern large-scale supercomputers and are possible even on intermediate sized PC-clusters (e.g. Oeser et al. 2006). A similar philosophy to our forward modelling approach was pursued recently by Jacobs & van den Berg (2011), who, however, did not compute the full wavefield but rather reflectivity profiles for the transition zone discontinuities. Furthermore, they restricted themselves to 2-D geodynamic models in combination with a thermodynamic model of mantle mineralogy. Our approach goes a long way further in that it is based on 3-D high-resolution MCMs and results in full seismic waveforms taking into account all possible waveform effects.

We start this paper with a detailed description of our methodology, including a short review of the mantle circulation model of Schuberth et al. (2009b), which we use for the simulations of seismic wave propagation. Particular issues related to the implementation of a geodynamic model into the software for seismic wave propagation are discussed in the appendix. Using a very large number of evenly distributed 'virtual' seismic stations, we obtain a data coverage suitable for comparison to global seismic observations even when simulating only 17 earthquakes. Owing to the large number of resulting synthetic seismograms, we use automated software for the measurement of traveltime variations, which we will describe briefly.

In this study, we will focus on *P*- and *S*-wave traveltime variations in a single frequency band and compare their standard deviation to that of observed traveltime variations. The analysis of traveltime variations themselves in different frequency bands and a discussion on frequency-dependent waveform effects will be postponed to a future publication. For the comparison to our predicted traveltime variations, we use the data set of Bolton & Masters (2001), which consists of ~40 000 delay times of long-period *P* and *S* waves measured in the same frequency band. This data set was chosen, because the simulation of 3-D global wave propagation for the respective frequencies is feasible for a reasonable number of earthquakes. More important, it is this data set that has been used to infer anti-correlation of bulk sound speed and shear wave velocity in the two large superplumes in the lowermost mantle (Masters et al. 2000), as mentioned earlier. Furthermore, besides the fact that the ratio $\partial \ln v_s / \partial \ln v_p$ derived directly from these data increases strongly with depth, Bolton & Masters (2001) made the intriguing observation that the standard deviation of delay times increases strongly as a function of ray-theoretical turning depth in case of *S* waves, while being almost constant throughout the lower mantle for *P* waves. Thus, the comparison of the statistics of our synthetic traveltime variations to that of these long-period observations presented here provides a further test of the hypothesis of a large core contribution to the mantle energy budget.
## 2 METHODOLOGY

As outlined in Section 1, we present a multidisciplinary approach that combines forward modelling techniques from geodynamics, mineral physics and seismology. 3-D high-resolution MCMs are used together with thermodynamic models of mantle mineralogy to predict the seismic heterogeneity of the mantle. The predicted seismic models are implemented in a spectral element code for the simulation of 3-D global wave propagation. In the following, we will discuss our modelling setup in detail.

### 2.1 Geodynamic input model

We use the mantle circulation model (MCM) M2 of Schuberth et al. (2009a,b), which compared best to tomographic shear velocity models in terms of spectral characteristics and, most important, the magnitude of velocity anomalies. Mantle general circulation was modelled with the parallel finite element code TERRA (Bunge & Baumgardner 1995; Bunge & Richards 1996; Bunge et al. 1996, 1997). A key improvement to earlier studies (e.g. Bunge et al. 2002) was the very high resolution of the mesh with more than 80 million finite elements. Sufficiently high numerical resolution is crucial to model global mantle flow at Earth-like convective vigour, which is one of the reasons for the good match to the strength of heterogeneity observed in tomographic models. Apart from the high resolution, Schuberth et al. (2009b) restricted themselves to a minimum number of assumptions: (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 postglacial rebound and geoid observations (e.g. Hager & Richards 1989; Paulson et al. 2007) and (3) isochemical whole mantle flow with a pyrolite composition.

In model M2, a very high core heat flow of 12 TW (~35 per cent of the outflow at the top of the mantle) was obtained by setting the temperature at the CMB to 4200 K. Viscosities of $10^{23}$, $10^{21}$ and $10^{23}$ Pa s were used for the lithosphere, the upper and the lower mantle, respectively, which are separated at depths of 100 and 650 km. Seismic heterogeneity was then predicted for a pyrolite composition by post-processing the temperature field of the isochemical MCM with the thermodynamic models of mantle mineralogy of Piazzoni et al. (2007) and of Stixrud & Lithgow-Bertelloni (2005b, 2007). In both models, equilibrium phase assemblages in the CFMAS (CaO – FeO – MgO – $Al_2O_3$ – $SiO_2$) system are computed by Gibbs free energy minimization. Here, we adopt the model of Stixrud & Lithgow-Bertelloni (2005b, 2007) again, but with the difference that we now use their seismic velocities corrected for the effect of anelasticity, which has been shown to be non-negligible (e.g. Karato 1993; Goes et al. 2004; Brodholt et al. 2007; Matas et al. 2007). The 1-D Q profile of PREM was used for the anelasticity correction, which, in the lower mantle, leads larger sensitivities to temperature for both $v_p$ and $v_s$ (on the order of 30 per cent). The resulting model of variations in $v_p$, $v_s$ and density, which we use throughout all the wave propagation simulations, will be denoted S09-M2-Q in the following.

As mentioned earlier, model S09-M2-Q shows realistic magnitudes of seismic velocity perturbations and explains the strength of heterogeneity in tomographic models well. It is limited, however, in that it is not able to correctly predict the exact location and pattern of structure in the lowermost mantle due to uncertainties in the initial and boundary conditions used in the simulation of mantle flow. The model exhibits obvious differences to tomographic models in terms

of the geographic distribution of fast and slow anomalies close to the CMB, as discussed in Schuberth et al. (2009b). For fast material, this is mainly due to increasing uncertainties in plate motion history models for earlier times and the limited time span that they typically cover. For example, the plate configuration and directions of motions given by the earliest plate stage of the reconstructions of Lithgow-Bertelloni & Richards (1998), which were used to construct model S09-M2-Q, result in fast material in the lowermost mantle under eastern Africa. More recent models of plate configurations dating further back in time and the use of better absolute plate motion reference frames provide a better geographic fit of MCMs to tomographic images also in the lower mantle (Shephard et al. 2012), which will help to reduce the impact of such limitations in future.

The location of plumes in the lower mantle poses a more fundamental problem, as models of plate motion history provide only constraints for the location of downwelling slab material, but there are no equivalent constraints for hot upwelling material at the CMB. In model S09-M2-Q, for example, a big upwelling exists in the southeast Pacific, just under Easter Island, a region of long-lived divergence of plates, resulting in upward flow in our model at this location. Especially for hot upwelling material, inverse modelling of mantle dynamics with the goal to constrain the initial condition will hopefully provide a way to improve the location of predicted structure in the lowermost mantle (e.g. Bunge et al. 2003).

A further limitation of model S09-M2-Q is that the absolute values of seismic velocities resulting from the conversion of absolute temperatures using the mineralogical model are different from those observed on Earth by up to 2 percent. This is a consequence of uncertainties in experimental mineralogical data (e.g. Kennett & Jackson 2009) and the unknown exact shape of the geotherm. This will result in rather different absolute traveltimes of seismic waves in our model compared to Earth. This, however, does not inhibit a comparison to observed data: understanding the lateral variations in seismic velocities in the mantle is all about analysing differences in the arrival time of seismic waves with respect to predictions from a 1-D reference model. In other words, we are currently not interested in explaining the absolute values of seismic velocities, but rather concentrate on their lateral variations as predicted by our geodynamic simulation. A number of studies have focused on the radial average properties of the mantle in the last years, but so far, interpretations have remained limited as the trade-offs between thermal and chemical effects have not allowed an unequivocal identification of the cause of seismic structure both for the upper and the lower mantle (Cammarano et al. 2003; Deschamps & Trampert 2004; Cammarano et al. 2005; Mattern et al. 2005; Matas et al. 2007; Cobden et al. 2008; Khan et al. 2008; Cobden et al. 2009; Styles et al. 2011). One choice of reference for our simulations could be one of the seismological reference models, which we could use together with the demeaned velocities predicted by the MCM. However, as will be discussed in Section 2.2 and the appendix, a better choice of reference are the average profiles of $v_p$, $v_s$ and density of the geodynamic model itself resulting in self-consistent velocity perturbations at the phase transitions.

As can be seen from the above considerations, exact comparisons of traveltime variations or full waveforms predicted from our model with real seismograms are out of reach at the current state of knowledge on mantle dynamics and mineral physics. In this study, we will try to show, however, that a comparison of the statistics of both the synthetic and observed data set is possible.
## 2.2 Simulation of 3-D global wave propagation

For solving the wave equation in (an-)elastic heterogeneous media, a variety of techniques is available nowadays. One prominent example is the SEM, which has proven to be well suited for 3-D wave propagation simulations due to its high accuracy, efficient parallelization, implicitly incorporated free-surface boundary conditions as well as its geometrical flexibility (e.g. Priolo & Seriani 1991; Seriani et al. 1992; Seriani & Priolo 1991, 1994; Seriani 1997; Komatitsch et al. 1999; Komatitsch & Tromp 1999; Komatitsch et al. 2000; Komatitsch & Tromp 2002a,b; Capdeville et al. 2003; Chaljub & Valette 2004; Komatitsch et al. 2004; Peter et al. 2011). One implementation of the SEM called 'SPECFEM3D' (available at: www.geodynamics.org) became a standard tool in computational seismology in recent years, which offers the possibility to incorporate a large range of physical phenomena in addition to the possibility to use 3-D heterogeneous velocity distributions (e.g. effects of topography/bathymetry, attenuation, anisotropy, finite sources, etc.). In this study, we use 'SPECFEM3D_GLOBE' (V5.0.1), the global version of this software.

In this study, we wish to analyse the direct body wave phases at frequencies of 0.1 Hz and below. The size of the numerical grid is thus chosen such as to generate synthetic seismograms with sufficient numerical accuracy down to a shortest period of 10 s (NEX = 432 resulting in ~1.3 billion gridpoints), and we use the 3-D distribution of seismic velocities predicted by model S09-M2-Q projected onto this grid. The implementation of our geodynamic model, in particular the issues related to the projection of such high-resolution models onto the even larger 'SPECFEM3D_GLOBE' grid are described in the appendix. There, we also discuss the advantages of using the temperature field for the projection (instead of seismic velocities computed on the grid of the convection simulation) and doing the conversion to seismic velocities afterwards on the 'SPECFEM3D_GLOBE' grid.

To obtain a reference seismogram for the cross-correlation measurements described in Section 2.4, we repeat all simulations with the average 1-D structure of S09-M2-Q. This 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 cancel out when later doing 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).

As we are interested in the seismic signature of a geodynamic hypothesis for the mantle, we try to isolate those wavefield effects that result from 3-D velocity variations in the mantle alone. Therefore, we assign a 1-D velocity model not only for the core, but also for the crust, for which we use the 'continental' version of 'PREM' (i.e. without the ocean layer). We further reduce the complexity in the wave propagation simulations to its minimum by switching off all other options available in 'SPECFEM3D_GLOBE' to compute realistic seismograms; that is, we do not include the effects of attenuation (dissipation of seismic energy due to anelasticity), anisotropy, topography/bathymetry, ellipticity, etc. At the current state of knowledge on mantle dynamics and mineral physics (cf. Section 2.1), we are not primarily interested in generating realistic seismograms that show a complexity comparable to real observations.

On the contrary, adding complexity to our simulations would require unnecessary additional post-processing of the synthetic seismograms, just the same as is needed in case of real data to extract the information on 3-D mantle heterogeneity (e.g. crustal corrections, estimation of $t^*$ and correction for dispersion related to intrinsic attenuation, estimation of the source time function, etc.). Typically, most of these post-processing techniques rely on the assumption of linearity of these effects anyway.

In our case, where we wish to study the wavefield effects resulting from variations in the isotropic elastic structure of the mantle, we thus opt to neglect the aforementioned phenomena in the first place. We treat the simulations for the 1-D model in exactly the same way as the 3-D simulations. Thus, including, for example, a 1-D Q structure or adding a 3-D crustal model to both the 1-D and the 3-D simulations should give the same results. We tested this for one earthquake adding the 3-D crustal model CRUST2.0 (Bassin et al. 2000) to both the 1-D and the 3-D simulations and found no significant difference in the resulting traveltime variations compared to those obtained without a 3-D crust.

The same will probably be true when incorporating attenuation with a 1-D Q structure, though neglecting attenuation here while taking the anelasticity correction into account in the temperature to seismic velocity conversion (see Section 2.1) is not entirely consistent. However, we wish to use the simulations also for future analyses, for example, to study the traveltime dispersion resulting from purely elastic structures. In addition, the effect of attenuation will probably change the 1-D and 3-D waveforms in almost the same way (i.e. second-order effects will be small), so that differences in the resulting traveltime variations are also expected to be insignificant. It remains to be seen whether this will change when adding 3-D variations in Q. The work of Tian et al. (2011) indicates that 3-D Q structure has a very minor effect on S-wave delay times compared to 1-D Q, and Sigloch et al. (2008) have already shown earlier that one can neglect 3-D Q completely for P waves. A similar result was obtained by Savage et al. (2010), who observed only a minor effect of 3-D Q on the traveltimes of seismic waves.

## 2.3 Earthquakes and seismic stations used for the simulations

Using the setup described in Section 2.2, we simulated seismic wave propagation on the global scale for 17 different earthquakes (see Table 1). The location and moment tensors of these events are plotted in Fig. 1 on top of the shear wave velocity perturbations of model S09-M2-Q at a depth of 50 km. The events consist of two subsets: (1) 12 events were selected from the Harvard CMT catalogue (www.globalcmt.org) such as to fall closest to the corners of an icosahedron with one of the corners located in Iceland. This selection of earthquakes aims at providing an almost even sampling of the Earth with seismic waves, (2) five additional events were selected from the catalogue to illuminate the prominent plume structure present in model S09-M2-Q in the southeast Pacific.

Calculating seismograms for this number of seismic sources and using only the existing seismic stations as receivers would result in a poor illumination of the geodynamically predicted seismic structure. Therefore, we additionally 'record' the seismic wavefield with a network of 42 250 equidistant 'virtual' seismic stations (~100 km apart), thereby drastically increasing the number of synthetic seismograms.

Figs 2 and 3 show examples of synthetic seismograms for our synthetic mantle structure. A movie of the wave propagation in the geodynamic model can be found in the Supporting Information
<table><caption>**Table 1.** Earthquakes simulated in this study. Names and location, etc., are taken from the Harvard global CMT catalogue (www.globalcmt.org).</caption><thead><tr><th>Nr.</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>Southwest 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>00909122006A</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_4.png" /><figcaption>**Figure 1.** Locations and Harvard moment tensor solutions (www.globalcmt.org) of the 17 earthquakes used in this study. The events are plotted on top of the shear wave velocity perturbations from model S09-M2-Q at a depth of 50 km. Moment magnitudes are in the range from 5.3 to 7.0 (see also Table 1).</figcaption></figure>

section online. Fig. 2 shows record sections of vertical and transverse component seismograms for a shallow event that happened at the central Mid-Atlantic Ridge ($M_w$ 6.2; event 8 in Table 1). To our knowledge, these are the first seismograms computed for a global 3-D MCM with full 3-D wave propagation techniques and accurate down to a shortest period of 10 s. On the vertical component seismograms, one can clearly see the direct *P* wave as well as its multiple surface reflections. Beyond 140°, the first visible phase is *PKP*. Between 60° and 70° epicentral distance, amplitudes of the *P* waves are larger than at shorter distance, which indicates focusing due to 3-D structure in our geodynamic model. The transverse component seismograms in Fig. 2, being mostly clean of *P*-wave energy, clearly show the direct *S* and the surface-reflected *SS* arrivals. Between 50° and 90° epicentral distance, the core reflected *ScS* phase can be distinguished. Energy visible after *SS* at distances beyond 120° is probably related to the arrival of the *ScSScS* wave. Owing to the fact that we neglect the attenuation of seismic waves, amplitudes decrease less rapidly with epicentral distance as one may anticipate from seismic observations.

Fig. 3 shows the three components of a seismic ‘recording’ at station PFO (Piñon Flat Observatory) for a deep earthquake

under western Brazil ($M_w$ 7.0; event 7 in Table 1). A number of secondary arrivals were identified in addition to the main seismic phases based on ray-theoretical arrival times [computed with the ‘TauP Toolkit’ (Crotwell *et al.* 1999) for the 1-D average velocity and density profiles of model S09-M2-Q]. Owing to the great depth of the earthquake (556 km), the depth phases are clearly separated from the main arrivals. Some of the arrivals show rather large amplitudes as a consequence of neglecting the intrinsic attenuation of the medium. Especially, the amplitude ratio between *S* and *P* waves is too large.

## 2.4 Automated measurement of traveltime variations

From our simulations, using 17 earthquakes and 42 250 stations, we obtain a large number of seismograms, for which we want to measure traveltime variations by cross-correlation. Obviously, an interactive measurement procedure as the one of Bolton & Masters (2001) is unfeasible in view of the more than 1.5 million traces that we have to process. Instead, we use the software of Zaroli *et al.* (2010), which performs an automated window selection and subsequent measurement of frequency-dependent traveltimes based on cross-correlation of the trace of interest with a reference seismogram. This software 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 just separating the seismograms into portions containing coherent seismic energy and portions containing mostly noise. Initially designed for measurements on seismic observations, this software also performs well in the absence of noise, as is the case for synthetic seismograms.

We applied the automatic multifrequency measurement technique of Zaroli *et al.* (2010) to a selection of the vertical (for *P*-wave measurements) and transverse component (for *S*-wave measurements) seismograms simulated in this study. A minimum epicentral distance 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. Only those measurements are retained, for which the cross-correlation coefficient was greater than 0.8. This removes about 0.2 per cent of the *P* wave and around 2 per cent of the *S*-wave measurements, which indicates that perturbations to the wavefield due to the 3-D heterogeneity in the MCM are small, especially in the case of *P*-wave propagation. This can be interpreted with a prevalence of single scattering events and only
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Epicentral distance plots for an earthquake at the central Mid-Atlantic Ridge ($M_w$ 6.2; event 8 in Table 1). The synthetic seismograms depict ground velocity and were computed using SPECFEM3D_GLOBE (Komatitsch & Tromp 2002a,b) for the 3-D mantle circulation model S09-M2-Q. The traces have been low-pass filtered at 0.09 Hz to remove numerical noise below 10 s period. Surface wave energy has been removed for better visualization of body wave phases and traces have been shifted along the time axis with a dynamic delay given by 0.04Δ s. The dissipation of seismic energy due to intrinsic anelasticity has not been taken into account in the simulation.</figcaption></figure>

a small contribution from multiple scattering. As our hypothetical mantle structure shows realistic spectral characteristics, a minor importance of multiple scattering of long-period waves may also be anticipated for Earth's mantle. This actually indicates that single scattering Born theory is appropriate for the interpretation of seismic observations.

Traveltime delays were measured in four different frequency bands using overlapping gaussian filters (15, 22, 34 and 51 s dominant period). As will be discussed in Section 2.5, we concentrate on the traveltime measurements at 15 s in this study. In total, we

obtained 363 893 measurements at 15 s for *P* and 350 387 measurements for *S*. It is worth noting that the broad-band (11–85 s) traveltime variations in our data set are very similar to those at 15 s (especially their statistics).

## 2.5 Cross-correlation at 15 s versus ‘first upswing’ measurements at 25 s

As outlined in Section 1, we wish to compare our synthetic data to the observations of Bolton & Masters (2001). Their measurements
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Example of a synthetic seismogram computed with the geodynamic model S09-M2-Q for a deep earthquake under western Brazil (*M*<sub>w</sub> 7.0; event 7 in Table 1). Shown are the three components of ground velocity ‘recorded’ at station PFO (Piñon Flat Observatory), which have been low-pass filtered at 0.1 Hz. Intrinsic attenuation was not taken into account in the simulation resulting in rather strong signals for some of the seismic phases. Red arrows and labels indicate the main seismic phases, as well as some of the secondary arrivals that we identified based on ray-theoretical arrival times. Owing to the great depth of the earthquake, the depth phases are clearly separated from the main arrivals. Note the occurrence of *P*-wave energy on the transverse component at times prior to the first *S*-wave arrival, which indicates *P* to *SH* scattering due to 3-D heterogeneity.</figcaption></figure>

were done on *P*- and *S*-waveforms of same frequency content (25 s dominant period) using a ‘first upswing’ philosophy; that is, only the first quarter of the waveform is considered and aligned with the same portion of a reference waveform. Introduced by Masters *et al.* (1996), this procedure is supposed to not only avoid problems related to depth phases or interfering core reflections, but to also highlight the higher frequency portion of the waveforms, thereby reducing the influence of wave front healing on the measurements. However, finite-frequency effects are to be expected just as with cross-correlation traveltimes, since the upswing delay time is influenced by scattered or diffracted energy that arrives after the minimum traveltime. Although no formal Fréchet derivative (i.e. the sensitivity of the measurements to changes in seismic velocities) exists for upswing delays, these estimates are expected to be very similar to cross-correlation measurements, owing to their similar nature.

This notion was confirmed by Hung *et al.* (2001), who compared ‘ground truth’ cross-correlation measurements (made on synthetics obtained with a pseudospectral method for a simple velocity perturbation included in an otherwise homogeneous medium), with predictions from banana-doughnut theory as well as with measurements based on the ‘first upswing’ philosophy. They found that traveltime predictions from banana-doughnut kernels with an effective period (*T*<sub>eff</sub>) of 15 s correspond best to the ‘first upswing’ measurements done at 25 s period. This work showed that the ‘first upswing’ measurements indeed highlight the higher frequency portion of the waveforms and guided our choice of using cross-correlation measurements at 15 s for the comparison of synthetic traveltime variations to the data of Bolton & Masters (2001). However, fig. 23 of Hung *et al.* (2001), displaying the effect of wave front healing as a

function of wavelengths propagated behind an anomaly, shows that the predictions of banana-doughnut theory do slightly differ from the ‘ground truth’ cross-correlation measurements, the latter actually suffering more from wave front healing at long propagation distances. The same figure also shows that the banana-doughnut kernels with $T_{\text{eff}}$ = 15 s predict slightly larger wave front healing for fast anomalies than is captured by the ‘first upswing’ measurements at long propagation distances, while predicting less at shorter distances. A similar bias was not observed in case of slow anomalies. Taken together, this suggests that the full waveform cross-correlation measurements at 15 s may be more affected by wave front healing compared to the measurements by Bolton & Masters (2001). In particular, this is probably more severe in case of *S* waves, for which the number of wavelengths propagated is larger than for *P* waves due to their smaller wavelength.

The exact difference between banana-doughnut predictions, ‘first upswing measurements’ and the ‘ground truth’ cross-correlation traveltime variations depends on the ratio of the size of the anomaly to the wavelength of the wave (Hung *et al.* 2001; Malcolm & Trampert 2011). As *P* and *S* waves have different wavelength but propagate through the same structure, it is not straightforward to estimate the difference in traveltimes and resulting standard deviation arising from the different measurement philosophies. As Hung *et al.* (2001) analysed the most simple case of a single heterogeneity inside a homogeneous medium only, we additionally tested whether in our more complex case, the different measurement procedures give results similar enough for our purposes. For this, we modified slightly the automated cross-correlation software described in Section 2.4 to select time windows containing only the first portion of the waveform and to measure the delay time at a dominant period
of 25 s. The traveltime difference with respect to the 1-D reference trace was then taken to be the time shift that minimizes the rms of the difference between the 3-D and 1-D traces. This automated procedure resembles the interactive ‘first upswing’ measurements of Bolton & Masters (2001) very well, and we applied it to our full set of synthetic seismograms. Comparison to the measurements done with the original version of the software (i.e. cross-correlation of full waveforms with 15 s dominant period) confirmed that the two measurement techniques give very similar results. Most important for the comparison to the observations discussed later, this is in particular true for the corresponding standard deviations of the delay times (see Section 3.2.2 and Fig. 7). However, in line with the considerations on wave front healing effects discussed earlier, we find that the ‘first upswing’ type of measurements consistently show a slightly larger standard deviation, especially in the case of S waves.

The similarity between the two measurement procedures means that we can use either of the two for our synthetic data set. We will focus our attention on the full waveform cross-correlation traveltime variations at 15 s period, as they have the advantage that their Fréchet derivative is well-established in the context of finite-frequency theory, and the corresponding banana-doughnut kernels allow for a direct and intuitive interpretation of the results. In addition, more and more seismic data sets exist nowadays, for which traveltimes have been measured in several different frequency bands using automatic algorithms based on cross-correlation of waveforms similar to the one used here (e.g. Sigloch et al. 2008; Zaroli et al. 2010). Thus, we will be able to also compare our measurements on the synthetic seismograms to such automatically measured observed data once a global combined P- and S-wave data sets of cross-correlation traveltime delays is available.

# 3 RESULTS

## 3.1 Maps of P- and S-wave traveltime variations

Before we concentrate on the statistics of our synthetic data, it is informative to look at the geographical distribution of traveltime variations computed for the synthetic mantle structure and to check whether the traveltimes carry any imprint of the MCM. Figs 4 and 5 show maps of P- and S-wave traveltime variations plotted at their corresponding station locations for four different earthquakes: Fig. 4 shows the results for events in Iceland (event 5 in Table 1) and on the Southwest Indian Ridge (event 14), whereas Fig. 5 shows events at the central Mid-Atlantic Ridge (event 8) and offshore the northeast coast of Russia (event 3). For the first three earthquakes, more positive than negative anomalies are visible, which is a direct consequence of the events being located in regions of slow seismic velocities (i.e. mid-ocean ridges).

Most importantly, the pattern of traveltime variations corresponds to the near-surface structure of model S09-M2-Q in all cases. In case of the positive traveltime anomalies, some spreading centres in the oceans can be distinguished: for the Iceland earthquake, the Mid-Atlantic Ridge is visible southwest of the event as well as the northern East-Pacific Rise southwest of the North American continent; for the Southwest Indian Ridge earthquake, the Mid-Atlantic Ridge is (here only weakly) visible between South America and Africa as well as the Central and Southeast Indian Ridge and the Pacific-Antarctic Ridge; and for the central Mid-Atlantic Ridge earthquake, the northern part of the ridge can be followed running northwards from the event towards the Arctic Sea. The southern

Mid-Atlantic Ridge is also visible as well as the East Pacific Rise, which runs from Mexico southwards until it hits the strongest positive anomaly located at the triple junction of the Pacific, Nazca and Antarctic plates. This large region of positive traveltime variations results from the slow seismic velocities of the big prominent plume structure in model S09-M2-Q mentioned in Section 2.1.

The largest negative traveltime anomalies in all maps of Figs 4 and 5 are clearly related to regions of subduction: for the Iceland event, the subduction of the Farallon Plate under North America is visible as three different slab segments. The subduction of the Pacific Plate under Eurasia is apparent from negative anomalies at the Aleutians, Japan and East Asia, and the Tethys subduction is visible from the Mediterranean to the Himalayas; for the Southwest Indian Ridge event, one can also see the imprint of the Tethys slab and the subduction along the Sumatra arc, as well as the southernmost end of the slab under South America; for the central Mid-Atlantic Ridge event, the whole length of the American continents shows negative traveltime anomalies related to subduction, as is the case for the Tethys subduction under the Mediterranean.

The fourth earthquake, the event offshore the northeast coast of Russia (bottom row in Fig. 5), is different from the others in that it is a deep subduction earthquake located in a region of fast seismic velocities. This leads to generally more negative traveltime variations compared to the other events. For example, the whole Eurasian continent and most of North America show negative anomalies compared to mostly positive anomalies in case of the Iceland or central Mid-Atlantic Ridge earthquakes. Still, the near-surface structure of model S09-M2-Q with fast seismic velocities is apparent as negative traveltime anomalies (Marianna trench, Fiji-Tonga trench, Sumatra arc and subduction under Alaska). However, owing to the great depth of this event, there are also regions in the oceans, which exhibit negative traveltime variations. This indicates a stronger influence of deeper structure on the traveltimes in this case.

For the event offshore the northeast coast of Russia, a ring-shaped pattern in the P-wave anomalies can be seen in the North Pacific, which seems to coincide with a nodal line in the radiation pattern of the earthquake source (cf. CMT solution in bottom row of Fig. 5). We attribute this pattern in the traveltimes to diffraction of seismic energy around heterogeneities in the model (in this case, probably the slab close to the source). Strong signals in the traveltime variations happen to appear along the nodal planes when energy from 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 that travels in the corresponding 1-D structure, which gives rise to large delay times when measured by cross-correlation. A similarly remarkable feature exists in the map for S-wave traveltimes: a linear band of strongly positive anomalies that runs through the Pacific from the event southwards passing by Hawaii. This again is interpreted as diffraction of energy around 3-D structure together with the different radiation pattern for S-wave energy. However, a more rigorous analysis of diffraction effects in 3-D geodynamic structures is needed to confirm our interpretation, which is beyond the scope of this study.

Both, the strong influence of the near-source and the near-receiver structure can be understood in the context of finite-frequency theory, together with the heterogeneity distribution in S09-M2-Q. The Fréchet kernels for traveltime anomalies typically show a very high sensitivity close to the endpoints of the propagation path, where their cross-sectional area becomes very small (e.g. Dahlen et al. 2000; Hung et al. 2001; Nolet 2008). In addition, seismic heterogeneity is strongest close to the surface (cf. fig. 5f of Schuberth et al. 2009b), where also the size of the anomalies tends to be smaller than
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Maps of traveltime variations in model S09-M2-Q, measured by cross-correlation of full waveform synthetic seismograms. Left-hand column: Traveltime variations of direct *P* waves for earthquakes in Iceland (top row; event 5 in Table 1) and at the Southwest Indian Ridge (bottom row; event 14 in Table 1). Right-hand column: Same for *S* waves. The traveltime anomalies are plotted at the location of their respective receiver. A minimum epicentral distance range of 30° is used to guarantee a clear separation of the direct phases from later arrivals and to avoid problems due to upper-mantle triplications. Note the different colour scales for *P* and *S* waves.</figcaption></figure>

in the lower mantle. The combination of both strong short-scale heterogeneity and a high sensitivity of the waves close to the surface causes the lithospheric structure to dominate the traveltime signals.

From Figs 4 and 5, the direct correspondence between $\partial \ln v_p$ and $\partial \ln v_s$ resulting from the isochemical nature of the geodynamic model is apparent from the strong resemblance of *P*-wave and *S*-wave traveltime patterns. Regions of negative and positive traveltime anomalies for *S* waves are almost in the same locations as for the *P* waves, but are larger by a factor of at least 2 (cf. colour scales

in Fig. 4). However, for the case of a deep earthquake (second event in Fig. 5), this is not strictly true anymore. Although the regions of fast seismic velocities related to subduction in the model appear as negative traveltime anomalies in both *P* and *S* (e.g. Marianna trench, Fiji-Tonga, Sumatra trench, Himalayas), there are some regions in the Indian Ocean, where the pattern of arrival times is different between the two types of waves. This can be taken as a first indication that *P* and *S* waves at the same frequency do actually not 'see' the exact same structure due to the size of their Fresnel zone (i.e. the
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Same as Fig. 4 for events at the central Mid-Atlantic Ridge (top row; event 8 in Table 1) and near the northeast coast of Russia (bottom row; event 3 in Table 1). The latter event is a deep earthquake (520 km depth) located in a region of fast seismic velocities, which results in more negative traveltime anomalies compared to the Mid-Atlantic Ridge event (e.g. Eurasia shows slightly negative anomalies in the bottom plots compared to positive anomalies in the maps at the top). The ring-shape pattern in *P*-wave anomalies visible in the North Pacific is attributed to diffraction of the waves around the slab close to the source.</figcaption></figure>

width of their banana-doughnut kernel) being different. The cross-sectional area of the *P*-wave kernels is almost twice as large as for *S* waves as a consequence of the higher *P*-wave velocities [the $v_p/v_s$ ratio is typically on the order of 1.8 in the mantle (Kennett et al. 1995)], and the width of the kernels is proportional to $\sqrt{\lambda L}$ (Nolet 2008). In other words, wave front healing affects *P* waves more than *S* waves of the same frequency [see also Malcolm & Trampert (2011) for a nice numerical study of waveform effects and related traveltime measurements]. This is probably also the reason why the

strong positive anomaly in the southeast Pacific seen for the central Mid-Atlantic ridge event (top row in Fig. 5) is larger and shows more detail in case of the *S* waves compared to the *P* waves.

## 3.2 Statistics of *P*- and *S*-wave traveltime variations

In the following, we will analyse the statistics of traveltime variations in the MCM and concentrate here especially on histograms and the variation of their width (i.e. standard deviation) with depth.
This is done to better understand how the seismic velocity pertur-
bations in model S09-M2-Q affect the distribution of traveltime
anomalies and to check whether the magnitude of traveltime varia-
tions in the geodynamic model are compatible with the statistics of
the observations.

3.2.1 Histograms of synthetic traveltime variations

Fig. 6 shows histograms of traveltime variations for our complete synthetic data set. The plots show the number of measurements binned as a function of delay time and turning depth on a logarithmic colour scale. The turning depth of the waves at 15 s period is taken to be the same as the ray-theoretical turning depth following seismological standard practice. Maximum delay times displayed in the histograms are ±7.5 s for *P* and ±15 s for *S* waves, which are similar to typical values chosen for seismic observations (Bolton & Masters 2001; Houser *et al.* 2008; Zaroli *et al.* 2010). On the left of each histogram, we show the corresponding values of the ‘scaled median average deviation’ (SMAD) of the histograms, which is a

<figure><img src="image_9.png" /><figcaption>Figure 6. Histograms of *P*- and *S*-wave traveltime variations as a function of ray-theoretical turning depth. The colour-scale gives the number of counts per bin on a logarithmic scale. The plots on the right-hand side show the scaled median average deviation (SMAD), which is a robust estimator for the standard deviation of the histograms (Kleiner & Graedel 1980). The grey solid and dashed lines indicate the location of the mean and maximum values of the histograms. Note the different scales of the *x*-axes for *P*- and *S*-wave traveltime variations as well as the different trends with depth of the SMAD curves.</figcaption></figure>

robust estimator for their standard deviation (Kleiner & Graedel
1980). We use the SMAD as a measure for the standard devia-
tion here, as it was also used by Bolton & Masters (2001) for the
observed traveltime variations.

The histogram of *P*-wave traveltime measurements shows a re-
markable trend of the maximum delay times with depth: in the
upper part of the lower mantle, the histograms have negative skew-
ness (i.e. larger negative variations than positive ones), whereas this
changes gradually to a positive skewness towards the CMB. This
trend is similar in case of the *S* wave but less apparent, as we chose
to limit the range of traveltime variations in the histogram plot to
±15 s, and some measurements actually fall beyond this limit.
From the SMAD curves on the right, however, it becomes evident
that there is a general difference between *P* and *S* waves in terms of
the spread of traveltime variations, which will be discussed in detail
in Section 3.2.2.

Owing to the logarithmic colour scale, one can clearly see from
the band of dark colours that most of the measurements fall within
the range of ±2.5 s for *P* waves and ±5 s for *S* waves. This can
readily be understood from the histograms of seismic velocity
perturbations in model S09-M2-Q, which show that most of the
mantle has rather small velocity variations (cf. fig. 6f of Schuberth
et al. 2009b); that is, the mantle is well mixed due to vigorous con-
vection at high Rayleigh numbers representative of the convective
regime of the Earth. A similar concentration of traveltime variations
around zero and a sudden drop of the number of measurements
towards larger values can also be seen in the observed data (cf.
fig. 7 of Bolton & Masters 2001). The values at which this drop
occurs are slightly larger than in case of our MCM. Note, however,
that our synthetic data set is free of any errors (such as earthquake
location errors, random errors, etc.), which contribute to the vari-
ance of the observed traveltime variations. Furthermore, no crustal
corrections have been applied in the synthetic case, which further
increase the variance in case of the observations (Bolton & Masters
2001). Therefore, we will perform a more detailed comparison of
the standard deviations of synthetic and observed *P*- and *S*-wave
traveltime variations in the following section. For this comparison,
we try to extract the signal from the SMAD curves of Bolton &
Masters (2001) that is solely due to the 3-D mantle structure.

3.2.2 Standard deviation of synthetic and observed traveltime variations

In the previous section, we have already introduced the SMAD,
which Bolton & Masters (2001) used to estimate the standard devi-
ation of observed *P*- and *S*-wave traveltime variations (SMAD and
standard deviation σ will be used interchangeably in the following).
Now we wish to analyse the standard deviation of traveltime varia-
tions in model S09-M2-Q in more detail and to compare it to that
part of the observed standard deviations only that can be attributed
to seismic heterogeneity in the mantle. As discussed in Section 2.5,
we concentrate on the synthetic traveltime measurements that were
done using cross-correlation of full waveforms at 15 s period in the
following, as the corresponding banana-doughnut kernels allow for
a direct interpretation of the results. Fig. 7 shows the SMAD curves
of the *P*- and *S*-wave histograms of Fig. 6 (i.e. geodynamically pre-
dicted) together with a range of values estimated from the data of
Bolton & Masters (2001) shown as shaded areas (SMAD curves for
the synthetic ‘first upswing’ type of measurements are shown for
completeness). To compute the range given by the observations, we
followed the philosophy of Bolton & Masters (2001), who tried to
isolate the part of the standard deviation σ₃⁻ₔ that is due to 3-D
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Comparison of the standard deviation of traveltime variations in model S09-M2-Q to that part of the standard deviation of the observations only that can be attributed to 3-D mantle heterogeneity alone. Intermediate and light shaded areas show the range of values inferred from the data of Bolton & Masters (2001). A, B and C denote measurements of different quality (see Section 3.2.2 for details). Blue lines: scaled median average deviation (SMAD) of the synthetic *P*-wave traveltime variations. Red lines: same for *S* waves. SMAD curves are shown for two types of measurements: (solid lines) full waveform cross-correlation measurements and (dashed lines) 'first upswing' measurements (see Section 2.5 for details). Note that we focus on the cross-correlation traveltime variations, as their corresponding sensitivity (i.e. banana-doughnut kernel) is well defined and allows for a direct and intuitive interpretation of the results.</figcaption></figure>

<table><caption>**Table 2.** Estimates of the standard deviation in traveltime variations due to location errors ($\sigma_X$) and additional random errors ($\sigma_N$) that contribute to the total variance of the observations [taken from Bolton & Masters (2001)]. $\sigma_X^l$ and $\sigma_X^u$ are the lower and upper estimates of $\sigma_X$, and A, B and C denote groups of traveltime measurements with different quality.</caption><thead><tr><th></th><th>σ<sub>X</sub><sup>l</sup></th><th>σ<sub>X</sub><sup>u</sup></th><th>σ<sub>N</sub><sup>A</sup></th><th>σ<sub>N</sub><sup>B</sup></th><th>σ<sub>N</sub><sup>C</sup></th></tr></thead><tbody><tr><th>P</th><td>0.6</td><td>1.2</td><td>0.5</td><td>0.9</td><td>1.3</td></tr><tr><th>S</th><td>1.6</td><td>2.5</td><td>0.5</td><td>1.1</td><td>2.2</td></tr></tbody></table>

mantle heterogeneity by separating the total standard deviation $\sigma_T$ into different contributions (depth-integrated values in their case).

$$ \sigma_T^2 = \sigma_{3-D}^2 + \sigma_X^2 + \sigma_N^2. \quad (1) $$

Here, $\sigma_X$ and $\sigma_N$ are the standard deviations related to location errors and additional random errors (e.g. incorrect instrument response modelling, unidentified station timing problems, incorrect crustal corrections, random picking errors, etc.), respectively. Bolton & Masters (2001) estimated $\sigma_{3-D}$ using a summary ray analysis and, based on a thorough analysis of event mislocation errors, they provided a range of values for $\sigma_X$ for both *P* and *S* waves. Finally, they solved for $\sigma_N$ and obtained three estimates for *P* and *S* each, depending on the quality assigned to their measurements (termed categories A, B and C). Their different estimates for $\sigma_X$ and $\sigma_N$ are summarized in Table 2.

For the comparison to our synthetic depth-dependent SMAD curves, we followed the philosophy of Bolton & Masters (2001),

but in the reverse sense, and computed the range of possible $\sigma_{3-D}$ in the observed data from the original (unprocessed) SMAD profiles given in fig. 10 of their paper by

$$ \sigma_{3-D}(d) = \sqrt{\sigma_T^2(d) - \sigma_X^2 - \sigma_N^2}, \quad (2) $$

where *d* is the ray-theoretical turning depth and $\sigma_X$ and $\sigma_N$ are the depth-independent estimates provided by Bolton & Masters (2001). The random error contribution $\sigma_N$ can be assumed to be depth-independent, and we can also assume that $\sigma_X$ varies only slightly as a function of depth (i.e. epicentral distance). The latter assumption is justified because the waves all have close to vertical take-off and incidence angles. We finally infer two ranges of values for $\sigma_{3-D}$ from the observations for both *P*- and *S*-wave measurements: for the wider range, labelled A–C in Fig. 7, the lower bounds are calculated by using the largest estimate for $\sigma_X$ and the estimate of category C for $\sigma_N$, whereas the upper bounds are calculated with the smallest $\sigma_X$ and the category A estimate for $\sigma_N$; for the narrower range, labelled <ABC> in Fig. 7, lower and upper bounds are calculated using the largest and smallest estimates of $\sigma_X$, respectively, together with the weighted mean for $\sigma_N$ (according to their relative occurrence).

From Fig. 7, one can clearly see that the standard deviations of *P*- and *S*-wave measurements inferred from the observations show rather different trends as a function depth. The range of values for *P* waves is centred around 1 s and remains almost constant with depth, whereas that of the *S* waves shows a strong increase with depth from around 2.2 s at 1000 km depth to 4.8 s close to the CMB, with the strongest increase happening below around 2000 km depth.

Looking at the standard deviations for our isochemical geodynamic model, one can see that the SMAD curve for synthetic *P*-wave traveltime variations also shows values of around 1 s throughout the mantle and falls within the narrower range given by the data, except near the very top of the lower mantle (between around 700 and 900 km depth). There, however, it still stays within the wider data bounds. The SMAD values for the synthetic *S*-wave delay times are slightly greater than 2 s at the top of the lower mantle and increase strongly to values of more than 4 s close to the CMB. The *S*-wave curve also falls within the narrower data range in most of the lower mantle (between 900 and 2200 km depth). Below 2200 km depth, it also increases strongly towards the CMB, but slightly less so compared to the observed traveltime variations. Below the depth of around 2500 km, the cross-correlation SMAD values are outside even the wider data range, but only by about 0.2 s (and only about ~0.7 s lower than the median of the data range).

It is a remarkable result that the synthetic traveltime variations from an isochemical MCM reproduce the general trend of both *P*- and *S*-wave data very well (i.e. constant with depth for *P* and strongly increasing with depth for *S*). This is not directly to be expected from the similarity of the seismic velocity variations $\partial \ln v_p$ and $\partial \ln v_s$ inherent in this type of model. In the following, we will discuss possible reasons for why the traveltime variations of *P* and *S* waves show such a different behaviour. Furthermore, we will discuss the significance of the fact that the *S*-wave variations in the geodynamic model have slightly lower standard deviation close to the CMB than the one of the observations.

## 4 DISCUSSION

In Section 3, we presented a first data set of synthetic long-period traveltime variations for an MCM and a comparison of their statistics with those of the observed delay times. Synthetic, in our case,
has a twofold meaning in that we computed synthetic waveforms for a synthetic Earth, resulting in seismograms that were obtained independently of seismic observations. The numerical modelling techniques used here are already at a rather high level of robustness today (Bunge 1996; Komatitsch & Tromp 2002a,b; Stixrude & Lithgow-Bertelloni 2005b; Peter et al. 2011). Still, our study represents only the first step into a new direction and there are a number of limitations that need careful consideration. Nevertheless, we hope that the work presented here demonstrates the potential of our approach, which may serve as a complementary tool to the geodynamic interpretation of tomographic models in future.

## 4.1 Seismic heterogeneity as seen by finite-frequency waves

The most important result of our study is that isochemical whole mantle flow with strong core heating and a pyrolite composition is largely consistent with the statistics of body wave data in terms of the magnitude of traveltime variations. In particular, the standard deviations of our *P*- and *S*-wave delay times as a function of depth show the same marked difference as is seen in the observations: *P*-wave traveltime variations with a turning point in the lower mantle have a standard deviation that stays almost constant with depth, whereas the standard deviation of the *S*-wave delay times shows an increase with turning depth, which is most prominent below around 2000 km depth. This is a rather intriguing result, as the perturbations in *P*- and *S*-wave velocities are perfectly correlated in our isochemical model, and the ratio $R = \partial \ln v_s / \partial \ln v_p$ is a very weak function of depth being ~1.8 at the top of the lower mantle and slightly above 2 close to the CMB.

The difference between the SMAD curves of *P*- and *S*-wave delay times observed not only in the data of Bolton & Masters (2001), but also in this synthetic study, suggests that the structure actually seen by the two wave types is not necessarily the same. In other words, the sensitivity of the traveltimes of *P* and *S* waves of same frequency to changes in the medium is not the same, because the Fréchet kernels for *P* waves have a cross-sectional area that is almost twice that of the *S*-wave kernels. Thus, as mentioned earlier, wave front healing may be more serious for *P* waves depending on the size of the anomalies and the distance that the wave propagated (e.g. Hung *et al.* 2001; Nolet 2008; Malcolm & Trampert 2011). The good match of the standard deviations of our synthetic traveltime variations to those of the seismic observations seems to indicate that the amount of wave front healing affecting our 15 s body waves during their passage through the MCM is to first-order correct. From this, one may conclude that in addition to the strength of heterogeneity, also the characteristic length scales of structure in model S09-M2-Q are of the right order of magnitude compared to Earth.

One possibility to confirm our interpretation that the difference in the statistics is related to finite-frequency effects could be to compare the traveltimes measured at two different frequencies chosen such that the volumes of sensitivity would be comparable. However, the full wavefield, such as the one computed here with spectral element techniques, is more complex than captured by finite-frequency theory. There, typically the near-field effects are neglected as well as *S*-to-*P* and *P*-to-*S* conversions, which lead to additional sensitivity close to the source and receiver in case of *S* waves (e.g. Favier *et al.* 2004). In case of *P* waves, the sensitivity of the traveltimes to variations in shear wave velocity may not be negligible any more at longer periods (Zhang & Shen 2008). A detailed analysis of these

effects is not straightforward and thus beyond the scope of this study.

Our results bear on the importance of chemical heterogeneity in the lower mantle: Fig. 7 shows that isochemical whole mantle flow with a pyrolite composition and high core heating is able to account for most of the variance of the observed traveltime variations of both *P* and *S* waves. Adding a compositional component to the seismic heterogeneity predicted by MCMs will probably increase the standard deviation of traveltime variations depending on their size (cf. discussion on wave front healing earlier). Iron enrichment in the lowermost mantle, for example, would result in even lower shear wave velocities (e.g. Wang & Weidner 1996; Jackson 1998), thereby probably increasing positive traveltime delays. A modest chemical heterogeneity component might thus be a candidate to increase the fit of the synthetic SMAD values in Fig. 7. However, there are a number of other factors that strongly influence the standard deviation of traveltime variations, which need careful consideration, first.

## 4.2 Factors influencing the standard deviation of synthetic traveltime variations

In Section 2.4, we already noted that cross-correlation traveltime measurements show slightly lower standard deviation compared to the ‘first upswing’ measurements. As can be seen from Fig. 7, a small part of the differences between our cross-correlation SMAD curves and those inferred from the observations can be explained by the different measurement techniques. Taking the error estimates that went into the calculation of $\sigma_X$ from Bolton & Masters (2001) at face value, their estimates for $\sigma_N$ should account for most of the errors contributing to the variance in the observed data including the crustal corrections. Thus, we are left with the following two factors that most likely influence our comparison in Fig. 7, which will be described below:

(1) Data coverage;

(2) Limitations of the mineralogical model.

(1) The fact that the data coverage is different between our synthetic data set and the measurements by Bolton & Masters (2001) is one factor that may contribute to the slight differences in the lowermost mantle between synthetic and observed SMAD values for *S* waves seen in Fig. 7. In this study, we have shown that it is possible to obtain an almost homogeneous illumination of our synthetic structure using a global distribution of 17 earthquakes together with a virtual network of 42 250 equidistant stations. In the comparison to the observations, we assumed that the ensemble of data used by Bolton & Masters (2001) is also representative for the global distribution of delay times. On Earth, however, both source and receiver locations are unequally distributed, resulting in laterally and radially inhomogeneous coverage. Also, while representative for its tectonic history, model S09-M2-Q differs from the Earth, as discussed in Section 2.1. Thus, apart from the fact that an exact replication of the real data distribution is computationally not feasible, the differences between our model and Earth’s mantle make such an effort superfluous. Our large synthetic data set allows, however, to test the effect of geographically biased sampling on the standard deviations by neglecting a part of the measurements: In oceanic regions, we randomly selected only one-third of the available virtual seismic stations. We repeated this random sampling 50 times and combined the remaining oceanic stations with the full set of stations located in continental regions. Fig. 8 shows that the
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Same as Fig. 7. Dark coloured areas show the variation of the standard deviation in our synthetic data due to a geographically biased coverage: 50 random subsets have been drawn from our data set, each time considering only 1/3 of the stations located in oceanic regions. Here, only SMAD curves calculated from the cross-correlation measurements are shown.</figcaption></figure>

resulting variation of SMAD values given by these 50 subsets of the data is rather small, and comparison to Fig. 7 reveals that the standard deviations are not much different in case of a more realistic coverage. We found that even neglecting all stations in oceanic regions gives similar standard deviations, which strengthens our finding that model S09-M2-Q is compatible with the variance of observed body wave traveltime variations.

(2) The conversion of temperatures to seismic velocities needed in our approach is based on the mineralogical model of Stixrude & Lithgow-Bertelloni (2005b, 2007). This model includes most of the chemical components (~98 per cent) and a large number of mineral phases presumably existing in the mantle. However, it does not incorporate some of the most recent developments in mineral physics: the post-perovskite phase transition, sodium oxides (Na₂O), the effect of water as well as the iron spin-transition.

Introducing sodium into the system will probably be most relevant for the upper mantle and will not change the seismic velocities dramatically in case of pyrolite (Stixrude & Lithgow-Bertelloni 2005a). The abundance of water in the lower mantle is to date still unclear (Karato 2011), but recent studies estimate it to be much lower than in the transition zone or upper mantle (Inoue et al. 2010). Introducing water into the mineralogical model may thus probably have the most notable influence on seismic observations in the transition zone by changing the pressure (i.e. depth) of phase boundaries (e.g. Ohtani & Litasov 2006) and likely also by increasing the sensitivity of the seismic velocities to temperature as well as increasing the intrinsic attenuation. The effects of water in the transition zone may thus influence the SMAD values at the top of the lower mantle (between 700 and 900 km depth; cf. Fig. 7). To date, however, experimental data on hydrous mantle phases is still sparse, and (to the best of our knowledge) no mineralogical model exists so far that includes water. The iron spin-transition may potentially affect the lower mantle, but it is still a matter of debate whether it is seismically visible (Speziale et al. 2007; Stackhouse et al. 2007; McCammon et al. 2008; Cammarano et al. 2010; Caracas et al. 2010; Antonangeli et al. 2011).

From the above considerations, it is clear that our conclusions may be affected most by the possible existence of post-perovskite. An additional phase transition close to the CMB will lead to strong lateral variations in seismic velocities, especially on shorter scales, due to the non-linear effects of temperature near phase boundaries; that is, the morphological part of the temperature sensitivity becomes large compared to the isomorphic contribution (Ricard et al. 2005; Stixrude & Lithgow-Bertelloni 2007). In case of the transition from perovskite to post-perovskite, the shear wave velocity variations will be large across the transition, whereas the P-wave velocity will change only slightly (e.g. Stackhouse et al. 2005; Hirose 2006). This difference in its effect on P- and S-wave velocities makes the incorporation of post-perovskite the most likely candidate to increase the standard deviation of our synthetic S-wave traveltime variations without pushing the SMAD of the P-wave delay times out of the range given by the data. This remains to be confirmed by future tests using mineralogical models, which include the thermodynamic and elastic properties of post-perovskite, such as, for example, those of Stixrude & Lithgow-Bertelloni (2011) and Zunino et al. (2011).

In addition to including post-perovskite, the 1-D Q profile used for the anelasticity correction in the temperature to seismic velocity conversion, and especially its variation with depth in the lower mantle, will also have some influence on the exact shape of the standard deviation curves. Aelasticity has a stronger effect on variations in shear wave velocity than on those of the compressional wave velocity, in particular at high temperatures. This effect will make the sensitivity of the S-wave velocity to hot anomalies larger near the CMB, which could also be causing the differences in the standard deviation curves. We tested this by running one additional simulation for seismic velocity variations calculated without the anelasticity correction. The difference between purely elastic and anelastic velocities has a stronger effect on the standard deviation of S wave compared to P-wave delay times, as expected. Most important for our interpretation, however, the striking difference in their behaviour with depth still exists even in case of using the purely elastic sensitivities. Nevertheless, the effects of Q on the statistics of body wave traveltimes will also have to be analysed in more detail in future studies.

The standard deviation of traveltime variations may potentially be further influenced by the effects of melting close to the CMB. However, the size of ultra-low velocity zones (ULVZs), which is on the order of kilometres to several tens of kilometres as mapped out using high-frequency data (Thomas et al. 2009; Rost et al. 2010), may be too small to affect waves at 15 s period.

Our MCM S09-M2-Q is far from representing the true structure of the Earth. However, except for the incorrect prediction of the location of structures in the lowermost mantle, it compares well to tomographic shear wave velocity models, especially in terms of the magnitude of velocity anomalies. This was the basis for this study, and testing the model S09-M2-Q directly against the statistics of observed traveltime variations confirms our earlier finding: isochemical whole mantle flow with strong core heating resulting from a large thermal gradient across D'' is compatible not only with tomographic S-wave models, but also with seismic observations for both P and S waves. One could imagine to further increase the CMB temperature to increase our synthetic SMAD values in the lowermost mantle, where they are slightly lower than those of the observations. The CMB temperature of 4200 K used to construct model S09-M2-Q is still somewhat below the upper bound of the range of suggested values (4380 ± 350 K, Terasaki et al. 2011). However, the magnitude of traveltime variations will depend not only on the strength of the heterogeneity, but also on its
characteristic length scales, as mentioned before, which in turn depend on the plate reconstruction used as surface boundary condition and the particular choice of rheology. Future studies will need to explore these factors in conjunction with updated mineralogical models including post-perovskite to further constrain the relative contributions of thermal anomalies and chemical variations to the observed seismic heterogeneity in the mantle.

Finally, we want to note that our predicted *P*- and *S*-wave structures correlate perfectly throughout the mantle and thus our model cannot explain any anticorrelation of $\partial \ln v_s$ and $\partial \ln v_\phi$, which has been observed in the tomographic model of Masters *et al.* (2000). Note, however, that the latter has been built using ray theory on the basis of the data set of Bolton & Masters (2001, i.e. measured on long-period *P* and *S* waves at the same frequency). This means that wavefield effects, such as wave front healing, were not taken into account. As we have seen, these effects are different for *P* and *S* waves of same frequency, rendering a direct comparison of their traveltimes problematic. It is interesting to note that also Kennett *et al.* (1998), using short-period ISC first onset picks for both *P* and *S* waves, find regions with significant anticorrelation between $\partial \ln v_\phi$ and $\partial \ln v_s$, which, however, are located not within, but rather outside the two large superplumes. Recently, Malcolm & Trampert (2011) have shown in a numerical study that wave front healing introduces a strong bias in tomographic inversions for $\partial \ln v_s$ and $\partial \ln v_\phi$, which may result in anticorrelation of the two even in case of perfect correlation of *P*- and *S*-wave heterogeneity in the input model. It remains to be seen, whether the strong increase of the ratio $R = \partial \ln v_s / \partial \ln v_p$ with depth in the mantle, as inferred directly from the data by Bolton & Masters (2001), can also be explained by wavefield effects.

## 5 CONCLUSIONS

We have presented synthetic seismograms and a corresponding data set of traveltime variations for a synthetic mantle structure derived from geodynamic considerations. Seismic heterogeneity is predicted from the temperature field of an MCM, which is converted to seismic velocities and density using a thermodynamically self-consistent model of mantle mineralogy for a pyrolite composition. Synthetic seismograms for periods down to 10 s were computed for the predicted structure using an SEM. Altogether, we simulated the seismic wavefield for 17 earthquakes distributed evenly over the globe. The wavefield of each earthquake is 'recorded' by a very large number of virtual seismic stations to achieve a relatively homogeneous illumination of our model even with a low number of seismic sources. In total, we obtain each ~350 000 cross-correlation traveltime measurements at a dominant period of 15 s for *P* and *S* waves. By combining forward modelling techniques from geodynamics, mineral physics and seismology, we have set up a new approach to assess geodynamic hypotheses directly against seismic observations, which may serve as a complementary tool to tomographic inversions. In this study, we used this approach to test the assumption of a large thermal gradient in $D''$ and corresponding high core heat flow against the statistics of long-period *P*- and *S*-wave traveltime variations. Comparing the geodynamic predictions to not only *S*, but also *P*-wave information allows a more robust assessment of the geodynamic model as was possible before. Interestingly, the data set of Bolton & Masters (2001), which we use for comparison to our synthetic data set, shows some peculiar behaviour of the statistics of traveltime delays as a function of turning depth of the waves: the standard deviation of *P*-wave traveltime varia-

tions stays almost constant with depth, whereas that of the *S*-wave traveltimes increases strongly with depth towards the CMB. This increase in case of *S* waves is particularly strong below a depth of around 2000 km. Such a difference between *P* and *S* waves cannot be explained in a ray-theoretical framework under the assumption of a chemically homogeneous mantle.

Using full wavefield simulations, we find, however, that the standard deviations of *P*- and *S*-wave traveltime variations in our geodynamic model show the same differing behaviour as the observations. This is a remarkable result in light of the isochemical nature of our MCM and highlights the importance of taking the correct physics of wave propagation into account in the interpretation of long-period seismic data. Most importantly, our comparison shows that isochemical whole mantle flow with strong core heating and a pyrolite composition is capable of explaining the statistics of seismic observations. The standard deviations of our synthetic *P*- and *S*-wave traveltimes do not only show different trends with depth, but are also matching those of the observations well in terms of their magnitude, when contributions from errors in the data are taken into account. Although this finding does not necessarily mean that there is no chemical heterogeneity present in the lower mantle, it shows that variations in chemical composition are not required by the data set studied here.

The remaining small differences between the predicted and observed standard deviations of *S*-wave traveltime variations close to the CMB are most likely related to the fact that post-perovskite is not included in the mineralogical model used here. To better constrain the relative importance of thermal anomalies and variations in chemical composition for generating seismic heterogeneity, the comparisons presented here need to be repeated with post-perovskite included in the mineralogical conversion before any remaining difference between geodynamic predictions and the data can be attributed to a contribution from chemical heterogeneity.

## ACKNOWLEDGMENTS

The authors are grateful to Diego Mercerat for fruitful discussions during the preparation of this work as well as to Yanick Ricard and an anonymous reviewer for constructive comments on the original version of the manuscript. BSAS was supported by 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). The authors thank the Leibniz Supercomputing Centre for access to computing resources on HLRBII and 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.

## REFERENCES

* Alfè, D., Price, G.D. & Gillan, M.J., 2002. Iron under Earth's core conditions: liquid-state thermodynamics and high-pressure melting curve from ab initio calculations, *Phys. Rev. B*, **65**(16), 165118-165129.

* Alfè, D., Gillan, M.J. & Price, G.D., 2007. Temperature and composition of the Earth's core, *Contemp. Phys.*, **48**(2), 63-80.

* Antonangeli, D. *et al.*, 2011. Spin crossover in ferropericlase at high pressure: a seismologically transparent transition? *Science*, **331**(6013), 64-67.

* Bassin, C., Laske, G. & Masters, G., 2000. The current limits of resolution for surface wave tomography in North America, *EOS, Trans. Am. geophys. Un.*, **81**, F897 Fall Meet. Suppl., Abstract S12A-03.
Boehler, R., 2000. High-pressure experiments and the phase diagram of lower mantle and core materials, *Rev. Geophys.*, **38**(2), 221–245.

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 mantles, *J. geophys. Res.*, **106**(B7), 13 527–13 540.

Boschi, L. & Dziewonski, A.M., 1999. High- and low-resolution images of the Earth's mantle: implications of different approaches to tomographic modeling, *J. geophys. Res.*, **104**(B11), 25 567–25 594.

Brodholt, J.P., Helffrich, G. & Trampert, J., 2007. Chemical versus thermal heterogeneity in the lower mantle: the most likely role of anelasticity, *Earth planet. Sci. Lett.*, **262**(3–4), 429–437.

Buffett, B.A., 2002. Estimates of heat flow in the deep mantle based on the power requirements for the geodynamo, *Geophys. Res. Lett.*, **29**(12), 1566.

Bunge, H.-P., 1996. Global mantle convection models, *Ph.D. thesis*, University of California, Berkeley.

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. & Davies, J. H., 2001. Tomographic images of a mantle circulation model, *Geophys. Res. Lett.*, **28**(1), 77–80.

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.

Bunge, H.-P., Ricard, Y. & Matas, J., 2001. Non-adiabaticity in mantle convection, *Geophys. Res. Lett.*, **28**(5), 879–882.

Bunge, H.-P., Richards, M. & Baumgardner, J., 2002. Mantle circulation models with sequential data-assimilation: inferring present-day mantle structure from plate motion histories, *Phil. Trans. R. Soc. Lond.*, **A**, **360**(1800), 2545–2567.

Bunge, H.-P., Hagelberg, C.R. & Travis, B. J., 2003. Mantle circulation models with variational data-assimilation: inferring past mantle flow and structure from plate motion histories and seismic tomography, *Geophys. J. Int.*, **2**(152), 280–301.

Cammarano, F., Goes, S., Vacher, P. & Giardini, D., 2003. Inferring upper-mantle temperatures from seismic velocities, *Phys. Earth planet. Inter.*, **138**(3–4), 197–222.

Cammarano, F., Deuss, A., Goes, S. & Giardini, D., 2005. One-dimensional physical reference models for the upper mantle and transition zone: combining seismic and mineral physics constraints, *J. geophys. Res.*, **110**, B01306, doi:10.1029/2004JB003272.

Cammarano, F., Marquardt, H., Speziale, S. & Tackley, P.J., 2010. Role of iron-spin transition in ferropericlase on seismic interpretation: a broad thermochemical transition in the mid mantle? *Geophys. Res. Lett.*, **37**(3), L03308, doi:10.1029/2009GL041583.

Capdeville, Y., Chaljub, E., Vilotte, J.P. & Montagner, J.P., 2003. Coupling the spectral element method with a modal solution for elastic wave propagation in global Earth models, *Geophys. J. Int.*, **152**, 34–67.

Caracas, R., Mainprice, D. & Thomas, C., 2010. Is the spin transition in Fe²⁺-bearing perovskite visible in seismology? *Geophys. Res. Lett.*, **37**(13), L13309, doi:10.1029/2010GL043320.

Chaljub, E. & Valette, B., 2004. Spectral element modelling of three-dimensional wave propagation in a self-gravitating Earth with an arbitrarily stratified outer core, *Geophys. J. Int.*, **158**, 131–141.

Cobden, L., Goes, S., Cammarano, F. & Connolly, J.A.D., 2008. Thermochemical interpretation of one-dimensional seismic reference models for the upper mantle: evidence for bias due to heterogeneity, *Geophys. J. Int.*, **175**(2), 627–648.

Cobden, L., Goes, S., Ravenna, M., Styles, E., Cammarano, F., Gallagher, K. & Connolly, J.A.D., 2009. Thermochemical interpretation of 1-D seismic data for the lower mantle: the significance of nonadiabatic thermal gradients and compositional heterogeneity, *J. geophys. Res.*, **114**, B11309, doi:10.1029/2008JB006262.

Crotwell, H.P., Owens, T.J. & Ritsema, J., 1999. The TauP toolkit: flexible seismic travel-time and ray-path utilities, *Seismol. Res. Lett.*, **70**(2), 154–160.

Dahlen, F.A., 2004. Resolution limit of traveltime tomography, *Geophys. J. Int.*, **157**(1), 315–331.

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

Dalton, C.A., Ekström, G. & Dziewonski, A.M., 2009. Global seismological shear velocity and attenuation: a comparison with experimental observations, *Earth planet. Sci. Lett.*, **284**(1–2), 65–75.

Davies, D. & Davies, J., 2009. Thermally-driven mantle plumes reconcile multiple hot-spot observations, *Earth planet. Sci. Lett.*, **278**(1–2), 50–54.

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

Deal, M. & Nolet, G., 1996. Nullspace shuttles, *Geophys. J. Int.*, **124**, 372–380.

Deschamps, F. & Trampert, J., 2004. Towards a lower mantle reference temperature and composition, *Earth planet. Sci. Lett.*, **222**(1), 161–175.

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

Dziewonski, A.M., Hager, B.H. & O’Connell, R.J., 1977. Large-scale heterogeneities in the lower mantle, *J. geophys. Res.*, **82**(2), 239–255.

Favier, N., Chevrot, S. & Komatitsch, D., 2004. Near-field influences on shear wave splitting and traveltime sensitivity kernels, *Geophys. J. Int.*, **156**, 467–482.

Fichtner, A., Kennett, B.L.N., Igel, H. & Bunge, H.-P., 2009. Full seismic waveform tomography for upper-mantle structure in the australasian region using adjoint methods, *Geophys. J. Int.*, **179**(3), 1703–1725.

Fichtner, A., Kennett, B.L., Igel, H. & Bunge, H.-P., 2010. Full waveform tomography for radially anisotropic structure: new insights into present and past states of the australasian upper mantle, *Earth planet. Sci. Lett.*, **290**(3–4), 270–280.

Glatzmaier, G.A. & Roberts, P.H., 1995. A 3-dimensional self-consistent computer-simulation of a geomagnetic-field reversal, *Nature*, **377**(6546), 203–209.

Goes, S., Cammarano, F. & Hansen, U., 2004. Synthetic seismic signature of thermal mantle plumes, *Earth planet. Sci. Lett.*, **218**(3–4), 403–419.

Grand, S., van der Hilst, R. & Widiyoranto, S., 1997. Global seismic tomography: A snapshot of mantle convection in the Earth, *GSA Today*, **7**, 1–7.

Gubbins, D., Alfé, D., Masters, G., Price, G.D. & Gillan, M., 2004. Gross thermodynamics of two-component core convection, *Geophys. J. Int.*, **157**(3), 1407–1414.

Hager, B.H. & Richards, M.A., 1989. Long-wavelength variations in Earths geoid—Physical models and dynamical implications, *Phil. Trans. R. Soc. Lond.*, **A**, **328**(1599), 309–327.

Hirose, K., 2006. Postperovskite phase transition and its geophysical implications, *Rev. Geophys.*, **44**, RG3001, doi:10.1029/2005RG000186.

Houser, C., Masters, G., Shearer, P. & Laske, G., 2008. Shear and compressional velocity models of the mantle from cluster analysis of long-period waveforms, *Geophys. J. Int.*, **174**(1), 195–212.

Hung, S.H., Dahlen, F.A. & Nolet, G., 2000. Frechet 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., Takeuchi, N., Geller, R.J., Megnin, C., Bunge, H.-P., Clévédé, E., Dalkolmo, J. & Romanowicz, B., 2000. The cosy project: verification of global seismic modeling algorithms, *Phys. Earth planet. Inter.*, **119**(1–2), 3–23.

Inoue, T., Wada, T., Sasaki, R. & Yurimoto, H., 2010. Water partitioning in the Earth’s mantle, *Phys. Earth planet. Inter.*, **183**(1–2), 245–251.
Ishii, M. & Tromp, J., 2001. Even-degree lateral variations in the Earth's mantle constrained by free oscillations and the free-air gravity anomaly, *Geophys. J. Int.*, **145**(1), 77–96.

Jackson, I., 1998. Elasticity, composition and temperature of the earth's lower mantle: a reappraisal, *Geophys. J. Int.*, **134**(1), 291–311.

Jacobs, M.H. & van den Berg, A.P., 2011. Complex phase distribution and seismic velocity structure of the transition zone: convection model predictions for a magnesium-endmember olivine-pyroxene mantle, *Phys. Earth planet. Inter.*, **186**(1–2), 36–48.

Karato, S., 1993. Importance of anelasticity in the interpretation of seismic tomography, *Geophys. Res. Lett.*, **20**, 1623–1626.

Karato, S., 2011. Water distribution across the mantle transition zone and its implications for global material circulation, *Earth planet. Sci. Lett.*, **301**(3-4), 413–423.

Kennel, M.B., 2004. KDTREE 2: Fortran 95 and C++ software to efficiently search for near neighbors in a multi-dimensional Euclidean spaces, arXiv:physics/0408067.

Kennett, B. & Jackson, I., 2009. Optimal equations of state for mantle minerals from simultaneous non-linear inversion of multiple datasets, *Phys. Earth planet. Inter.*, **176**(1–2), 98–108.

Kennett, B.L.N., Engdahl, E.R. & Buland, R., 1995. Constraints on seismic velocities in the Earth from travel-times, *Geophys. J. Int.*, **122**(1), 108–124.

Kennett, B.L.N., Widiyantoro, S. & van der Hilst, R.D., 1998. Joint seismic tomography for bulk sound and shear wave speed in the Earth's mantle, *J. geophys. Res.*, **103**(B6), 12469–12493.

Khan, A., Connolly, J.A.D. & Taylor, S.R., 2008. Inversion of seismic and geodetic data for the major element chemistry and temperature of the earth's mantle, *J. geophys. Res.*, **113**, B09308, doi:10.1029/2007JB005239.

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

Komatitsch, D. & Tromp, J., 1999. Introduction to the spectral-element method for 3-D seismic wave propagation, *Geophys. J. Int.*, **139**, 806–822.

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

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

Komatitsch, D., Vilotte, J.P., Vai, R., Castillo-Covarrubias, J.M. & Sánchez-Sesma, F.J., 1999. The spectral element method for elastic wave equations: application to 2D and 3D seismic problems, *Int. J. Numer. Methods Eng.*, **45**, 1139–1164.

Komatitsch, D., Barnes, C. & Tromp, J., 2000. Simulation of anisotropic wave propagation based upon a spectral element method, *Geophysics*, **65**(4), 1251–1260.

Komatitsch, D., Liu, Q., Tromp, J., Süss, P., Stidham, C. & Shaw, J.H., 2004. Simulations of ground motion in the Los Angeles basin based upon the spectral-element method, *Bull. seism. Soc. Am.*, **94**, 187–206.

Komatitsch, D., Vinnik, L.P. & Chevrot, S., 2010. Shdiff/svdiff splitting in an isotropic earth, *J. geophys. Res.*, **115**(B7), B07312, doi:10.1029/2009JB006795.

Kuang, W.L. & Bloxham, J., 1997. An Earth-like numerical dynamo model, *Nature*, **389**(6649), 371–374.

Kustowski, B., Ekstrom, G. & Dziewonski, A.M., 2008. Anisotropic shear-wave velocity structure of the Earth's mantle: a global models, *J. geophys. Res.*, **113**(B6), B06306, doi:10.1029/2007JB005169.

Lay, T., 2008. Sharpness of the D'' discontinuity beneath the Cocos plate: implications for the perovskite to post-perovskite phase transition, *Geo-phys. Res. Lett.*, **35**, L03304, doi:10.1029/2007GL032465.

Li, X.D. & Romanowicz, B., 1996. Global mantle shear velocity model developed using nonlinear asymptotic coupling theory, *J. geophys. Res.*, **101**(B10), 2245–2272.

Li, C., van der Hilst, R.D., Engdahl, E.R. & Burdick, S., 2008. A new global model for *P* wave speed variations in Earth's mantle, *Geochem. Geophys. Geosyst.*, **9**(5), Q05018, doi:10.1029/2007GC00180.

Lin, F.-C. & Ritzwoller, M.H., 2011. Apparent anisotropy in inhomogeneous isotropic media, *Geophys. J. Int.*, **186**(3), 1205–1219, doi:10.1111/j.1365-246X.2011.05100.x.

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.

Masters, G., Johnson, S., Laske, G., Bolton, H. & Davies, J.H., 1996. A shear-velocity model of the mantle, *Phil. Trans. R. Soc. Lond., A*, **354**(1711), 1385–1410.

Masters, G., Laske, G., Bolton, H. & Dziewonski, A.M., 2000. The relative behavior of shear velocity, bulk sound speed, and compressional velocity in the mantle: implications for chemical and thermal structure, in *Earth's Deep Interior: Mineral Physics and Tomography From the Atomic to the Global Scale*, pp. 63–87, American Geophysical Union, Washington DC.

Matas, J., Bass, J., Ricard, Y., Mattern, E. & Bukowinski, M.S.T., 2007. On the bulk composition of the lower mantle: predictions and limitations from generalized inversion of radial seismic profiles, *Geophys. J. Int.*, **170**, 764–780.

Mattern, E., Matas, J., Ricard, Y. & Bass, J., 2005. Lower mantle composition and temperature from mineral physics and thermodynamic modelling, *Geophys. J. Int.*, **160**(3), 973–990.

McCammon, C. et al., L., 2008. Stable intermediate-spin ferrous iron in lower-mantle perovskite, *Nat. Geosci.*, **1**(10), 684–687.

Mégnin, C., Bunge, H.-P., Romanowicz, B. & Richards, M., 1997. Imaging 3-D spherical convection models: what can seismic tomography tell us about mantle dynamics? *Geophys. Res. Lett.*, **24**(11), 1299–1302.

Mittelstaedt, E. & Tackley, P.J., 2006. Plume heat flow is much lower than CMB heat flow, *Earth planet. Sci. Lett.*, **241**(1–2), 202–210.

Montelli, R., Nolet, G., Dahlen, F.A. & Masters, G., 2006. A catalogue of deep mantle plumes: new results from finite-frequency tomography, *Geochem. Geophys. Geosyst.*, **7**, Q11007, doi:10.1029/2006GC001248.

Ni, S.D., Tan, E., Gurnis, M. & Helmberger, D.V., 2002. Sharp sides to the African superplume, *Science*, **296**, 1850–1852.

Nimmo, F., Price, G.D., Brodholt, J. & Gubbins, D., 2004. The influence of potassium on core and geodynamo evolution, *Geophys. J. Int.*, **156**(2), 363–376.

Nolet, G., 2008. *A Breviary of Seismic Tomography*, Cambridge Univ. Press, New York, NY, ISBN 978-0-521-88244-6.

Nolet, G., Karato, S.I. & Montelli, R., 2006. Plume fluxes from seismic tomography, *Earth planet. Sci. Lett.*, **248**, 685–699.

Oeser, J., Bunge, H.P. & Mohr, M., 2006. Cluster design in the earth sciences—Tethys, in *Proceedings of High Performance Computing and Communications*, 4208, pp. 31–40.

Ohtani, E. & Litasov, K.D., 2006. The effect of water on mantle phase transitions, *Rev. Mineral. Geochem.*, **62**(1), 397–420.

Panning, M. & Romanowicz, B., 2006. A three dimensional radially anisotropic model of shear velocity in the whole mantle, *Geophys. J. Int.*, **167**, 361–379.

Paulson, A., Zhong, S.J. & Wahr, J., 2007. Inference of mantle viscosity from GRACE and relative sea level data, *Geophys. J. Int.*, **171**(2), 497–508.

Peter, D. et al., 2011. Forward and adjoint simulations of seismic wave propagation on fully unstructured hexahedral meshes, *Geophys. J. Int.*, **186**(2), 721–739.

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.

Priolo, E. & Seriani, G., 1991. A numerical investigation of Chebyshev spectral element method for acoustic wave propagation, in *Proc. 13th IMACS Conf. on Comp. Appl. Math., v. 2*, Dublin, pp. 551–556.
Ricard, Y., Richards, M., Lithgow-Bertelloni, C. & LeStunff, Y., 1993. A geodynamic model of mantle density heterogeneity, *J. geophys. Res.*, **98**(B12), 21 895–21 909.

Ricard, Y., Mattern, E. & Matas, J., 2005. Synthetic tomographic images of slabs from mineral physics, in *Earth's Deep Mantle: Structure, Composition, and Evolution*, pp. 285–302, American Geophysical Union, Washington, D.C.

Ritsema, J. & van Heijst, H.J., 2002. Constraints on the correlation of P- and S-wave velocity heterogeneity in the mantle from P, PP, PPP and PKPab traveltimes, *Geophys. J. Int.*, **149**(2), 482–489.

Ritsema, J., van Heijst, H.J. & Woodhouse, J.H., 2004. Global transition zone tomography, *J. geophys. Res.*, **109**(B2), B02302, doi:10.1029/2003JB002610.

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.

Rost, S., Garnero, E.J. & Stefan, W., 2010. Thin and intermittent ultralow-velocity zones, *J. geophys. Res.*, **115**(B6), B06312, doi:10.1029/2009JB006981.

Savage, B., Komatitsch, D. & Tromp, J., 2010. Effects of 3D attenuation on seismic wave amplitude and phase measurements, *Bull. seism. Soc. Am.*, **100**(3), 1241–1251.

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/2009GC00240.

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.

Seriani, G., 1997. A parallel spectral element method for acoustic wave modeling, *J. Comput. Acoust.*, **5**, 53–69.

Seriani, G. & Priolo, E., 1991. High-order spectral element method for acoustic wave modeling, *SEG Expanded Abstracts*, **10**, 1561–1564.

Seriani, G. & Priolo, E., 1994. Spectral element method for acoustic wave simulation in heterogeneous media, *Finite Elem. Anal. Des.*, **16**, 337–348.

Seriani, G., Priolo, E., Carcione, J.M. & Padovani, E., 1992. High-order spectral element method for elastic wave modeling, *SEG Expanded Abstracts*, **11**, 1285–1288.

Shephard, G., Bunge, H.-P., Schuberth, B., Müller, R., Talsma, A., Moder, C. & Landgrebe, T., 2012. Testing absolute plate reference frames and the implications for the generation of geodynamic mantle heterogeneity structure, *Earth planet. Sci. Lett.*, **317**, 204–217, doi:10.1016/j.epsl.2011.11.027.

Sigloch, K., McQuarrie, N. & Nolet, G., 2008. Two-stage subduction history under North America inferred from multiple-frequency tomography, *Nat. Geosci.*, **1**(7), 458–462.

Simmons, N.A., Forte, A.M. & Grand, S.P., 2007. Thermochemical structure and dynamics of the African superplume, *Geophys. Res. Lett.*, **34**, L02301, doi:10.1029/2006GL028009.

Simmons, N.A., Forte, A.M. & Grand, S.P., 2009. Joint seismic, geodynamic and mineral physical constraints on three-dimensional mantle heterogeneity: implications for the relative importance of thermal versus compositional heterogeneity, *Geophys. J. Int.*, **177**(3), 1284–1304.

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.

Speziale, S., Lee, V.E., Clark, S.M., Lin, J.F., Pasternak, M.P. & Jeanloz, R., 2007. Effects of Fe spin transition on the elasticity of (Mg, Fe)O magnesiowüstites and implications for the seismological properties of the earth's lower mantle, *J. geophys. Res.*, **112**(B10), B10212, doi:10.1029/2006JB004730.

Stackhouse, S., Brodholt, J.P., Wookey, J., Kendall, J.-M. & Price, G.D., 2005. The effect of temperature on the seismic anisotropy of the perovskite and post-perovskite polymorphs of MgSiO₃, *Earth planet. Sci. Lett.*, **230**(1–2), 1–10.

Stackhouse, S., Brodholt, J.P. & Price, G.D., 2007. Electronic spin transitions in iron-bearing MgSiO₃ perovskite, *Earth planet. Sci. Lett.*, **253**(1–2), 282–290.

Steinberger, B. & Holme, R., 2008. Mantle flow models with core-mantle boundary constraints and chemical heterogeneities in the lowermost mantle, *J. geophys. Res.*, **113**(B5), B05403, doi:10.1029/2007JB005080.

Steinle-Neumann, G., Stixrude, L., Cohen, R.E. & Gulseren, O., 2001. Elasticity of iron at the temperature of the Earth's inner core, *Nature*, **413**(6851), 57–60.

Stixrude, L. & Lithgow-Bertelloni, C., 2005a. Mineralogy and elasticity of the oceanic upper mantle: origin of the low-velocity zone, *J. geophys. Res.*, **110**, B03204, doi:10.1029/2004JB002965.

Stixrude, L. & Lithgow-Bertelloni, C., 2005b. 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.

Styles, E., Davies, D.R. & Goes, S., 2011. Mapping spherical seismic into physical structure: biases from 3-D phase-transition and thermal boundary-layer heterogeneity, *Geophys. J. Int.*, **184**(3), 1371–1378, doi:10.1111/j.1365-246X.2010.04914.x.

Su, W.-J. & Dziewonski, A.M., 1997. Simultaneous inversion for 3-D variations in shear and bulk velocity in the mantle, *Phys. Earth planet. Inter.*, **100**(1–4), 135–156.

Su, W.-J., Woodward, R.L. & Dziewonski, A.M., 1994. Degree 12 model of shear velocity heterogeneity in the mantle, *J. geophys. Res.*, **99**(B4), 6945–6980.

Tape, C., Liu, Q., Maggi, A. & Tromp, J., 2009. Adjoint tomography of the southern california crust, *Science*, **325**(5943), 988–992.

Terasaki, H., Kamada, S., Sakai, T., Ohtani, E., Hirao, N. & Ohishi, Y., 2011. Liquidus and solidus temperatures of a Fe-O-S alloy up to the pressures of the outer core: implication for the thermal structure of the Earth's core, *Earth planet. Sci. Lett.*, **304**(3–4), 559–564.

Thomas, C., Kendall, J.-M. & Helffrich, G., 2009. Probing two low-velocity regions with PKP b-caustic amplitudes and scattering, *Geophys. J. Int.*, **178**(1), 503–512.

Tian, Y., Sigloch, K., Nolet, G. & Laske, G., 2011. Structure of North American mantle constrained by simultaneous inversion of multiple-frequency SH, SS, and Love waves, *J. geophys. Res.*, **116**, B02307, doi:10.1029/2010JB007704.

To, A., Romanowicz, B., Capdeville, Y. & Takeuchi, N., 2005. 3D effects of sharp boundaries at the borders of the African and Pacific superplumes: observation and modeling, *Earth planet. Sci. Lett.*, **233**(1–2), 137–153.

Trampert, J., Deschamps, F., Resovsky, J. & Yuen, D., 2004. Probabilistic tomography maps chemical heterogeneities throughout the lower mantle, *Science*, **306**(5697), 853–856.

Tromp, J., Tape, C. & Liu, Q.Y., 2005. Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels, *Geophys. J. Int.*, **160**(1), 195–216.

van der Hilst, R.D. & Karason, H., 1999. Compositional heterogeneity in the bottom 1000 kilometers of Earth's mantle: toward a hybrid convection model, *Science*, **283**(5409), 1885–1888.

van der Hilst, R.D., Widiyantoro, S. & Engdahl, E.R., 1997. Evidence for deep mantle circulation from global tomography, *Nature*, **386**(6625), 578–584.

van der Hilst, R.D., de Hoop, M.V., Wang, P., Shim, S.H., Ma, P. & Tenorio, L., 2007. Seismostratigraphy and thermal structure of Earth's core-mantle boundary region, *Science*, **315**, 1813–1817.

Wang, Y. & Wen, L., 2004. Mapping the geometry and geographic distribution of a very low velocity province at the base of the Earth's mantle, *J. geophys. Res.*, **109**(B10), B10305, doi:10.1029/2003JB002674.
Wang, Y.B. & Weidner, D.J., 1996. $(\partial\mu/\partial t)_P$ of the lower mantle, *Pure appl. Geophys.*, **146**(3–4), 533–549.

Wen, L.X., Silver, P.G., James, D. & Kuehnel, R., 2001. Seismic evidence for a thermo-chemical boundary at the base of the Earth's mantle, *Earth planet. Sci. Lett.*, **189**, 141–153.

Wielandt, E., 1987. On the validity of the ray approximation for interpreting delay times, in *Seismic Tomography*, D. Reidel, Norwell, MA.

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.

Zhang, Z. & Shen, Y., 2008. Cross-dependence of finite-frequency compressional waveforms to shear seismic wave speeds, *Geophys. J. Int.*, **174**(3), 941–948.

Zhong, S.J., 2006. Constraints on thermochemical convection of the mantle from plume heat flux, plume excess temperature, and upper mantle temperature, *J. geophys. Res.*, **111**(B4), B04409, doi:10.1029/2005JB003972.

Zunino, A., Connolly, J.A.D. & Khan, A., 2011. Precalculated phase equilibrium models for geophysical properties of the crust and mantle as a function of composition, *Geochem. Geophys. Geosyst.*, **12**, Q04001, doi:10.1029/2010GC003304.

## APPENDIX: IMPLEMENTATION OF MCMS IN WAVE PROPAGATION SOFTWARE

Computing whole waveform synthetic seismograms for the whole globe at relevant frequencies became feasible on a routine basis in recent years thanks to rapidly growing computational resources. However, it has long been unclear how to introduce geodynamic models into seismological forward simulations in an efficient and consistent manner. One part of the problem is related to the size of the numerical grids: the large number of model parameters from the mantle convection simulation need to be projected onto the even larger grid used for the simulation of wave propagation. This task has to be performed every time the underlying geodynamic hypotheses or the mineralogical/compositional model used to compute the seismic velocities change. Even more important, owing to the large meshes needed to resolve the wavefield on a global scale and the related high storage requirements, standard practice on modern supercomputers may change towards recomputing the mesh and assign the velocity model 'on the fly'. This means that the interpolation of parameters from one grid to the other may have to be repeated for every simulation.

The second part of the problem is related to the different nature of temperature and seismic velocity perturbations. A naive approach would be to first compute seismic velocity and density perturbations from the temperatures of the geodynamic model on the grid of the convection simulation, and then try to do an interpolation of the elastic parameters onto the grid used for simulation of wave propagation. However, as shown in Schuberth et al. (2009b), the seismic velocities show very strong heterogeneity in the vicinity of phase transitions due to their strong non-linear dependence on temperature (Ricard et al. 2005; Stixrude & Lithgow-Bertelloni 2007). To correctly capture this behaviour would require complicated high-order interpolation schemes. Temperature, on the other hand, is a much

smoother function of space, for which simple linear interpolation is sufficient. Thus, a more consistent way of implementing geodynamically derived velocity models is to first project the temperature field onto the wave propagation grid and to do the conversion to seismic velocities there. This way, the non-linearity of the temperature sensitivity of the seismic velocities and the complicated nature of resulting perturbations in the vicinity of phase transitions is taken into account consistently and represented correctly in the wave propagation simulations.

In this study, we are faced with the task to project the roughly 80 million temperature values of model S09-M2-Q onto the grid of 'SPECFEM3D_GLOBE' consisting of ~1.3 billion points in our case, for which efficient algorithms have to be used. As we are interested in a linear interpolation of temperature only, the main task is to find the nearest points of the geodynamic grid for each node of the 'SPECFEM3D_GLOBE' grid. We take advantage of the fact that the grid used in the geodynamic simulations of Schuberth et al. (2009b) uses the same lateral discretization in each radial layer (a triangular mesh derived from the regular icosahedron with 655 360 nodes in each of the 128 radial layers). The problem of searching the nearest neighbours can therefore be separated for the radial and lateral part, thereby dramatically reducing the computational effort. The search in radial direction is trivial, and we make use of an open-source implementation of a tree-search algorithm called `kdtree2` (Kennel 2004) for the lateral part. The temperature values at the three nearest vertices in both enclosing radial layers are then interpolated by bilinear interpolation using barycentric coordinates, followed by a simple linear interpolation in radial direction to obtain the temperature value at the current node of interest in the spectral element grid. This procedure is repeated for every gridpoint in the mantle, which makes up the largest part of the numerical grid. For the core, we use the seismic velocities and densities of the 1-D reference model 'PREM' (Dziewonski & Anderson 1981).

## SUPPORTING INFORMATION

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

**Movie.** This movie shows the wave propagation in the geodynamic model for an event in the Fiji Islands region (event 4 in Table 1). The wavefield is depicted by green and magenta colours together with the shear wave velocity variations in model S09-M2-Q, for which vertical cross-sections and iso-surfaces are shown on a blue to brownish colour scale ranging from -2 to 2 per cent. Surface topography is also shown for parts of the globe for geographic reference.

Please note: Wiley-Blackwell are 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 article.