Terrestrial ice bodies are important regulators of climate and sea level variations. They influence the water cycle, provide fresh water and energy for human society, and contribute to the living basis of numerous ecosystems. Understanding the structure and dynamics of land ice requires knowledge of its mass density, which is essential for ice core climatology and estimates of mass balance components, such as mass loss, ice discharge and surface melt. We combine densely sampled fiber-optic sensing data from strong serendipitous anthropogenic sources with Hamiltonian Monte Carlo sampling to extract direct seismic constraints on firn density (i.e. the transitional layer between fresh snow and glacial ice). Our approach avoids biases introduced by subjective regularization choices, does not require empirical scaling relations from seismic wave speeds to density, and provides reliable uncertainty estimates. We demonstrate that high-quality surface-wave overtone data can directly constrain density to around 100 m depth. Commonly used scaling relations from seismic wave speeds to density, however, fail to reproduce resolvable details of glacial density structure, and they tend to deviate from direct constraints on the order of ±10 %. Consequently, ice mass inferred from seismic wave speed may be incorrect by a similar amount.

Seismic tomography, a critical tool for studying Earth's interior structure and dynamics, has revealed positive seismic wave speed anomalies in the mantle that are commonly interpreted as slabs, the remnants of subducted lithosphere. However, classical travel-time tomography relies on the inversion of travel times of a few easily identifiable body wave phases along ray paths or volumetric sensitivity kernels, which is strongly dependent on the geometry of seismic sources and receivers. Since both of these are primarily clustered on modern convergent plate boundaries, the resulting tomographic resolution is highly variable across the mantle. Full-waveform inversion (FWI) attempts to reduce this dependence by fitting whole seismograms, thereby including many reflected and refracted body wave phases to enhance the volumetric sensitivity of the inversion.Here, we analyse a new global tomographic model constructed using FWI. The mantle structure imaged in this model reveals significantly more positive seismic wave speed anomalies in the mantle when compared to travel-time tomography, particularly in regions with low seismic activity and limited station coverage. Notably, FWI detects positive wave speed anomalies with slab-like morphologies at ~1000 km depth beneath the Pacific Ocean that fall outside the coverage of classical travel-time tomography. We demonstrate the sensitivity of FWI to wave speed anomalies below the western Pacific using forward wavefield modelling. Importantly, we find that these newly imaged positive wave speed anomalies do not correspond to reconstructed subduction zones in existing global plate reconstructions.Our work challenges the widespread assumption that positive wave speed anomalies (exclusively) represent subducted slabs, highlighting potential gaps in either global plate reconstructions or the current understanding of the nature of seismic anomalies in the mantle.

Distributed Acoustic Sensing (DAS) technology is currently used to monitor seismic activity, offering a unique spatially-dense representation of the along-the-cable strain wavefield. Traditional seismic networks typically rely on the timing of specific seismic phases to estimate source locations. In this context, DAS arrays may fail to provide accurate traveltimes because of spatially-heterogeneous waveforms. The motivations are (but not limited to) the directional sensitivity, the heterogeneous cable ground-coupling and the enhanced sensitivity to lateral variations in the medium elastic properties. The resulting fluctuations in signal-to-noise ratios of the dense DAS channels pose significant challenges in the automatic picking of body phases, e.g., P-wave Absolute Arrival Times (P-ARTs). Consequently, the complex distribution of the estimated traveltimes impacts the accuracy of event locations, especially if incorrect assumptions on error statistics (e.g., Normal distribution) are considered. In this study, we address this issue by exploiting the intrinsic DAS measurements' spatial density and testing selected P-wave Differential Arrival Times (P-DATs) for source location. We estimate P-DATs for all the possible DAS channel pairs by identifying the time delay corresponding to the peak of each cross-correlation function. Subsequently, we select P-DATs based on two criteria: interchannel distance and cross-correlation index value. This procedure is often employed to reduce the risk of mixing delay times from coherent and incoherent waveforms. As a first test, using a probabilistic inversion (Hamiltonian Monte Carlo method), we demonstrate how the selected P-DATs provide a better constraint on the event's azimuthal direction compared to P-ARTs. Then, as a second experiment, we move from a subjective selection of P-DATs. To do so, we test a fully-automated and data-driven covariance matrix weighting procedure, in a probabilistic inversion scheme. Specifically, we compute posterior probability distributions for both the physical parameters (event location) and hyperparameters related to data features (interchannel distance and cross-correlation index thresholds). In this scheme, the hyperparameters define each weight along the diagonal of the covariance matrix. These tests offer useful insights into the utilization of P-DATs for event location with DAS. Moreover, we provide an automatic approach to avoid subjective biases based on pre‐conceptions in the a-priori data selection.

Mitigating the energy requirements of artificial intelligence requires novel physical substrates for computation. Phononic metamaterials have vanishingly low power dissipation and hence are a prime candidate for green, always-on computers. However, their use in machine learning applications has not been explored due to the complexity of their design process. Current phononic metamaterials are restricted to simple geometries (e.g., periodic and tapered) and hence do not possess sufficient expressivity to encode machine learning tasks. A non-periodic phononic metamaterial, directly from data samples, that can distinguish between pairs of spoken words in the presence of a simple readout nonlinearity is designed and fabricated, hence demonstrating that phononic metamaterials are a viable avenue towards zero-power smart devices. Elastic neural networks composed of phononic metamaterials respond differently to different spoken commands, passively solving a speech classification problem. Their design harnesses the vanishingly low power dissipation of elastic waves, combined with the high expressivity and efficient simulation of metamaterials. This capability can be leveraged to build smart sensors that detect events without standby power consumption.image

Determining Earth’s structure is paramount to unravel its interior dynamics. Seismic tomography reveals positive wave speed anomalies throughout the mantle that spatially correlate with the expected locations of subducted slabs. This correlation has been widely applied in plate reconstructions and geodynamic modelling. However, global travel-time tomography typically incorporates only a limited number of easily identifiable body wave phases and is therefore strongly dependent on the source-receiver geometry. Here, we show how global full-waveform inversion is less sensitive to source-receiver geometry and reveals numerous previously undetected positive wave speed anomalies in the lower mantle. Many of these previously undetected anomalies are situated below major oceans and continental interiors, with no geologic record of subduction, such as beneath the western Pacific Ocean. Moreover, we find no statistically significant correlation positive anomalies as imaged using full-waveform inversion and past subduction. These findings suggest more diverse origins for these anomalies in Earth’s lower mantle, unlocking full-waveform inversion as an indispensable tool for mantle exploration.

The use of the probabilistic approach to solve inverse problems is becoming more popular in the geophysical community, thanks to its ability to address nonlinear forward problems and to provide uncertainty quantification. However, such strategy is often tailored to specific applications and therefore there is a lack of a common platform for solving a range of different geophysical inverse problems and showing potential and pitfalls. We demonstrate a common framework to solve such inverse problems ranging from, e.g, earthquake source location to potential field data inversion and seismic tomography. Within this approach, we can provide probabilities related to certain properties or structures of the subsurface. Thanks to its ability to address high-dimensional problems, the Hamiltonian Monte Carlo (HMC) algorithm has emerged as the state-of-the-art tool for solving geophysical inverse problems within the probabilistic framework. HMC requires the computation of gradients, which can be obtained by adjoint methods, making the solution of tomographic problems ultimately feasible. These results can be obtained with "HMCLab", a tool for solving a range of different geophysical inverse problems using sampling methods, focusing in particular on the HMC algorithm. HMCLab consists of a set of samplers and a set of geophysical forward problems. For each problem its misfit function and gradient computation are provided and, in addition, a set of prior models can be combined to inject additional information into the inverse problem. This allows users to experiment with probabilistic inverse problems and also address real-world studies. We show how to solve a selected set of problems within this framework using variants of the HMC algorithm and analyze the results. HMCLab is provided as an open source package written both in Python and Julia, welcoming contributions from the community.

SUMMARYIce streams are major contributors to ice sheet mass loss and sea level rise. Effects of their dynamic behaviour are imprinted into seismic properties, such as wave speeds and anisotropy. Here, we present results from a distributed acoustic sensing (DAS) experiment in a deep ice-core borehole in the onset region of the Northeast Greenland Ice Stream, with focus on phenomenological and methodological aspects. A series of active seismic surface sources produced clear recordings of the P and S wavefield, including internal reflections, along a 1500 m long fibre-optic cable that was placed into the borehole. The combination of nonlinear traveltime tomography with a firn model constrained by multimode surface wave data, allows us to invert for P and S wave speeds with depth-dependent uncertainties on the order of only 10 m s−1, and vertical resolution of 20–70 m. The wave speed model in conjunction with the regularly spaced DAS data enable a straightforward separation of internal upward reflections followed by a reverse-time migration that provides a detailed reflectivity image of the ice. While the differences between P and S wave speeds hint at anisotropy related to crystal orientation fabric, the reflectivity image seems to carry a pronounced climatic imprint caused by rapid variations in grain size. Further improvements in resolution do not seem to be limited by the DAS channel spacing. Instead, the maximum frequency of body waves below ∼200 Hz, low signal-to-noise ratio caused by poor coupling, and systematic errors produced by the ray approximation, appear to be the leading-order issues. Among these, only the latter has a simple existing solution in the form of full-waveform inversion. Improving signal bandwidth and quality, however, will likely require a significantly larger effort in terms of both sensing equipment and logistics.

: The continuously increasing quantity and quality of seismic waveform data carry the potential to provide images of the Earth’s internal structure with unprecedented detail. Harnessing this rapidly growing wealth of information, however, constitutes a formidable challenge. While the emergence of faster supercomputers helps to accelerate existing algorithms, the daunting scaling properties of seismic inverse problems still demand the development of more efficient solutions. The diversity of seismic inverse problems – in terms of scientific scope, spatial scale, nature of the data, and available resources – precludes the existence of a silver bullet. Instead, efficiency derives from problem adaptation. Within this context, this chapter describes a collection of methods that are smart in the sense of exploiting specific properties of seismic inverse problems, thereby increasing computational efficiency and usable data volumes, sometimes by orders of magnitude. These methods improve different aspects of a seismic inverse problem, for instance, by harnessing data redundancies, adapting numerical simulation meshes to prior knowledge of wavefield geometry, or permitting long-distance moves through model space for Monte Carlo sampling.

We present a distributed acoustic sensing (DAS) experiment at Grímsvötn, Iceland. This is intended to investigate volcano-microseismicity at Grímsvötn specifically, and to assess the suitability of DAS as a subglacial volcano monitoring tool in general. In spring 2021, we trenched a 12 km long fiber-optic cable into the ice sheet around and within the caldera, followed by nearly one month of continuous recording. An image processing algorithm that exploits spatial coherence in DAS data detects on average ~100 events per day, almost 2 orders of magnitude more than in the regional earthquake catalog. A nonlinear Bayesian inversion reveals the presence of pronounced seismicity clusters, containing events with magnitudes between −3.4 and 1.7. Their close proximity to surface volcanic features suggests a geothermal origin. In addition to painting a fine-scale picture of seismic activity at Grímsvötn, this work confirms the potential of DAS in subglacial volcano monitoring.

The M series of chips produced by Apple have proven a capable and power-efficient alternative to mainstream Intel and AMD x86 processors for everyday tasks. Additionally, the unified design integrating the central processing and graphics processing unit, have allowed these M series chips to excel at many tasks with heavy graphical requirements without the need for a discrete graphical processing unit (GPU), and in some cases even outperforming discrete GPUs. In this work, we show how the M series chips can be leveraged using the Metal Shading Language (MSL) to accelerate typical array operations in C++. More importantly, we show how the usage of MSL avoids the typical complexity of CUDA or OpenACC memory management, by allowing the central processing unit (CPU) and GPU to work in unified memory. We demonstrate how performant the M series chips are on standard one-dimensional and two-dimensional array operations such as array addition, SAXPY and finite difference stencils, with respect to serial and OpenMP accelerated CPU code. The reduced complexity of implementing MSL also allows us to accelerate an existing elastic wave equation solver (originally based on OpenMP accelerated C++) using MSL, with minimal effort, while retaining all CPU and OpenMP functionality. The resulting performance gain of simulating the wave equation is near an order of magnitude for specific settings. This gain attained from using MSL is similar to other GPU-accelerated wave-propagation codes with respect to their CPU variants, but does not come at much increased programming complexity that prohibits the typical scientific programmer to leverage these accelerators. This result shows how unified processing units can be a valuable tool to seismologists and computational scientists in general, lowering the bar to writing performant codes that leverage modern GPUs.

We present the results of an experiment with Distributed Acoustic Sensing (DAS) on Grímsvötn in Iceland. DAS is a novel detection method that samples the strain wavefield due to ground motion along a fibre-optic cable with high temporal (kHz) and spatial (m) resolution. Consequently, it has the potential to increase our understanding of physical volcanic processes. We deployed a 12 km long fibre-optic cable for one month (May 2021) on Grímsvötn, Iceland’s most active volcano, which is completely covered by the large Vatnajökull ice sheet. The cable was trenched 50 cm into the ice, following the caldera rim and ending near the central point of the caldera on top of a subglacial lake. A large number of hammer blow experiments allow us to estimate the Rayleigh wave dispersion curves, and thickness of the ice layer on top of the volcanic rock. We have discovered previously undetected levels of seismicity, with up to several hundreds of local events per day, using an automated earthquake detection algorithm that is based on image processing techniques. First arrival picks are identified with an automated cross-correlation based algorithm, developed specifically for complex and local events recorded with DAS. The first arrival times, combined with a probabilistic interpretation and the Hamiltonian Monte Carlo algorithm, allow us to estimate event locations and their respective uncertainties, even in the absence of a detailed velocity model. The detection and localisation of the recorded events paints a differentiated picture of Grímsvötn’s volcano-seismicity. The preliminary results of our experiment highlight the potential of DAS for studies of active volcanoes covered by glaciers, and we hope that this research will contribute to the fields of volcano monitoring and hazard assessment.

We present a probabilistic approach to constrain the density distribution in the Earth based on surface wave dispersion. Despite its outstanding importance in studies of the Earth’s thermo-chemical state and dynamics, 3D density variations remain poorly known, thereby posing one of the major challenges in geophysics. Since the sensitivity of most seismic data to density is small compared to sensitivity with respect to seismic velocities, regularisation in traditional deterministic inversion tends to bias the recovered density image significantly. To avoid this issue, we propose to solve a regularisation-free Bayesian inference problem using the Hamiltonian Monte Carlo Markov Chain algorithm. In the interest of simplicity, we consider anisotropic stratified media, where dispersion curves and corresponding sensitivity kernels can be computed semi-analytically. Exploiting derivative information for efficient sampling, Hamiltonian Monte Carlo approximates the posterior probability density of all model parameters, namely the P-wave velocities vPV and vPH , the S-wave velocities vSV and vSH , the anisotropy parameter η, and, of course, density ρ. The proposed method forms the foundation of an open-source tool box that can be used to assess the unbiased ability of surface wave dispersion data, characterised in terms of frequency and modal content, to constrain density variations and their trade-offs with other Earth model parameters.

We present `psvWave', a basic numerical finite difference solver for Python and C++, specifically targeted at seismologists. The solver is based on the well-established staggered grid approaches developed for the P-SV elastic wave equation. Although its functionality is limited (solely moment tensor sources, only Ricker wavelets source time functions), it does possess the ability to perform adjoint simulations, and its performance has so far allowed the development of Bayesian sampling for Full-Waveform Inversion using the Hamiltonian Monte Carlo algorithm. We present this as an open source project, and invite anyone to contribute.

The Hamiltonian Monte Carlo method (HMC) is gaining popularity in the geophysical community to fully address nonlinear inverse problems and related uncertainty quantification. We present here an application of HMC to invert seismic data in the acoustic approximation in the context of reflection seismology. We address a 2-D problem, in the form of a vertical cross section where both source and receivers are located near the surface of the model. To solve the forward problem we utilise the finite-difference method with PML absorbing boundary conditions. The observed data are represented by a set of shotgathers. The crucial aspect for a successful application of the HMC lies in the capability of performing gradient computations in an efficient manner. To this end, we use the adjont state method to compute the gradient of the misfit functional, which has a computational cost of only about twice that of the forward computation, a very efficient strategy. From the collection of samples characterising the posterior distribution obtained with the HMC, we can derive quantities of interest using statistical analysis and assess uncertainties. We illustrate an application of this methodology on a synthetic test mimicking the setup encountered in exploration problems.

Uncertainty quantification is an essential part of many studies in Earth science. It allows us, for example, to assess the quality of tomographic reconstructions, quantify hypotheses and make physics-based risk assessments. In recent years there has been a surge in applications of uncertainty quantification in seismological inverse problems. This is mainly due to increasing computational power and the ‘discovery’ of optimal use cases for many algorithms (e.g., gradient-based Markov Chain Monte Carlo (MCMC). Performing Bayesian inference using these methods allows seismologists to perform advanced uncertainty quantification. However, oftentimes, Bayesian inference is still prohibitively expensive due to large parameter spaces and computationally expensive physics. Simultaneously, machine learning has found its way into parameter estimation in geosciences. Recent works show that machine learning both allows one to accelerate repetitive inferences [e.g. Shahraeeni & Curtis 2011, Cao et al. 2020] as well as speed up single-instance Monte Carlo algorithms using surrogate networks [Aleardi 2020]. These advances allow seismologists to use machine learning as a tool to bring accurate inference on the subsurface to scale. In this work, we propose the novel inclusion of adjoint modelling in machine learning accelerated inverse problems. The aforementioned references train machine learning models on observations of the misfit function. This is done with the aim of creating surrogate but accelerated models for the misfit computations, which in turn allows one to compute this function and its gradients much faster. This approach ignores that many physical models have an adjoint state, allowing one to compute gradients using only one additional simulation. The inclusion of this information within gradient-based sampling creates performance gains in both training the surrogate and the sampling of the true posterior. We show how machine learning models that approximate misfits and gradients specifically trained using adjoint methods accelerate various types of inversions and bring Bayesian inference to scale. Practically, the proposed method simply allows us to utilize information from previous MCMC samples in the algorithm proposal step. The application of the proposed machinery is in settings where models are extensively and repetitively run. Markov chain Monte Carlo algorithms, which may require millions of evaluations of the forward modelling equations, can be accelerated by off-loading these simulations to neural nets. This approach is also promising for tomographic monitoring, where experiments are repeatedly performed. Lastly, the efficiently trained neural nets can be used to learn a likelihood for a given dataset, to which subsequently different priors can be efficiently applied. We show examples of all these use cases. Lars Gebraad, Christian Boehm and Andreas Fichtner, 2020: Bayesian Elastic Full‐Waveform Inversion Using Hamiltonian Monte Carlo. Ruikun Cao, Stephanie Earp, Sjoerd A. L. de Ridder, Andrew Curtis, and Erica Galetti, 2020: Near-real-time near-surface 3D seismic velocity and uncertainty models by wavefield gradiometry and neural network inversion of ambient seismic noise. Mohammad S. Shahraeeni and Andrew Curtis, 2011: Fast probabilistic nonlinear petrophysical inversion. Mattia Aleardi, 2020: Combining discrete cosine transform and convolutional neural networks to speed up the Hamiltonian Monte Carlo inversion of pre‐stack seismic data.

Constraints on the 3-D density structure of Earth’s mantle provide important insights into the nature of seismically observed features, such as the Large Low Shear Velocity Provinces (LLSVPs) in the lower mantle under Africa and the Pacific. The only seismic data directly sensitive to density variations throughout the entire mantle are normal modes: whole Earth oscillations that are induced by large earthquakes (Mw > 7.5). However, their sensitivity to density is weak compared to the sensitivity to velocity and different studies have presented conflicting density models of the lower mantle. For example, Ishii & Tromp (1999) and Trampert et al. (2004) have found that the LLSVPs have a larger density than the surrounding mantle, while Koelemeijer et al. (2017) used additional Stoneley-mode observations, which are particularly sensitive to the core-mantle boundary region, to show that the LLSVPs have a lower density. Recently, Lau et al. (2017) have used tidal tomography to show that Earth's body tides prefer dense LLSVPs.A large number of new normal-mode splitting function measurements has become available since the last density models of the entire mantle were published. Here, we show the models from our inversion of these recent data and compare our results to previous studies. We find areas of high as well as low density at the base of the LLSVPs and we find that inside the LLSVPs density varies on a smaller scale than velocity, indicating the presence of compositionally distinct material. In fact, we find low correlations between the density and velocity structure throughout the entire mantle, revealing that compositional variations are required at all depths inside the mantle.

We propose methods to efficiently explore the generalized nullspace of (non-linear) inverse problems, defined as the set of plausible models that explain observations within some misfit tolerance. Owing to the random nature of observational errors, the generalized nullspace is an inherently probabilistic entity, described by a joint probability density of tolerance values and model parameters. Our exploration methods rest on the construction of artificial Hamiltonian systems, where models are treated as high-dimensional particles moving along a trajectory through model space. In the special case where the distribution of misfit tolerances is Gaussian, the methods are identical to standard Hamiltonian Monte Carlo, revealing that its apparently meaningless momentum variable plays the intuitive role of a directional tolerance. Its direction points from the current towards a new acceptable model, and its magnitude is the corresponding misfit increase. We address the fundamental problem of producing independent plausible models within a high-dimensional generalized nullspace by autotuning the mass matrix of the Hamiltonian system. The approach rests on a factorized and sequentially preconditioned version of the L-BFGS method, which produces local Hessian approximations for use as a near-optimal mass matrix. An adaptive time stepping algorithm for the numerical solution of Hamilton's equations ensures both stability and reasonable acceptance rates of the generalized nullspace sampler. In addition to the basic method, we propose variations of it, where autotuning focuses either on the diagonal elements of the mass matrix or on the macroscopic (long-range) properties of the generalized nullspace distribution. We quantify the performance of our methods in a series of numerical experiments, involving analytical, high-dimensional, multimodal test functions. These are designed to mimic realistic inverse problems, where sensitivity to different model parameters varies widely, and where parameters tend to be correlated. The tests indicate that the effective sample size may increase by orders of magnitude when autotuning is used. Finally, we present a proof of principle of generalized nullspace exploration in viscoelastic full-waveform inversion. In this context, we demonstrate (1) the quantification of inter- and intraparameter trade-offs, (2) the flexibility to change model parametrization a posteriori, for instance, to adapt averaging length scales, (3) the ability to perform dehomogenization to retrieve plausible subwavelength models and (4) the extraction of a manageable number of alternative models, potentially located in distinct local minima of the misfit functional.

Mitigating the energy requirements of artificial intelligence requires novel physical substrates for computation. Phononic metamaterials have a vanishingly low power dissipation and hence are a prime candidate for green, always-on computers. However, their use in machine learning applications has not been explored due to the complexity of their design process: Current phononic metamaterials are restricted to simple geometries (e.g. periodic, tapered), and hence do not possess sufficient expressivity to encode machine learning tasks. We design and fabricate a non-periodic phononic metamaterial, directly from data samples, that can distinguish between pairs of spoken words in the presence of a simple readout nonlinearity; hence demonstrating that phononic metamaterials are a viable avenue towards zero-power smart devices.