Abstract Seismic full-waveform inversion (FWI) is a powerful tool for monitoring subsurface changes during carbon capture and storage operations, but its ill-posed nature makes uncertainty quantification (UQ) essential for reliable interpretation. Sampling-based Bayesian methods, such as Markov chain Monte Carlo, provide rigorous UQ but are computationally demanding. In conventional FWI, elastic properties are assigned to densely discretized space-filling cells, resulting in a high-dimensional parameterization that makes large-scale elastic FWI computationally infeasible. To address this challenge, a rock-physics-guided parameter-reduction strategy that compactly represents the CO2 plume geometry using cubic splines controlled by a limited number of nodes was proposed. This parsimonious parameterization not only significantly reduced the number of model parameters and the forward simulations required for effective UQ using sampling methods but also showed potential to improve the practicality and efficiency of other types of UQ methods. Numerical experiments on a cross-well synthetic scenario and a field-scale synthetic case based on the Aquistore storage site in Saskatchewan, Canada, demonstrated that the method efficiently reconstructed the plume shape and its extent and that it converged to consistent posterior distributions across multiple Markov chains.

To maximize the utility of seismic imaging and inversion results, we need to compute not only a final image but also quantify the uncertainty in the image. Although the most thorough approach to quantify the uncertainty is to use a method such as Markov chain Monte Carlo, which systematically samples the entire posterior distribution, this is often inefficient, and not all applications require a full representation of the posterior. We use normalizing flows (NFs), a machine learning technique to perform uncertainty quantification (UQ) in full-waveform inversion (FWI), specifically for time-lapse data. As with any machine learning algorithm, the NF learns only the mapping from the part of the prior spanned by the training data to the distribution of final models spanned by the training data. Here, we make use of this property to perform UQ efficiently by learning a mapping from the prior to the distribution that characterizes the model perturbations within a specific range. Our approach involves using a range of starting models paired with final models from a standard FWI as training data. Although this does not capture the full posterior of the FWI problem, it enables us to quantify the uncertainties associated with updating from an initial to a final model. Because our target is to perform UQ for time-lapse imaging, we use a local wave-equation solver that allows us to solve the wave equation in a small subset of our entire model, thereby keeping computational costs low. Numerical examples demonstrate that incorporating the training step for NF provides a distribution of model perturbations, which is dependent on a designated prior, to quantify the uncertainty of FWI results.

We consider application of full-waveform inversion (FWI) to radio-frequency electromagnetic (EM) data. Radio-frequency imaging (RIM) is a cross-borehole technique to image EM subsurface properties from measurements of transmitted radio-frequency waves. It is used in coal seam imaging, ore exploration and various engineering and civil engineering applications. RIM operates at frequencies from 50 kHz to several tens of MHz. It differs from other geophysical EM methods, because the frequency band includes the transition between the wave propagation and diffusion regimes. RIM data are acquired in 2-D cross-hole sections in a reciprocal manner. Traditionally, radio-frequency data are inverted by straight-ray tomography because it is inexpensive and easy to implement. It is argued that due to attenuation, the sensitivity of the transmitted electric field is the strongest within the first Fresnel zone of the ray connecting the transmitter and receiver. While straight-ray tomography is a simple method to implement and fast, the nonlinearity in the relationship between model parameters and data is often strong enough to warrant nonlinear inversion techniques. FWI is an iterative high-resolution technique, in which the physical properties are updated to minimize the misfit between the measured and modelled wavefields. Full-waveform techniques have been used and extensively studied for the inversion of seismic data, and more recently, they have been applied to the inversion of ground penetrating radar data. Nonlinear inversion methods for RIM data are less advanced. Their use has been hindered by the high cost of full-wave modelling and the high conductivity contrasts of many RIM targets, and, to some extent, by the limitations of the measuring instruments. We present the first application of this methodology to perform simultaneous conductivity and permittivity inversion of RIM data. We implement the inversion in the frequency domain in two dimensions using Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) optimization. We analyse the sensitivity of the data to the model parameters and the parameter trade-off and validate the proposed methodology on a synthetic example with moderate conductivity variations and localized highly conductive targets. We then apply the FWI methodology to a field data set from Sudbury, Canada. For the field data set, we determine the most appropriate pre-processing steps that take into account specific peculiarities of RIM: the insufficient prior information about the subsurface and the limitations of the measuring equipment. We show that FWI is applicable under the conditions of RIM and is robust to imperfect prior knowledge: we obtain satisfactory model recoveries starting from homogeneous initial models in all of our examples. Just as other methods, FWI underestimates large conductivity contrasts due to the loss of sensitivity of the transmitted electric field to the conductivity variations as the conductivity increases above a certain level. The permittivity inside high conductors cannot be recovered, however, recovering permittivity variations in the resistive zones helps obtain better focused conductivity images with fewer artefacts. Overall, FWI produces cleaner, less noisy and higher resolution reconstructions than the methods currently used in practice.

We develop a comprehensive study involving three different types of machine learning (unsupervised, supervised, and semisupervised, which we emphasize) for bedrock-lithology classification using a publicly available data set from New South Wales, Australia. The goal of this work is to demonstrate (1) the value each different type of machine learning can provide and (2) which machine learning type(s) may be preferable under different circumstances. Training data are characteristically limited for geoscience problems, which makes supervised techniques susceptible to overfitting; we explore if semisupervised methods can perform better in these circumstances. Using the geophysical data and geologic map provided for the study area, we compare the performance of two supervised methods (the Light Gradient Boosting Machine and eXtreme Gradient Boosting) with one semisupervised algorithm (label propagation [LP]) in three scenarios with varied limited a priori lithologic constraints (i.e., the training data). Hyperparameter tuning is an essential component of supervised and semisupervised techniques, and the default procedure is to choose the hyperparameter combination with the largest mean cross-validation score. However, we use a new hyperparameter selection strategy that simultaneously uses the mean and standard deviation scores, and we test this new tactic for supervised and semisupervised methods. The results indicate (1) that the new hyperparameter selection technique can slightly improve the performance for supervised and semisupervised methods by 1%–2% compared with the standard selection approach and (2) that LP can outperform the two supervised methods by up to 10%, but it depends on how the training data are distributed. As for the unsupervised analysis, the clusters indicate heterogeneous regions that correlate well with the high-entropy areas in the supervised and semisupervised results. The clustering provides complementary results to the other two types of machine learning and is a source of supporting evidence for suggesting where more in-depth field mapping may be needed.

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.

Abstract When two waves interact within a rock sample, the interaction strength depends strongly on the sample’s microstructural properties, including the orientation of the sample layering. The study that established this dependence on layering speculated that the differences were caused by cracks aligned with the layers in the sample. To test this, we applied a uniaxial load to similar samples of Crab Orchard Sandstone and measured the nonlinear interaction as a function of the applied load and layer orientation. We show that the dependence of the nonlinear signal changes on applied load is exponential, with a characteristic load of 11.4–12.5 MPa that is independent of sample orientation and probe wavetype (P or S); this value agrees with results from the literature, but does not support the cracks hypothesis.

<p>Due to the nonlinearity of inversion as well as the noise in the data, seismic inversion results certainly have uncertainties. Whether quantifying these uncertainties is useful depends at least in part on the computational cost of computing them.  Bayesian techniques dominate uncertainty quantification for seismic inversion.  The goal of these methods is to estimate the probability distribution of the model parameters given the observed data. The Markov Chain Monte Carlo algorithm is widely employed for approximating the posterior distribution. However, generating the posterior samples by combining the prior and the likelihood is intractable for large problems and challenging for smaller problems. We apply a machine learning method called normalizing flows, which consists of a series of invertible and differentiable transformations, as an alternative to the sampling-based methods. In our work, the normalizing flows method is combined with full waveform inversion(FWI) using a numerically exact local solver to quantify the uncertainty of time-lapse changes. We integrate uncertainty quantification(UQ) and FWI by estimating UQ on the images generated by FWI making it computationally practical. In this way, a reasonable posterior probability distribution is directly predicted and produced by transforming from a normal distribution, measuring the amount and spread of variation in FWI images by sample mean and standard deviation. In our numerical results, the method for calculating the posterior distribution of the model is verified to be practical and advantageous in terms of effectiveness.</p>

Reverse time migration (RTM), as a state-of-the-art imaging technique, provides outstanding imaging capabilities due to its use of a full wave equation. Least-squares RTM (LSRTM) seeks the solution of a linearized wave equation via the minimization of a data misfit term; however, the quality of the results de-creases when the assumptions of the method are not satisfied. This occurs, for example, when we use an erroneous velocity model or inadequate physics for inverting the data. In such cases, appropriate regularization is required to mitigate these shortcomings and stabilize the LSRTM solution. However, even for structurally simple earth models, the reflectivity images are complicated and may not be explained properly by particular regularization methods such as the Tikhonov, total variation (TV), or sparse regularization. Reflectivity images can be thought of as the difference between two structurally simpler components: a piecewise-constant component (the true squared slowness) and a smooth component (the background model). We have used a combined Tikhonov-TV regularizer to regular-ize these components separately, leading to an effective regulari-zation for the reflectivity image. Because the background model is known in advance, this combined regularization reduces to a shifted TV regularization for which the associated optimization problem is solved efficiently using a new implementation of the Bregmanized operator splitting algorithm applied to the shifted TV method and the usual TV method. We determine the perfor-mance of our method with a set of numerical examples. The results confirm that our shifted regularization increases the robustness of LSRTM and allows us to estimate high-quality reflectivity images and properly update the background velocity model.

SUMMARY Dynamic nonlinear elasticity in rocks may play an important role in earth processes, such as earthquake nucleation. In order to understand how nonlinear elasticity occurs within the shallow crust, experiments are required that simulate the in situ conditions of intact crustal rocks. Additionally, exploring the behaviour of nonlinear elasticity in response to changes in external parameters (e.g. temperature and wave frequency) acts as a means to further illuminate the complex mechanisms which give rise to nonlinear elasticity in rocks. In this study, we perform dynamic acoustoelastic testing (DAET) experiments on an intact cataclasite from the damage zone of the Alpine Fault, New Zealand. By performing pump-probe DAET experiments inside a temperature-controlled chamber, we are able to investigate a rich variety of nonlinear behaviour as a function of temperature. We find that the magnitude of average softening, cubic nonlinearity, and hysteresis tend to increase as temperature increases from 20 to 110 °C. In contrast, quadratic nonlinearity decreases with increasing temperature. These observations support the hypothesis that at least two distinct mechanisms control nonlinear phenomena in rocks. Nonlinear parameters show little to no dependence on frequency over the 200–600 Hz pump range, although values of the nonlinear parameter α are found to be nearly two orders of magnitude smaller than those determined using ultrasonic perturbations. Additionally, an analysis using different time windows shows that the surface waves of the ultrasonic probe sense greater nonlinearity compared to the direct P- wave due to differences in the polarization and propagation paths. As well as providing further insight into the mechanisms responsible for nonlinear elasticity in rocks, our experiments show that nonlinear softening will increase as temperature increases in the damage zones of faults. This has potential implications for understanding earthquake nucleation processes.

Determining if uncertainty quantification is worth it or not is closely related to how that uncertainty is computed and the associated computational cost. For seismic imaging, it is typically done using Markov chain Monte Carlo algorithms (McMC). Solving an inverse problem using McMC means exploring and characterizing the ensemble of all plausible models through more or less point-wise random walk in the data misfit landscape. This is typically done using Bayes’ theorem via the computation of a posterior probability density function. Even though this can sound naively simple, it can come with a significant computational burden given the dimension of the problem to be solved and the expense of the forward solver. This is because as the number of dimensions grow, there are exponentially more possible guesses the algorithm can make, while only a few of these models will be accepted as plausible. More advanced uncertainty quantification methods such as Hamiltonian Monte Carlo (HMC) could be beneficial because they can handle higher dimensions because efficient sampling of the model space through pseudo-mechanical trajectories in the data misfit landscape is expected. In order for an HMC algorithm to efficiently sample the model space of interest and provide meaningful uncertainty estimates, three hyper-parameters need to be tuned for trajectory design: the Leapfrog steps L, the Leapfrog stepsize ε, and the Mass Matrix M. There has been already work showing how one can choose L and ε; however designing the appropriate M is far more challenging. We consider a time-lapse seismic scenario and use a local acoustic solver for fast forward solutions. We then use Singular value decomposition, in the vicinity of the true model, to transform our time-lapse optimal model to a system of normal coordinates and use only a few of the eigenvalues and eigenvectors of the Hessian as oscillators. By doing so, we can efficiently understand the impact of the initial conditions and the choice of M and gain insight on how to design M in the standard system. This gives us an intuitive way to understand the mass matrix, allowing us to determine whether gains from the HMC algorithm are worth the cost of determining the parameters.

Summary We present a new target-enclosing full waveform inversion (FWI) method based on a interferometric objective function, where the inversion is driven by the misfit between the wavefields reconstructed in a subdomain of interest via the convolution and correlation representation formulas. The method is fully local in the sense that it does not depend on the reconstruction of the physical properties outside the local domain, and only requires a kinematic velocity model estimate for redatuming. We compare the proposed method to another fully local FWI method based on the convolution representation formula and to full-model surface-data FWI. We demonstrate the potential of the proposed interferometric full waveform inversion method to achieve higher resolution images at a comparable or lower cost than the other methods.

Pore geometry is an important parameter in reservoir characterization that affects the permeability of reservoirs and can also be a controlling factor on the impact of pressure and saturation on reservoirs elastic properties. We have used selective laser sintering 3D printing technology to build physical models to experimentally investigate the impacts of pore aspect ratio on compressional (P-) and shear wave (S-wave) velocities and amplitude variation with offset (AVO). We printed six models to study the effects of the pore aspect ratio of prolate and oblate pore structures on elastic properties and AVO signatures. We found that the P-wave velocity is reduced by decreasing the pore aspect ratio (flatter pore structure), whereas the S-wave velocity is less sensitive to the pore aspect ratio. This effect is reduced when the samples are water saturated. We developed new experimental and processing techniques to extract realistic AVO signatures from our experimental data and found that the pore aspect ratio has similar effects on AVO as fluid compressibility. This indicated that not considering the pore aspect ratio in AVO analysis can lead to misleading interpretations. We also found that these effects are reduced in water-saturated samples.

A comparison of three different isotropic non-linear elastic models uncovers subtle but important differences in the acoustoelastic responses of a material slab that is subjected to dynamic deformations during a pump-probe experiment. The probe wave deformations are small and are superimposed on larger underlying deformations using three different models: Landau–Lifshitz (using its fourth-order extension), compressible neo-Hookean model (properly accounting for volumetric deformations), and an alternative neo-Hookean model (fully decoupled energies due to distortional isochoric and volumetric deformations). The analyses yield elasticity tensors and respective expressions for the propagation speeds of P-wave and S-wave probes for each model. Despite having many similarities, the different models give different predictions of which probe wave types will have speeds that are perturbed by different pump wave types. The analyses also show a conceptual inconsistency in the Landau–Lifshitz model, that a simple shear deformation induces a stress and a shear wave probe speed that depend on the second-order elastic constant λ , which controls resistance to volumetric changes and thus should not be present in the expressions for shear stress and shear wave probe speeds. Thus, even though the Landau–Lifshitz model is widely used, it may not always be the best option to model experimental data.

Summary Radio-frequency imaging (RIM) is a cross-borehole technique to image electromagnetic subsurface properties from measurements of radio-frequency waves. RIM operates at mid-range frequencies and has most applications in mining. Traditionally, RIM problem has been solved by straight ray tomography. Recently, an inverse scattering method has been proposed demonstrating the potential for higher-resolution images by incorporating more exact physics into the inversion process. We present an application of full waveform inversion (FWI) to conductivity imaging with RIM data. FWI is a high resolution technique, in which the physical property is updated iteratively to minimize the misfit between the measured and modelled wavefields. The full waveform modelling with Maxwell’s equations is efficiently implemented in the frequency domain. The model update is calculated by the L-BFGS method, where the gradient is evaluated by the adjoint state technique. We show that the resolution of a half-width of the first Fresnel zone is achievable to correctly recover the shape and location of conductive targets. Large conductivity contrasts are underestimated due to attenuation of the wavefields in highly conductive zones. The method can be extended to include electric permittivity inversion.

Least-squares reverse time migration (LSRTM) is a leading tool in imaging complex geological structures. This is because it includes the full-wave equation via reverse time migration (RTM), and solves the linearized wave equation by the least-squares (LS) minimization. To stabilize the LSRTM solution and mitigate the shortcomings that arise from linearization, simplification, approximation in simulation, and imperfection of data acquisition an appropriate regularization must be used. However, the conventional regularization methods such as Tikhonov, total variation (TV), or sparse regularization perform suboptimally, even in recovering structurally simple earth models. Here we introduce a new extension of Tikhonov-TV compound regularization, called shifted TV, to regularize the unknown reflectivity image. This new regularization implemented in the frame of a nonlinear migration allows for high-resolution and stable imaging in the presence of rough migration velocity models. We demonstrate the performance of our method using a short offset academic-style streamer data generated in the crustal-scale GO 3D OBS synthetic subduction zone model. The results confirm that the shifted regularization increases the robustness of LSRTM with an imperfect background velocity model and allows us to estimate high-quality reflectivity images of complex geological setting despite the limited streamer length. Note: This paper was accepted into the Technical Program but was not presented at IMAGE 2022 in Houston, Texas.

Summary We present a new target-enclosing full-waveform inversion method based on a new objective function. This objective function is based on the mismatch between wavefields reconstructed in the target domain with the convolution and correlation representation formulas, using data from the boundary of the target domain. The proposed method requires only kinematic knowledge of the subsurface model, particularly the overburden for redatuming, and no reconstruction of the model outside of the target area. In this sense it is truly local. We show a synthetic example with exactly redatumed wavefields and compare it with conventional FWI of surface data. The potential of the proposed method is demonstrated both in terms of the quality of target recovery and reduction in computational cost.

The Jizhong depression contains several geothermal reservoirs that are characterized by localized low-velocity anomalies. In this article, full-waveform inversion (FWI) is used to characterize these anomalies and determine their extent. This is a challenging problem because the reservoirs are quite small and the available data have usable frequencies only down to 5 Hz. An accurate-enough starting model is carefully built by using an iterative travel time tomography method combined with a cycle-skipping assessment method to begin the inversion at 5 Hz. A multiscale Laplace–Fourier-domain FWI with a layer-stripping approach is implemented on the starting model by gradually increasing the maximum offset. The result of overlapping the recovered velocity model on the migrated seismic profile shows a good correlation between the two results. The recovered model is assessed by ray tracing, synthetic seismogram modeling, checkerboard testing and comparisons with nearby borehole data. These tests indicate that low-velocity anomalies down to a size of 0.3 km × 0.3 km at a maximum depth of 2 km can be recovered. Combined with the well log data, the resulting velocity model allows us to delineate two potential geothermal resources, one of which was previously unknown.

Full waveform inversion (FWI) is beginning to be used to characterize weak seismic events at different scales, an example of which is microseismic event (MSE) characterization. However, FWI with unknown sources is a severely underdetermined optimization problem, and hence requires strong prior information about the sources and/or the velocity model. The frequency-domain wavefield reconstruction inversion method (WRI) has shown promising results to mitigate the nonlinearity of the FWI objective function that is generated by cycle-skipping. WRI relies on the reconstruction of data-assimilated wavefields, which approach the true wavefields near the receivers, a helpful feature when the source is added as an additional optimization variable. We present an adaptation of a recently proposed version of WRI based on the alternating direction method of multipliers (ADMM) that first finds the location of the MSEs and then reconstructs the wavefields and the source signatures jointly. Finally, the subsurface model is updated to focus the MSEs at their true location. The method does not require prior knowledge of the number of MSEs. The inversion is stabilized by sparsifying regularizations separately tailored to the source location and velocity model subproblems. The method is tested on the Marmousi model using one MSE and two clusters of MSEs with two different initial velocity models, an accurate one and a rough one, as well as with added noise. In all cases, the method accurately locates the MSEs and recovers their source signatures.

Full-waveform inversion (FWI) is a high-resolution and computationally intensive imaging technique to reconstruct unknown parameters in the computational model in which the waves propagate; however, an accurate model of only part of this medium is required for some applications. To decrease the computational burden of such problems, target-oriented FWI was proposed where the redatumed data on the part of the medium or localized solvers for the wave equation are used. On the other hand, the classical formulation of FWI suffers from non-linearity and ill-posedness, which makes FWI sensitive to the initial model, the low-frequency content of the data, and limited illumination. In this study, we propose a localized version of the alternating direction method of multipliers (ADMM)-based FWI method, which was proposed to solve these problems in classical FWI. In our localized FWI or LWI, the medium is decomposed into a few subdomains, where some of them are updated, and the others are kept fixed based on an adaptation of multi-block ADMM, which is a powerful algorithm for solving inverse problems with decomposition and block separability. Numerical tests on the Marmousi model for a time-lapse application confirm the computational efficiency and robustness against background velocity model errors.

We present a novel methodology for exploring 4D seismic data in the context of monitoring subsurface resources. Data‐space exploration is a key activity in scientific research, but it has long been overlooked in favor of model‐space investigations. Our methodology performs a data‐space exploration that aims to define structures in the covariance matrix of the observational errors. It is based on Bayesian inferences, where the posterior probability distribution is reconstructed through trans‐dimensional (trans‐D) Markov chain Monte Carlo sampling. The trans‐D approach applied to data‐structures (termed ”partitions”) of the covariance matrix allows the number of partitions to freely vary in a fixed range during the McMC sampling. Due to the trans‐D approach, our methodology retrieves data‐structures that are fully data‐driven and not imposed by the user. We applied our methodology to 4D seismic data, generally used to extract information about the variations in the subsurface. In our study, we make use of real data that we collected in the laboratory, which allows us to simulate different acquisition geometries and different reservoir conditions. Our approach is able to define and discriminate different sources of noise in 4D seismic data, enabling a data‐driven evaluation of the quality (so‐called “repeatability”) of the 4D seismic survey. We find that: (a) trans‐D sampling can be effective in defining data‐driven data‐space structures; (b) our methodology can be used to discriminate between different families of data‐structures created from different noise sources. Coupling our methodology to standard model‐space investigations, we can validate physical hypothesis on the monitored geo‐resources.