Taming Snakes and Dragons in the Earth
By Alison Malcolm and Andreas Fichtner
As most of recorded human history had it, the Earth’s interior was ruled by giant snakes and dragons. Only with pioneering studies of surface-exposed rocks (J. Hutton, 1726–1797) did our planet come to be perceived as a complex superorganism that could be understood without the invocation of mystical creatures. The initially descriptive geosciences have gradually been complemented by physical models of plate tectonics, mantle convection, and seismic wave propagation, among many others. Over the past few decades, the tools of mathematics have been increasingly exploited to improve our understanding of these models, as well as to shape their development. As with many fields, continuously increasing computer power is producing new tools with which we can explore our planet. In recent years, computational science and applied mathematics have evolved—from mere tools, they have become major drivers of geoscientific progress. The ongoing mathematization of Earth system studies was clear in Boston at this year’s SIAM Conference on Computational Science and Engineering, where cutting-edge research was presented in numerous sessions on such diverse topics as 4D imaging, flow in porous media, adjoint methods, ice sheet modeling, uncertainty quantification, and climate simulation.
Seismic modeling and inversion are two of the principal Earth-science beneficiaries of advances in computational mathematics. New methods that harness the power of modern HPC allow us to address key problems related to these topics, including the assimilation of massive data sets, 3D frequency-domain wave propagation and imaging, nonlinear inversion with multiple scattering, multiscale seismic wavespeed tomography, uncertainty quantification, and the search for optimal misfit functionals for full-waveform inversion.
As discussed at CSE 2013, today’s seismic acquisition systems produce immense quantities of waveform data that are increasingly difficult to handle and to model. To reduce the computational burden of seismic waveform inversion, F.J. Herrmann has used ideas from compressive sensing to design an iterative method in which updates are computed from random data subsets via curvelet-domain sparsity promotion. In a closely related effort to minimize the number of forward simulations in waveform inversion, K. van den Doel and colleagues used experimental design techniques to develop optimal source encoding strategies that model multiple seismic sources simultaneously, without loss of information about subsurface structure.
Whether the seismic inverse problem should be solved in the time or the frequency domain continues to be a subject of active debate. The first quantitative and objective comparison of the time- and frequency-domain approaches was presented at CSE 2013 by R. Brossier, whose group concluded that 3D full-waveform inversion in the frequency domain is feasible and competitive. Additional impetus to frequency-domain modeling was provided by B. Engquist, who developed “sweeping preconditioners” for the Helmholtz equation, for which the number of operations scales approximately linearly in the number of unknowns.
Owing to the extreme computational requirements, uncertainty quantification in full-waveform inversion has long been considered infeasible. G. Stadler described a practical solution to this problem, based on second-order adjoints, which he and colleagues used to construct a low-rank approximation of the Hessian that carries information on the posterior covariance. The approximated posterior distribution can then be sampled to provide a range of plausible Earth models (see Figure 1).
Figure 1. Sample from the posterior probability density of local wave speed at 70 km depth, resulting from a synthetic full-waveform inversion. Receivers were located only in the northern hemisphere. The gray lines in the graphs are the prior point marginals, the blue lines the posterior point marginals. Red dots indicate the point value of the sample. In the northern hemisphere, the posterior variance is reduced significantly compared to the prior variance. Only small variance reductions are observed in the southern hemisphere, where no receivers are located.
Full-waveform inversion is generally posed as an optimization problem in which the difference between modeled and observed data is to be minimized. A key challenge in this approach comes from the non-convexity of this objective function, which leads to the requirement of an initial guess that is close to the optimal solution. More specifically, this initial model must accurately represent the low-frequency parts of the model. While researchers who study the Earth at a large scale (100s of km) generally have a reasonably good estimate for such a model, such models are not generally available at smaller scales (1-10 km).
The emphasis then shifts to the related problems of estimating appropriate starting models, and of developing objective functions that are less susceptible to local minima. Along these lines, W. Symes illustrated how linear modeling can be extended to exploit the redundancy in the data and how this extension can be merged into full-waveform inversion to recover the low-frequency part of the model and thus avoid local minima. L. Demanet and V. Jugnon discussed a method that also aims to address this problem, using correlations between different parts of the data set. Both approaches seem to indicate that it is in new ways of combining complementary techniques and ideas that we can expect to overcome this problem.
A similar thread was visible in other sessions addressing geophysical problems. We heard from A. Malcolm and A. Richardson and from S. Fomel that, in some cases, adaptations of older imaging techniques can be superior to current techniques, although in different contexts. Malcolm and Richardson showed how including multiply-scattered events separately in an approximate method can sometimes result in a better image than processing all events simultaneously, as is done in more complete methods. S. Fomel stepped even further back, to the simplest of all imaging techniques, known as time migration, and showed how it can be improved via techniques of modern geometry.
Unlike seismic exploration, for which abundant regularly sampled data are available, global seismology suffers from a strongly heterogeneous distribution of sources and receivers over a wide range of spatio-temporal scales (1–10,000 km, 1–1000 s). Assimilating these data into a consistent Earth model requires the development of multiscale methods for seismic inversion, as in the work presented by A. Fichtner, in which multiple nested inverse problems on various scales are solved simultaneously (see Figure 2).
Figure 2. Multiscale full-waveform inversion with resolution analysis. (a) Shear velocity at 75 km depth beneath Eurasia. (b) Zoom-in view of the crustal velocity structure beneath Anatolia at 20 km depth. (c) Resolution length in N–S direction at 75 km depth, derived with second-order adjoints. AM, Armorican Massif; APB, Alghero–Provençal Basin; BM, Bohemian Massif; CSVF, Central Slovakian Volcanic Field; EEC, East European Craton; HS, Hellenic Slab; KM, Kirsehir Massif; MC, Massif Central; MM, Menderes Massif; PB, Pannonian Basin; RG, Rhine Graben; TS, Tyrrhenian Sea; TTL, Tornquist–Teisseyre Line.
Progress in the mathematical and computational Earth sciences is extraordinarily rapid. Problems like the assimilation of massive seismic data sets into 3D full-waveform inversions would have been considered a luxury only a few years ago. The snakes and dragons of Hutton’s time are definitely tamed. Nevertheless, numerous challenges remain: Researchers in seismic imaging and inversion will have to abandon computationally less expensive acoustic approximations to fully exploit the wealth of information contained in seismic recordings. The transition to elastic modeling and inversion needs to go hand in hand with the development of powerful wave equation solvers that honor the complexities of the real Earth, including irregular topography, fluid–solid interfaces, anisotropy, and poro-elasticity. With increasing complexity in modeling and inversion come new challenges for uncertainty quantification. Current methods rely on local approximations of the misfit functional, and they are effective for single-parameter problems. Future methods for uncertainty quantification must be able to detect the presence of multiple local minima, and to operate in multi-parameter inversions that constrain not only isotropic wave speeds, but also anisotropy, attenuation, and density. Such methods, in turn, will alleviate many problems related to local minima of the misfit function.
Alison Malcolm is an assistant professor of geophysics in the Department of Earth, Atmospheric, and Planetary Sciences at MIT. Andreas Fichtner is an assistant professor in the Department of Earth Sciences at ETH Zurich.