SUMMARY We present the theory for and applications of Hamiltonian Monte Carlo (HMC) solutions of linear and nonlinear tomographic problems. HMC rests on the construction of an artificial Hamiltonian system where a model is treated as a high-dimensional particle moving along a trajectory in an extended model space. Using derivatives of the forward equations, HMC is able to make long-distance moves from the current towards a new independent model, thereby promoting model independence, while maintaining high acceptance rates. Following a brief introduction to HMC using common geophysical terminology, we study linear (tomographic) problems. Though these may not be the main target of Monte Carlo methods, they provide valuable insight into the geometry and the tuning of HMC, including the design of suitable mass matrices and the length of Hamiltonian trajectories. This is complemented by a self-contained proof of the HMC algorithm in Appendix A. A series of tomographic/imaging examples is intended to illustrate (i) different variants of HMC, such as constrained and tempered sampling, (ii) the independence of samples produced by the HMC algorithm and (iii) the effects of tuning on the number of samples required to achieve practically useful convergence. Most importantly, we demonstrate the combination of HMC with adjoint techniques. This allows us to solve a fully nonlinear, probabilistic traveltime tomography with several thousand unknowns on a standard laptop computer, without any need for supercomputing resources.

No abstract

No abstract

[1] We jointly invert local fundamental-mode and higher-order surface-wave phase-velocities for radial models of the thermo-chemical and anisotropic physical structure of the Earth's mantle to ∼1000 km depth beneath the North American continent. Inversion for thermo-chemical state relies on a self-consistent thermodynamic method whereby phase equilibria and physical properties (P-, S-wave velocity and density) are computed as functions of composition (in the Na 2 O-CaO-FeO-MgO-Al 2 O 3 -SiO 2 model system), pressure and temperature. We employ a sampling-based strategy to solve the non-linear inverse problem relying on a Markov Chain Monte Carlo method to sample the posterior distribution in the model space. A range of models fitting the observations within uncertainties are obtained from which any statistics can be estimated. To further refine sampled models we compute geoid anomalies for a collection of these and compare with observations, exemplifying a posteriori filtering through the use of additional data. Our thermo-chemical maps reveal the tectonically stable older eastern parts of North America to be chemically depleted (high Mg#) and colder (>200°C) relative to the active younger regions (western margin and oceans). In the transition zone the thermo-chemical structure decouples from that of the upper mantle, with a relatively hot thermal anomaly appearing beneath the cratonic area that likely extends into the lower mantle. In the lower mantle no consistent large-scale thermo-chemical heterogeneities are observed, although our results do suggest distinct upper and lower mantle compositions. Concerning anisotropy structure, we find evidence for a number of distinct anisotropic layers pervading the mantle, including transition zone and upper-most lower mantle.

[1] Here we discuss the nature of velocity heterogeneities seen in seismic tomography images of Earth's mantle whose origins and relation to thermochemical variations are yet to be understood. We illustrate this by inverting fundamental-mode and higher-order surface-wave phase velocities for radial models of the thermochemical and anisotropic structure of the mantle to 450 km depth. Dispersion data are linked to thermochemical parameters through a thermodynamic formalism for computing mantle mineral phase equilibria and physical properties. The inverse problem is solved using a probabilistic inference approach whereby robust uncertainty estimates are obtained. We find that both compositional and thermal anomalies are required if observations are to be satisfied. Mantle thermochemical variations extend to 250 km depth beneath western and central Australia and are characterized by increased Mg/Fe and Mg/Si values relative to surrounding mantle. Correlated herewith are thermal variations that closely follow surface tectonics. We also observe a strong contribution to lateral variations in structure and topography across the "410 km" seismic discontinuity from thermochemically induced phase transformations that appear much stronger than lateral variations immediately above and below. Inside the transition zone, a general decoupling of structure relative to that of the upper mantle occurs driven by a reversal in Mg/Si, while thermal anomalies are smoothed out. Comparison of presently derived shear-wave tomography models with other regional models is encouraging. Radial anisotropy is strongest at 150/200 km depth beneath oceanic/continental areas, respectively, and appears weak and homogeneous below. Finally, geoid anomalies are computed for a subset of sampled model and compared to observations.

The problem of inferring information about the Earth can be described as a data integration problem, where the solu tion is a probability distribution that combines all available information. The theory is conceptually simple, but appli cation in practice can be challenging. Probabilistic data integration requires that the information at hand can be quantified in the form of a probability distribution, either (a) directly through specification of an analytical descrip tion of a probability distribution or (b) indirectly through algorithms that can sample an often unknown probability distribution. Once all information has been quantified, efficient numerical algorithms are needed for inferring infor mation from the combined probability distribution. In the following, methods for probabilistic characterization of different kinds of geo-information are presented. Then a number of methods that allow inferring information from the probability distribution that combines all available information will be discussed. Straight forward application of classic sampling algorithms such as the rejection sampler and the Metropolis algorithm will in most cases lead to computationally intractable problems. However, a number of methods exist that can turn an otherwise intractable data integration problem into a manageable one.

We present a method to explore the effective nullspace of nonlinear inverse problems without Monte Carlo sampling. This is based on the construction of an artificial Hamiltonian system where a model is treated as a high‐dimensional particle. Depending on its initial momentum and mass matrix, the particle evolves along a trajectory that traverses the effective nullspace, thereby producing a series of alternative models that are consistent with observations and their uncertainties. Variants of the nullspace shuttle enable hypothesis testing, for example, by adding features or by producing smoother or rougher models. Furthermore, the Hamiltonian nullspace shuttle can serve as a tunable hybrid between deterministic and probabilistic inversion methods: Choosing random initial momenta, it resembles Hamiltonian Monte Carlo; requiring misfits to decrease along a trajectory, it transforms into gradient descent. We illustrate the concept with a low‐dimensional toy example and with high‐dimensional nonlinear inversions of seismic traveltimes and magnetic data, respectively.

We present a new methodology for inverting P‐to‐S receiver function (RF) waveforms directly for mantle temperature and composition. This is achieved by interfacing the geophysical inversion with self‐consistent mineral phase equilibria calculations from which rock mineralogy and its elastic properties are predicted as a function of pressure, temperature, and bulk composition. This approach anchors temperatures, composition, seismic properties, and discontinuities that are in mineral physics data, while permitting the simultaneous use of geophysical inverse methods to optimize models of seismic properties to match RF waveforms. Resultant estimates of transition zone (TZ) topography and volumetric seismic velocities are independent of tomographic models usually required for correcting for upper mantle structure. We considered two end‐member compositional models: the equilibrated equilibrium assemblage (EA) and the disequilibrated mechanical mixture (MM) models. Thermal variations were found to influence arrival times of computed RF waveforms, whereas compositional variations affected amplitudes of waves converted at the TZ discontinuities. The robustness of the inversion strategy was tested by performing a set of synthetic inversions in which crustal structure was assumed both fixed and variable. These tests indicate that unaccounted‐for crustal structure strongly affects the retrieval of mantle properties, calling for a two‐step strategy presented herein to simultaneously recover both crustal and mantle parameters. As a proof of concept, the methodology is applied to data from two stations located in the Siberian and East European continental platforms.