Full waveform inversion (FWI) is a powerful technique for estimating high-resolution subsurface velocity models by minimizing the discrepancy between modelled and observed seismic data. However, the oscillatory nature of seismic waveforms makes point-wise discrepancy measures highly prone to cycle skipping, especially when the initial velocity model is inadequate. To address this challenge, various alternative misfit functions have been proposed in the literature, each with unique strengths and limitations. Dynamic time warping (DTW) is a popular technique in signal processing for aligning time series using dynamic programming. While a differentiable variant of DTW has been recently proposed, its use in FWI is hindered by high-frequency artefacts in the adjoint source and the substantial computational cost of gradient evaluations. In this study, we propose a neural network-based approach to learn the time shifts that align two time series in a supervised manner. The trained network is then utilized to compare traces from observed and modelled seismic data, offering a stable and computationally efficient alternative to DTW. Furthermore, the inherent differentiability of neural networks via backpropagation enables seamless integration into the FWI framework as a misfit function. We validate this approach on two synthetic datasets, namely the Marmousi model and the Chevron blind test dataset, demonstrating in both cases a similar convergence behaviour to that of Soft-DWT whilst drastically reducing the computational time of the adjoint source calculation.

Determining the dynamics of fluid-rock interactions is key to deciphering processes in the lithosphere and their importance in industrial applications, including the storage of CO2 and hydrogen, as well as in the development of geothermal energy technologies. Many geological occurrences show that when fluids interact with rocks, they can form their own pathways by developing a temporary network of connected pores as a result of mineral replacement reactions. Nonetheless, following the reaction's cessation, most pores become isolated retarding fluid flow. Here, we investigate transient porosity generation phenomena by conducting 4D (3D plus time) synchrotron tomography experiments using a salt analogue system. To capture the rapid evolution of reaction-induced pore networks, we acquire a full tomography volume every minute. After segmenting this extensive dataset using a deep convolutional neural network, the dynamics of the reaction are quantified by spatiotemporal correlation functions. These high-order statistics enable us to evaluate the evolution of the structural and morphological properties of the induced pore network during the salt replacement reaction. Inspired by image editing techniques in computer vision, we then train a leading-edge generative model, StyleGAN2-ADA (Karras et al., 2020[P(1] ), utilizing the 4D tomography dataset. Our results demonstrate that these advanced generative models can accurately simulate the microstructural evolution observed in the experiment. Next, we delve into the potential of deploying the trained model to reconnect isolated pores observed in data sets of natural mineral replacement reactions. We finally apply a voxel-based finite element method to simulate fluid flow in the salt and natural system. Our deep-learning method not only pioneers in determining transient fluid-rock interaction phenomena, but also, for the first time, enables a direct estimation of reaction-induced permeabilities.Putnis, A. Mineral replacement reactions. Rev. mineralogy geochemistry 70, 87–124 (2009).Karras, T. et al. Training generative adversarial networks with limited data. Adv. Neural Inf. Process. Syst. 33, 12104–12114 (2020).

Laboratory stick-slip experiments are a simple analogue for the earthquake cycle. The acoustic emissions (AE) of these experiments have been shown to contain hidden patterns. Machine Learning (ML) can extract these patterns and information on the fault state can be inferred (e.g. shear stress and time to failure). Two different ML approaches have been used in the past: 1) ensemble tree models, which are relatively easy to evaluate why they made a certain prediction, but only look at a snapshot in time and 2) deep neural networks using Long Short-Term Memory (LSTM), which have the ability to find patterns in the temporal changes in the signal, but act more as a black-box model, so the final predictions are hard to evaluate. Here we introduce an additional step in the workflow that can be used to allow the ensemble tree models information about the temporal changes of the input features. Furthermore, it is able to quantify and visualize whether a pattern is repetitive or not. Like earlier studies we start by calculating (statistical) features using a rolling window on the AE. The features are not directly used as the input of the model, but are placed in a larger Hankel matrix, where the consecutive time windows are the rows of the matrix. Using Principal Component Analysis (PCA) and Uniform Manifold Approximation and Projection (UMAP) we create an embedded version of this array that holds temporal information of features calculated in the previous step. Visual inspection of these embeddings shows that some features map to very distinct patterns that are repetitive over the majority of the stick-slip cycles. The advantage of this method is that an inverse mapping is easily available, allowing for an interpretable embedding of the data.

SUMMARY In this work we present a novel, experimentally efficient set-up for performing non-contacting laser vibrometry on geologic materials and their analogues. We show it is possible to acoustically monitor a granular material experiment in real time compared to the typical timescale of analogue modelling experiments. We acquire non-contacting waveform data with consistently high signal-to-noise ratio. Compared to previously used standard contacting transducers, the novel joint use of sources and receivers that are both laser-based resulted in measured signals with improved waveforms and temporal bandwidths. These data acquisition improvements, in our case where surface waves are prominent in the data, enable enhanced multichannel surface wave processing, for example, in terms of reliable dispersion curve estimates. We find, given the high waveform fidelity of our acquisition system, that the observed surface waves are highly sensitive to relatively small changes in the medium’s elastic properties, making them a demonstrably reliable to monitor any processes that affect elasticity in these models in near real time. As a demonstration, we continuously monitor a scaled analogue model containing granular glass beads. By continuously monitoring—that is, performing repeatable active-source acousto-seismic surveys at short time-lapse intervals—over a period of 10 hr, we find that an increase of relative humidity of 10 per cent can lead to as much as a factor of two increase in surface wave group velocities. Finally, we discuss future applications of the developed method by considering surface wave inversion for fault and stress monitoring during the deformation of a model.

Summary Point-estimate statistics for inverse problems, such as the maximum a-posteriori estimator, require the solution of a problem that is typically ill-conditioned for seismic imaging applications, with important implications in terms of computational complexity. Ill-conditioning is arguably an even bigger challenge for uncertainty quantification, which aims at a more comprehensive characterisation of the posterior distribution. One classical strategy to assuage these issues is the multiscale approach, where the original problem is broken down into a hierarchical sequence of sub-problems with increasing computational complexity. Each sub-problem describes scale-dependent features of the overall solution for any chosen scale-dependent decomposition. This gives rise to an efficient iterative method that progressively builds a solution from "coarse" to "fine" scales. We propose to leverage recent developments in machine-learning-based variational inference for uncertainty quantification that uses a wavelet-based generative model of the posterior distribution. The architectural design is based on a normalising flow that generates the scales of a sample sequentially and conditionally based on the coarser scales. As a first application of this framework, we study post-stack seismic inversion here.

Summary A fundamental aspect of modern seismic data processing involves the reconstruction of wavefields in areas where missing sources or receivers result in data gaps. Despite recent developments in data-driven multidimensional processing, limitations persist due to incomplete sampling in source, receiver, or both coordinates. Similarly, dense acquisition systems, though costly, are confronted with physical and environmental limitations, hindering the recording of well-sampled data. A recent innovation leveraging low-rank completion techniques offers superior seismic data restoration, potentially enabling data kernel compression and unlocking previously restricted modern processing methods in 3D field scenarios. By introducing a fast column/row-wise cyclic shearing of the data array as stored in the regular acquisition domain, we explore an alternative low-dimensional domain for seismic reconstruction that is designed to align the prominent energy contribution distributed along diagonal entries while preserving the inherent local proximity among source-receiver pairs. This approach avoids the implicit zero-padding in conventional midpoint-offset domain matrix completion schemes. We show the simultaneous reconstruction of missing sources and receivers in synthetic and field data by accommodating additional lateral constraints by means of regularization.

Constraining strain localization and the growth of shear fabrics within brittle fault zones at sub-seismic slip rates is important for understanding fault strength and frictional stability. We conducted direct shear experiments on simulated sandstone-derived fault gouges at an effective normal stress of 40 MPa, a pore pressure of 15 MPa, and a temperature of 100 degrees C. Using a passive strain marker and X-ray Computed Tomography, we analyzed the spatial distribution of deformation in gouges deformed in the strain-hardening, subsequent strain-softening, and then steady-state regimes at displacement rates of 1, 30, and 1,000 mu m/s. We developed a machine-learning-based automatic boundary detection method to recognize the shear fabrics and quantify displacement partitioning between each fabric element. Our results show fabrics oriented along R1 and Y (including boundary) shears are the two major fabric elements. At rates of 1 and 30 mu m/s, the relative amount of displacement on R1 shears is displacement dependent, increasing to similar to 20% of the total displacement up to the strain-softening stage, then decreasing to similar to 10%-18% at the steady state. This trend is absent at the high rate where similar to 18% of the displacement occurs on R1 shears throughout all investigated stages. At all rates, the relative amount of displacement on Y shears increases linearly with displacement to a total of larger than 50% at the steady state. Our study provides constraints on the development of the active slip zone, which is an important factor controlling heating and weakening associated with small-magnitude earthquakes with limited displacement (mm-dm), such as induced seismicity. In the past few decades, several studies have focused on the mechanical behavior of simulated fault zones to understand earthquake nucleation. However, minor attention has been paid to the microstructural characterization of fault-zone gouges due to the difficulties in quantifying deformation. An understanding of brittle fault-zone fabrics and their development provides crucial constraints on the mechanical strength and stability of faults. We conducted laboratory experiments on simulated sandstone-derived fault gouges with a passive strain marker under the conditions relevant to earthquake nucleation to explore the development and evolution of shear zone fabrics and their relations to fault strength. We combined X-ray Computed Tomography (XCT) and a custom-designed machine-learning-based automatic boundary detection method to analyze the spatial distribution of gouge deformation and to quantify displacement partitioning between deformation features. The results show that our samples have similar evolution of the shear fabrics, partitioning of displacement, and mechanical response with increasing shear strain at all tested nucleation velocities. The evolution of the mechanical behavior from strain-hardening, to softening, to steady-state stages is related to the transformation of R1 to shear-parallel shear bands. Up to 50% of the total displacement can be accommodated within shear-parallel shear bands, facilitating strain weakening of the materials. We performed 3-D analyses on the evolution of shear zone fabrics within sandstone-derived fault gouges utilizing the X-ray CT technique Our samples show similar evolution of shear fabrics, slip partitioning, and mechanical response with shear strain at all tested velocities Up to 50% of the total imposed displacement can be accommodated within shear-parallel shear bands during earthquake nucleation

Summary Multi-dimensional deconvolution (MDD) is a highly desired inversion method that can in principle address multiple challenges in a variety of seismic applications ranging from ocean-bottom data processing, to advanced target-oriented imaging and monitoring. Recent studies show that time-domain MDD, with efficient operator schemes and physics-based constraints, can significantly outperform its more commonly used frequency-domain counterpart in terms of robustness and waveform fidelity when retrieving target responses. Building on this time-domain framework, we propose an alternative MDD scheme based on inversion of point-spread-functions (PSFs), i.e., by recasting the original problem from up/down multidimensional deconvolution into a deblurring problem in the time domain. Due to the symmetry properties of the PSFs in deblurring-based MDD, we show that a relatively simple approach to restricting numerical integration yields inversion results with relatively little loss of accuracy for operators with notably reduced dimensions, which is not the case for the conventional MDD scheme. We illustrate our approach and findings with a numerical example of depth-domain redatuming for target-oriented imaging in a complex subsalt setting. Our approach and results show promise for enabling large-scale time-domain MDD calculations, as well as for bringing more complex MDD applications, such as elastic MDD, into industrial practice.

In complex geologic settings and in the presence of sparse acquisition systems, seismic migration images manifest as nonstationary blurred versions of the unknown subsurface model. Thus, image-domain deblurring is an important step to produce interpretable and high-resolution models of the subsurface. Most deblurring methods focus on inverting seismic images for their underlying reflectivity by iterative least-squares inversion of a local Hessian approximation; this is obtained by either direct modeling of the so-called point-spread functions (PSFs) or by a migration-demigration process. In this work, we adopt a novel deep-learning (DL) framework, based on invertible recurrent inference machines (i-RIMs), which allows approaching any inverse problem as a supervised learning task informed by the known modeling operator (convolution with PSFs in our case): our algorithm can directly invert migrated images for impedance perturbation models, assisted with the prior information of a smooth velocity model and the modeling operator. Because i-RIMs are constrained by the forward operator, they implicitly learn to shape/regularize output models in a training-data-driven fashion. As such, the resulting deblurred images indicate great robustness to noise in the data and spectral deficiencies (e.g., due to limited acquisition). The key role played by the i-RIM network design and the inclusion of the forward operator in the training process is supported by several synthetic examples. Finally, using field data, we find that i-RIM-based deblurring has great potential in yielding robust, high-quality relative impedance estimates from migrated seismic images. Our approach could be of importance toward future DL-based quantitative reservoir characterization and monitoring.

It has previously been shown that features calculated using rolling windows of acoustic emissions (AE) emitted from laboratory shear experiments give a robust basis to predict shear stress and time to failure (TTF) using machine learning (ML). Early studies used ensemble tree methods like random forests and gradient boosted trees for learning and inference. The upside of these methods is that they allow for relatively easy analysis of the behaviour of learned models. It is possible to go through every decision inside the trees to evaluate why the model made a certain prediction, based on the features fed into the model. These predictions, however, are only based on local time-window information. The recent AE inside the time window hold no information of the temporal evolution of the features. More recently, progress has been made in improving the quality of predictions by making use of for example Deep Neural Networks (DNN’s) containing techniques like long short-term memory (LSTM) which have the advantage of learning from temporal evolutions. These more advanced ML techniques act however, more like a black-box model and as a result it is more complicated, to nearly impossible, to explain why a certain prediction is made. Here dimension reduction techniques like Principal Component Analysis (PCA) and Uniform Manifold Approximation and Projection (UMAP) are used to visualize the temporal evolution of features. The embeddings of features can both be used for visual inspection and as the input for ML models. The embedded features are fed to ensemble learning methods like a random forest as additional input to predict shear stress and TTF. This technique provides the relatively simple ML method with temporal information, while still enabling the interpretation of predictions made by the model. This approach is applied to both earlier studied data of the double-direct shear apparatus located at the Pennsylvania State University and new data obtained by experiments in the Rotary Shear Apparatus located at the Utrecht University High Pressure and Temperature Laboratory. Visual inspection of the embedded features shows that different shear stress measurements map to distinct patterns, meaning that it is not only the local time windows of the feature that hold predictive information, but also their global temporal evolution, and that the patterns repeat for multiple laboratory earthquake cycles. Feeding the embedded features of variance and kurtosis (the features previously shown to be most influential on the final prediction) into a ML model additionally to the traditional time windows results in highly improved predictions (R2-scores in the order of 0.9), nearly as good as the state-of-the-art LSTM models. From this it can be concluded that the temporal evolution of features must be considered, because it holds additional information critical to making the best predictions possible. Besides improving the predictions, the embedded features also provide a quantitative data-driven framework to study the temporal evolution of features.

We present and demonstrate our new application of a geophysical seismic technique to acoustically characterise and image layers with different impedance contrast in analogue models. A high-powered pulsed laser in combination with a mirror galvanometer is used to generate a powerful acoustic shockwave at any point of the surface of the analogue model. Reflections, refractions, and diffractions of the acoustic source wave, induced by internal structures inside an analogue model, produce vibrations of the top surface of a model, which are measured by laser vibrometer.Using our setup, we acquire seismic receiver gathers in less than a minute. Interpretation of the gathers allowed to identify the presence of internal reflecting and refracting material interfaces. In a series of test models, we determined the speed of both P-waves and surface waves in a multitude of brittle analogue materials. In uniform layered models we performed 1D inversion using the gathered waveform data. The results are validated by simulating the test experiments in a finite-difference solver. The novel method will be developed further, aiming to determine stress build-up in the material prior to fault formation or activity.

Marchenko-type integrals typically relate so-called focusing functions and Green’s functions via the reflection response measured on the open surface of a volume of interest. Originating from one dimensional inverse scattering theory, the extension to two and three dimensions set in motion various new developments regarding imaging in complex materials. This extension, however, is based on wavefield decomposition inside the volume and a truncated medium state, i.e. a version of the medium that is reflection-free underneath the focusing location, suggesting that evanescent, refracted and diving waves cannot be included in the representation. We elaborate on a new derivation for Marchenko-like integrals that (i) extends the concept of wavefield focusing by using a generalised homogeneous Green’s function, (ii) is based on partial differential equations and thereby allows for additional insights and a new physical intuition for Marchenko equations, (iii) unifies wavefield focusing for open and closed boundary systems, (iv) does not require wavefield decomposition or a truncated medium state, thus including the full wavefield Green’s function, (v) enables using forward modelling to obtain, e.g., Marchenko-type, time-compact focusing functions. We place a particular focus on the latter point, illustrating and investigating how to solve the underlying partial differential equations for various types of focusing functions. This paves the way for a deeper understanding of focusing functions as well as advanced full wavefield Marchenko schemes. While the derivations are generally presented for the 3D case, we show numerical examples in 1D.

Full waveform inversion is a high-resolution subsurface imaging technique, in which full seismic waveforms are used to infer subsurface physical properties. We present a novel, target-enclosing, full-waveform inversion framework based on an interferometric objective function. This objective function exploits the equivalence between the convolution and correlation representation formulas, using data from a closed boundary around the target area of interest. Because such equivalence is violated when the knowledge of the enclosed medium is incorrect, we propose to minimize the mismatch between the wavefields independently reconstructed by the two representation formulas. The proposed method requires only kinematic knowledge of the subsurface model, specifically the overburden for redatuming, and does not require prior knowledge of the model below the target area. In this sense it is truly local: sensitive only to the medium parameters within the chosen target, with no assumptions about the medium or scattering regime outside the target. We present the theoretical framework and derive the gradient of the new objective function via the adjoint-state method and apply it to a synthetic example with exactly redatumed wavefields. A comparison with FWI of surface data and target-oriented FWI based on the convolution representation theorem only shows the superiority of our method both in terms of the quality of target recovery and reduction in computational cost.

A crucial step in seismic data processing consists in reconstructing the wavefields at spatial locations where faulty or absent sources and/or receivers result in missing data. Several developments in seismic acquisition and interpolation strive to restore signals fragmented by sampling limitations; still, seismic data frequently remain poorly sampled in the source, receiver, or both coordinates. An intrinsic limitation of real-life dense acquisition systems, which are often exceedingly expensive, is that they remain unable to circumvent various physical and environmental obstacles, ultimately hindering a proper recording scheme. In many situations, when the preferred reconstruction method fails to render the actual continuous signals, subsequent imaging studies are negatively affected by sampling artefacts. A recent alternative builds on low-rank completion techniques to deliver superior restoration results on seismic data, paving the way for data kernel compression that can potentially unlock multiple modern processing methods so far prohibited in 3D field scenarios. In this work, we propose a novel transform domain revealing the low-rank character of seismic data that prevents the inherent matrix enlargement introduced when the data are sorted in the midpoint-offset domain and develop a robust extension of the current matrix completion framework to account for lateral physical constraints that ensure a degree of proximity similarity among neighbouring points. Our strategy successfully interpolates missing sources and receivers simultaneously in synthetic and field data.

Marchenko multiple elimination methods remove all orders of overburden-generated internal multiples in a data-driven way. In the presence of thin beds, however, these methods have been shown to underperform. This is because the underlying inverse problem requires the information about short-period internal multiple (SPIM) imprint on the inverse transmission to be correctly constrained. This has been addressed in 1.5D media with energy conservation and minimum-phase reconstruction. Extending the applications to two dimensions and, hence, making the step toward field data were believed to be hampered by the need for a multidimensional minimum-phase reconstruction which is (1) not unique and (2) no algorithm has been found to perform this in practice on band-limited data. Here, we address both of these problems with an approach that includes solving the Marchenko equation with a trivial constraint, evaluating the energy conservation condition of its solutions to find the spatially dependent error syndrome, using the 1.5D minimum-phase reconstruction for each shot gather to find the spatially dependent constraint, and finally using that inside another run of the Marchenko equation solver to find a much-improved result. We find that the method works because in 2D media the expression of SPIMs in the inverse transmission coda is approximately 1.5D. We then investigate a class of models and synthetic data sets to verify where the 1.5D approximation starts breaking down. Our analysis indicates that this approach could perform very well in settings with moderate lateral variations, which also is where the (short-period) internal multiples are most difficult to differentiate from primary reflections.

The key to most subsurface processes is to determine how structural and topological features at small length scales, i.e., the microstructure, control the effective and macroscopic properties of earth materials. Recent progress in imaging technology has enabled us to visualise and characterise microstructures at different length scales and dimensions. However, one limitation of these technologies is the trade-off between resolution and sample size (or representativeness). A promising approach to this problem is image reconstruction which aims to generate statistically equivalent microstructures but at a larger scale and/or additional dimension. In this work, a stochastic method and three generative adversarial networks (GANs), namely deep convolutional GAN (DCGAN), Wasserstein GAN with gradient penalty (WGAN-GP), and StyleGAN2 with adaptive discriminator augmentation (ADA), are used to reconstruct two-dimensional images of two hydrothermally rocks with varying degrees of complexity. For the first time, we evaluate and compare the performance of these methods using multi-point spatial correlation functions—known as statistical microstructural descriptors (SMDs)—ultimately used as external tools to the loss functions. Our findings suggest that a well-trained GAN can reconstruct higher-order, spatially-correlated patterns of complex earth materials, capturing underlying structural and morphological properties. Comparing our results with a stochastic reconstruction method based on a two-point correlation function, we show the importance of coupling training/assessment of GANs with higher-order SMDs, especially in the case of complex microstructures. More importantly, by quantifying original and reconstructed microstructures via different GANs, we highlight the interpretability of these SMDs and show how they can provide valuable insights into the spatial patterns in the synthetic images, allowing us to detect common artefacts and failure cases in training GANs.

SUMMARYFull waveform inversion and least-squares reverse time migration are the leading technologies for imaging with seismic waves. Both of them usually rely (in one way or another) on a single-scattering approximation, i.e. the Born approximation, to compute gradients and obtain an updated model. This approximation linearises the relation between modelled data and model by ignoring multiple scattering. We propose to use the Marchenko integral, an equation originating from inverse scattering theory, to obtain an alternative linear equation. Using the Marchenko method we can retrieve Green’s functions, including all orders of scattering, for virtual sources anywhere within the volume of interest – without prior knowledge of the high-wavelength model variations that induce scattering. Plugging these estimated Green’s functions into the Lippmann–Schwinger integral delivers a Marchenko-linearised relation between the full waveform data and the model. We present this new linearisation strategy and illustrate its advantages and disadvantages by comparing numerical results for different inversion kernels. Our new linearisation is exact, i.e. it does not exclude any orders of scattering, however, it relies on the quality of the Marchenko-derived Green’s functions. These Marchenko-based Green’s functions require an estimate of the first arrivals of the Green’s functions – commonly obtained by modelling in a background medium. Although these first arrival estimates strongly bias our results for inaccurate background models, we find the Marchenko-linearisation to deliver overall slightly better inverted models than the single-scattering approximation.