Simulation of 3-D seismic wave propagation in a synthetic Earth
Long-standing questions in the study of Earth's deep interior are about the origin of seismic heterogeneity and the nature of flow in the mantle. Understanding the dynamic behaviour is important as the flow drives plate tectonics and controls the way the Earth looses its heat. Thus, it is a crucial factor in tectonic modelling or in simulations of the geodynamo and the thermal evolution of the Earth. A better understanding of these aspects is of great societal importance. For example, the continuous drift of tectonic plates relative to each other results in a build up of stress at the plate boundaries. This stress can eventually exceed the yield stress of rock thus leading to (often disastrous) earthquakes.
Most of our knowledge on deep Earth structure and flow comes from the analysis of recordings of seismic waves that travel through the Earth after large
earthquakes. Similar to medical tomography, seismic tomography allows us to “see” the present day elastic structure of Earth’s mantle in 3-D. During the last two decades, a variety of such 3-D models were built from different datasets. However, seismic tomography can only provide a limited resolution; that is, only a blurred and distorted image of Earth's structure can be obtained. This is a consequence of errors in the data and the non-uniform distribution of earthquakes as well as seismic receiver locations, which mainly cover continental regions. A second problem inherent to any tomographic inversion is that an infinite number of models will
fit the data equally well, a problem known
as non-uniqueness of the solution.
Geochemical observations also offer important constraints on Earth’s structure and evolution, providing compel-ling evidence for compositional heterogeneity within the mantle. Geophysical techniques, including seismology, can provide detailed information on such heterogeneity. However, despite substantial improvements in seismolog-ical methods and a dramatic increase in the amount and quality of available data, interpretations of images from seismic studies remain hampered by trade-offs between thermal and chemical effects. As a consequence, the distribution of chemical heterogeneity within Earth’s mantle and, accordingly, its role in governing mantle dynamics, remains poorly understood.
Imaging Earth’s interior structure at higher resolution and the convergence towards consistent models serving as a 3-D seismic reference is currently one of the most pressing challenges of modern global seismology today. So far, tomographic studies only used the arrival times of a very limited set of seismic phases from all of the information that is contained in a seismogram. In addition, most studies were based on ray theory, thus completely neglecting the wave character of seismic disturbances.
In order 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. High expectations to gain new insight currently lie within numerical simulations
of wave propagation through complex three-dimensional structures.
Figure 1: Snapshots of the three-dimensional wavefield in one of our geodynamic models. 3-D global wave propagation was simulated for an earthquake in the Fiji Islands region using a spectral element technique. The wavefield is depicted by green and magenta colours together with the shear wave
velocity variations in the model, 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 
Modern computational tools for seismic
wave propagation incorporate a large range of physical phenomena and are able to produce synthetic datasets that show a complexity comparable to real observations. Also, computing whole waveform synthetic seismograms at relevant frequencies became feasible on a routine basis in recent years thanks to rapidly growing computational resources. However, it has long been not clear how to introduce geodynamic considerations into seismological forward simulations in an efficient and consistent manner, and how to benefit from expensive large-scale simulations for our understanding of deep Earth structure and dynamics. This was the motivation to develop a novel method, in which we generate synthetic 3-D mantle structures based on dynamic flow calculations that serve as input models in the simulation of seismic wave propagation.
Here, we present the results of our new multi-disciplinary approach that combines forward modelling techniques from geodynamics, mineral physics and seismology. The thermal state of Earth's mantle at present-day geologic time is predicted by 3-D high-resolution mantle circulation models using a finite-element method. The temperature field is then mapped to seismic velocities. For this task, we take advantage of recent progress in describing the state of dynamic Earth models in terms of elastic properties through thermodynamically self-consistent models of mantle mineralogy. The predicted models of seismic heterogeneity are then implemented in a spectral element code for the simulation of 3-D global wave propagation. Using state-of-the-art techniques to solve the wave equa-tion in 3-D heterogeneous media, this approach allows us to capture the correct physics of wave propagation.
Both, the geodynamic as well as the seismic simulations require large-scale high-performance calculations. The computational resources provided through the supercomputing platform HLRB II at the Leibniz Supercomputing Centre (LRZ) allowed for the first time to simulate seismic wave propagation in synthetic Earths; that is, we are now able to compute synthetic seismograms completely independent of seismic observations. This means that we can test geodynamic hypotheses directly against seismic observations, which may serve as a complementary tool to tomographic inversions. More specifically, it is for the first time possible to study frequency-dependent waveform effects, such as wavefront healing and focusing/defocusing in mantle structures with realistic length-scales; that is, in a physically consistent manner.
One specific question that we addressed with our joint forward modelling approach is the origin of two large regions of strongly reduced seismic velocities in the lowermost mantle (the so-called African and Pacific super-plumes). Several seismological observations are typically taken as an indication that the superplumes are being caused by compositional variations and that they are “piles” of material with higher density than normal mantle rock. However, a large number of recent geodynamic, mineralogical and also seismological studies argue for a strong thermal gradient across the core-mantle boundary (CMB) that might provide an alternative explanation through the resulting large lateral temperature variations. We tested the hypothesis whether the presence of such a strong thermal gradient in isochemical whole mantle flow is compatible with geophysical observations.
We have computed synthetic seismograms and a corresponding dataset of traveltime variations for one of our synthetic mantle models assuming a pyrolite composition for the mantle mineralogy. Synthetic seismograms for periods down to 10 seconds were computed for the predicted structure using 486 cores on HLRB II. Altogether, we simulated the seismic wavefield for 17 earthquakes distributed evenly over the globe. The wavefield of each earthquake was "recorded" by a very large number of virtual seismic stations in order to achieve a relatively homogeneous illumination of our model even with a low number of seismic sources. In total, we obtained each ~350,000 cross-correlation traveltime measurements at a dominant period of 15 s for compressional (P) and shear (S) waves.
The traveltimes of observed seismic waves, to which we compared our synthetic dataset, show a peculiar behav-iour of their statistics as a function of the turning depth of the waves in the mantle: The standard deviation of P-wave traveltime variations stays almost constant with depth, while that of the S-wave traveltimes increases strongly towards the CMB. This increase in case of S-waves is particularly strong below a depth of around 2000 km (cf. Fig. 2). Such a difference between P- and S-waves is not expected in a chemically homogeneous mantle in the framework of ray-theory, which forms the backbone of seismology.
Figure 2: Comparison of the standard deviation (SMAD = scaled median average deviation) of traveltime variations in our geodynamic model to that of the observations. Intermediate and light shaded areas show the range of values inferred from the data . Blue lines: simulated P-wave traveltime variations. Red lines: same for S-waves. Solid and dashed lines show SMAD curves for two different measurement techniques .
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 peculiar behaviour as the observations. This is a remarkable result in light of the isochemical nature of our mantle circulation model and highlights the importance of taking the correct physics of wave propagation into account in the interpretation of long-period seismic data. Most important, 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. While this finding does not necessarily mean that there is no chemical heterogeneity present in the lower mantle, it shows that complex large-scale variations in chemical composition are not required by the dataset studied here.
On-going Research / Outlook
One particular aspect that we currently strive to clarify in the elastic representation of geodynamic heterogeneity is the suspected existence of a mineral phase named post-perovskite (the high pressure polymorph of magnesium silicate MgSiO3). The elastic properties of the mantle will be affected by the existence of this phase, which would cause strong shear wave velocity variations across the phase transition while the velocity of compressional waves would change only slightly. To better constrain the relative importance of thermal anomalies and variations in chemical composition for generating seismic heterogeneity, the comparison presented here needs to be repeated with post-perovskite included in the mineralogical conversion.
Seismic datasets 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. In this study, we have focused on a single frequency band. In future, it will become important to understand also the multi-frequency da-tasets from a forward modelling perspective complementary to using them for tomographic inversions. A limitation that we still face is that memory requirements are huge in order to simulate seismic wave propagation up to the highest frequencies used in seismological studies on a global scale (around 1 Hz; that is, 1 s shortest period). SuperMUC, the successor system to HLRB II will certainly help to bring us closer to the goal of covering the whole range of frequencies relevant for deep Earth studies.
Bernhard Schuberth was supported by a Marie Curie Intra European Fellow-
ship within the 7th European Community
Framework Programme [FP7/2007-2013] under grant agreement no. 235861. Guust Nolet and Christophe Zarolie re-
ceived support from the European Research
Council (ERC Advanced grant 226837). The authors thank the Leibniz Super-
computing 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.
 Schuberth, B.S.A., Zaroli, C. & Nolet,
J. Int., 188(3), 1393-1412, 
 Bolton, H. & Masters, G., J. Geophys.
Res., 106(B7), 13,527–13,540, 
• Bernhard Schuberth
• Christophe Zaroli
• Guust Nolet
Observatoire de la
sité de Nice