Monitoring carbon dioxide (CO2) injected and stored in subsurface reservoirs is critical for avoiding failure scenarios and enables real-time optimization of CO2injection rates. Bayesian sequential data assimilation (DA) is a statistical method for combining information over time from multiple sources to estimate a hidden state, such as the spread of the subsurface CO2plume. Existing literature in the seismic-CO2monitoring domain uses small physical systems with unscalable DA algorithms, ignores the CO2flow dynamics, or simulates seismic data without the wave equation or with unrealistic survey designs. We improve upon existing DA literature in this domain by applying the scalable ensemble Kalman filter (EnKF) DA algorithm to a high-dimensional CO2reservoir using two-phase flow dynamics and time-lapse full waveform seismic data with a realistic surface-seismic survey design. We show this DA method is more accurate compared to using either the seismic data or the fluid physics alone. Furthermore, we show the stability of this method by testing a range of values for the EnKF hyperparameters and give guidance on their selection for seismic CO2reservoir monitoring.

Seismic imaging in challenging sedimentary basins and reservoirs requires acquiring, processing, and imaging very large volumes of data (tens of terabytes). To reduce the cost of acquisition and the time from acquiring the data to producing a subsurface image, novel acquisition systems based on compressive sensing (CS), low-rank matrix recovery (LRMR), and randomized sampling have been developed and implemented. These approaches allow practitioners to achieve dense wavefield reconstruction from a substantially reduced number of field samples. However, designing acquisition surveys suited for this new sampling paradigm remains a critical and challenging role in oil, gas, and geothermal exploration. Typical random designs studied in the LRMR and CS literature are difficult to achieve by standard industry hardware. For practical purposes, a compromise between stochastic and realizable samples is needed. In this article, we propose a deterministic and computationally cheap tool to alleviate randomized acquisition design, prior to survey deployment and large-scale optimization. We consider universal and deterministic matrix completion results in the context of seismology, where a bipartite graph representation of the source–receiver layout allows for the respective spectral gap (SG) to act as a quality metric for wavefield reconstruction. We provide realistic scenarios to demonstrate the utility of the SG as a flexible tool that can be incorporated into existing survey design workflows for successful seismic data acquisition via low-rank and sparse signal recovery.

Compressive sensing has driven a resurgence of sparse recovery algorithms with $\ell_1$-norm minimization. While these minimizations are relatively well understood for small underdetermined, possibly inconsistent systems, their behavior for large overdetermined and inconsistent systems has received much less attention. Specifically, we focus on large systems where computational restrictions call for algorithms that use randomized subsets of rows that are touched a limited number of times. In that regime, $\ell_1$-norm minimization algorithms exhibit unwanted fluctuations near the desired solution, and the linear Bregman iterations are no exception. These fluctuations result in uncertainty about the recovery results, forcing increased effort such as longer run times or other additional search efforts. We explain this observed lack of performance in terms of chatter, a well-known phenomenon observed in nonsmooth dynamical systems, where intermediate solutions wander between different states, stifling convergence. By identifying chatter as the culprit, we then reduce it by modifying the Bregman iterations with adaptive elementwise step lengths combined with potential support detection via threshold crossing. We demonstrate the performance of our algorithm on carefully selected stylized examples that mimic real large scale problems and on a realistic seismic imaging problem involving millions of unknowns and matrix-free matrix-vector products that involve expensive wave-equation solves.

Stencil computations are a key part of many high-performance computing applications, such as image processing, convolutional neural networks, and finite-difference solvers for partial differential equations. Devito is a framework capable of generating highly optimized code given symbolic equations expressed in Python , specialized in, but not limited to, affine (stencil) codes. The lowering process—from mathematical equations down to C++ code—is performed by the Devito compiler through a series of intermediate representations. Several performance optimizations are introduced, including advanced common sub-expressions elimination, tiling, and parallelization. Some of these are obtained through well-established stencil optimizers, integrated in the backend of the Devito compiler. The architecture of the Devito compiler, as well as the performance optimizations that are applied when generating code, are presented. The effectiveness of such performance optimizations is demonstrated using operators drawn from seismic imaging applications.

Adapting the cloud for high-performance computing (HPC) is a challenging task, as software for HPC applications hinges on fast network connections and is sensitive to hardware failures. Using cloud infrastructure to recreate conventional HPC clusters is therefore in many cases an infeasible solution for migrating HPC applications to the cloud. As an alternative to the generic lift and shift approach, we consider the specific application of seismic imaging and demonstrate a serverless and event-driven approach for running large-scale instances of this problem in the cloud. Instead of permanently running compute instances, our workflow is based on a serverless architecture with high throughput batch computing and event-driven computations, in which computational resources are only running as long as they are utilized. We demonstrate that this approach is very flexible and allows for resilient and nested levels of parallelization, including domain decomposition for solving the underlying partial differential equations. While the event-driven approach introduces some overhead as computational resources are repeatedly restarted, it inherently provides resilience to instance shut-downs and allows a significant reduction of cost by avoiding idle instances, thus making the cloud a viable alternative to on-premise clusters for large-scale seismic imaging.

Large-scale parameter estimation problems are among some of the most computationally demanding problems in numerical analysis. An academic researcher’s domain-specific knowledge often precludes that of software design, which results in inversion frameworks that are technically correct but not scalable to realistically sized problems. On the other hand, the computational demands for realistic problems result in industrial codebases that are geared solely for high performance, rather than comprehensibility or flexibility. We propose a new software design for inverse problems constrained by partial differential equations that bridges the gap between these two seemingly disparate worlds. A hierarchical and modular design reduces the cognitive burden on the user while exploiting high-performance primitives at the lower levels. Our code has the added benefit of actually reflecting the underlying mathematics of the problem, which lowers the cognitive load on the user using it and reduces the initial startup period before a researcher can be fully productive. We also introduce a new preconditioner for the 3D Helmholtz equation that is suitable for fault-tolerant distributed systems. Numerical experiments on a variety of 2D and 3D test problems demonstrate the effectiveness of this approach on scaling algorithms from small- to large-scale problems with minimal code changes.

We propose an extended full-waveform inversion formulation that includes general convex constraints on the model. Though the full problem is highly nonconvex, the overarching optimization scheme arrives at geologically plausible results by solving a sequence of relaxed and warm-started constrained convex subproblems. The combination of box, total variation, and successively relaxed asymmetric total variation constraints allows us to steer free from parasitic local minima while keeping the estimated physical parameters laterally continuous and in a physically realistic range. For accurate starting models, numerical experiments carried out on the challenging 2004 BP velocity benchmark demonstrate that bound and total variation constraints improve the inversion result significantly by removing inversion artifacts, related to source encoding, and by clearly improved delineation of top, bottom, and flanks of a high-velocity high-contrast salt inclusion. The experiments also show that for poor starting models these two constraints by themselves are insufficient to detect the bottom of high-velocity inclusions such as salt. Inclusion of the one-sided asymmetric total variation constraint overcomes this issue by discouraging velocity lows to buildup during the early stages of the inversion. To the best of the authors' knowledge the presented algorithm is the first to successfully remove the imprint of local minima caused by poor starting models and band-width limited finite aperture data.

Seismic imaging can be formulated as a linear inverse problem where a medium perturbation is obtained via minimization of a least-squares misfit functional. The demand for higher resolution images in more geophysically complex areas drives the need to develop techniques that handle problems of tremendous size with limited computational resources. While seismic imaging is amenable to dimensionality reduction techniques that collapse the data volume into a smaller set of “super-shots”, these techniques break down for complex acquisition geometries such as marine acquisition, where sources and receivers move during acquisition. To meet these challenges, we propose a novel method that combines sparsity-promoting (SP) solvers with random sub-set selection of sequential shots, yielding a SP algorithm that only ever sees a small portion of the full data, enabling its application to very large-scale problems. Application of this technique yields excellent results for a complicated synthetic, which underscores the robustness of sparsity promotion and its suitability for seismic imaging.

Sparco is a framework for testing and benchmarking algorithms for sparse reconstruction. It includes a large collection of sparse reconstruction problems drawn from the imaging, compressed sensing, and geophysics literature. Sparco is also a framework for implementing new test problems and can be used as a tool for reproducible research. Sparco is implemented entirely in Matlab, and is released as open-source software under the GNU Public License.