---
title: "De-risking geological carbon storage from high resolution time-lapse seismic to explainable leakage detection"
author: |
    Ziyi Yin, Huseyin Tuna Erdinc, Abhinav Prakash Gahlot, Mathias Louboutin, Felix J. Herrmann \
    Georgia Institute of Technology \
bibliography:
    - paper.bib
---

## Abstract:

Geological carbon storage represents one of the few truly scalable technologies capable of reducing the CO~2~ concentration in the atmosphere. While this technology has the potential to scale, its success hinges on our ability to mitigate its risks. An important aspect of risk mitigation concerns assurances that the injected CO~2~ remains within the storage complex. Amongst the different monitoring modalities, seismic imaging stands out with its ability to attain high resolution and high fidelity images. However, these superior features come, unfortunately, at prohibitive costs and time-intensive efforts potentially rendering extensive seismic monitoring undesirable. To overcome this shortcoming, we present a methodology where time-lapse images are created by inverting non-replicated time-lapse monitoring data jointly. By no longer insisting on replication of the surveys to obtain high fidelity time-lapse images and differences, extreme costs and time-consuming labor are averted. To demonstrate our approach, hundreds of noisy time-lapse seismic datasets are simulated that contain imprints of regular CO~2~ plumes and irregular plumes that leak. These time-lapse datasets are subsequently inverted to produce time-lapse difference images used to train a deep neural classifier. The testing results show that the classifier is capable of detecting CO~2~ leakage automatically on unseen data and with a reasonable accuracy.

## Introduction

For various reasons, seismic monitoring of geological carbon storage (GCS) comes with its own set of unique challenges. Amongst these challenges, the need for low-cost highly repeatable, high resolution, and high fidelity images ranks chiefly. While densely sampled and replicated time-lapse surveys---which rely on permanent reservoir monitoring systems---may be able to provide images conducive to interpretation and reservoir management, these approaches are often too costly and require too much handholding to be of practical use for GCS.

To overcome these challenges, we replace the current paradigm of costly replicated acquisition, cumbersome time-lapse processing, and interpretation, by a joint inversion framework mapping time-lapse data to high fidelity and high resolution images from sparse non-replicated time-lapse surveys. We demonstrate that we arrive at an imaging framework that is suitable for automatic detection of pressure-induced CO~2~ leakage. Rather than relying on meticulous 4D workflows where baseline and monitoring surveys are processed separately to yield accurate and artifact-free time-lapse differences, our approach exposes information that is shared amongst the different vintages by formulating the imaging problem in terms of an unknown fictitious common component, and innovations of the baseline and monitor survey(s) with respect to this common component. Because the common component is informed by all time-lapse surveys, its image quality improves when the surveys bring complementary information, which is the case when the surveys are not replicated. In turn, the enhanced common component results in improved images for the baseline, monitor, and their time-lapse difference. Joint inversion also leads to robustness with respect to noise, calibration errors, and time-lapse changes in the background velocity model.

To showcase the achievable imaging gains and how these can be used in a GCS setting where CO~2~ leakage is of major consideration, we create hundreds of time-lapse imaging experiments involving CO~2~ plumes whose behavior is determined by the two-phase flow equations. To mimic irregular flow due to pressure-induced opening of fractures, we increase the permeability in the seal at random locations and pressure thresholds. The resulting flow simulations are used to generate time-lapse datasets that serve as input to our joint imaging scheme. The produced time-lapse difference images are subsequently used to train and test a neural network that as an explainable classifier determines whether the CO~2~ plume behaves regularly or shows signs of leakage.

Our contributions are organized as follows. First, we discuss the time-lapse seismic imaging problem and its practical difficulties. Next, we introduce the joint recovery model that takes explicit advantage of information shared by multiple surveys. By means of a carefully designed synthetic case study involving saline aquifers made of Blunt sandstone in the Southern North Sea, we demonstrate the uplift of the joint recovery model and how its images can be used to train a deep neural network classifier to detect erroneous growth of the CO~2~ plume automatically. Aside from determining whether the CO~2~ plume behaves regularly or not, our network also provides class activation mappings that visualize areas in the image on which the network is basing its classification. 

## Seismic monitoring with time-lapse imaging

To keep track of CO~2~ plume development during geological carbon storage (GCS) projects, multiple time-lapse surveys are collected. Baseline surveys are acquired before the supercritical CO~2~ is injected into the reservoir. These baseline surveys, denoted by the index ``j=1``, are followed by one or more monitor surveys, collected at later times and indexed by ``j=2,\cdots,n_v`` with ``n_v`` the total number of surveys.

Seismic monitoring of GCS brings its own unique set of challenges that stem from the fact that its main concern is (early) detection of possible leakage of CO~2~ from the storage complex. To be successful with this task, monitoring GCS calls for a time-lapse imaging modality that is capable of

- detecting weak time-lapse signals associated with small rock-physics changes induced by CO~2~ leakage
- attaining high lateral resolution from active-source surface seismic data to detect vertically moving leakage
- handling an increasing number of seismic surveys collected over long periods of time (``\sim50-100`` years)
- reducing costs drastically by no longer insisting on replication of time-lapse surveys to attain high degrees of repeatability 
- lowering the cumulative environmental imprint of active source acquisition

### Monitoring with the joint recovery model

To meet these challenges, we choose a linear imaging framework where observed linearized data for each vintage is related to perturbations in the acoustic impedance via

```math #eq-lin-model
\mathbf{b}_j=\mathbf{A}_j\mathbf{x}_j+\mathbf{e}_j\quad \text{for}\quad j=1,2,\cdots,n_v.
```

In this expression, the matrix ``\mathbf{A}_j`` stands for the linearized Born scattering operator for the ``j\,\mathrm{th}`` vintage. Observed linearized data, collected for all shots in the vector ``\mathbf{b}_j``, is generated by applying the ``\mathbf{A}_j``'s to the (unknown) impedance perturbations denoted by ``\mathbf{x}_j`` for ``j=1,2,\cdots, n_v``. The task of time-lapse imaging is to create high resolution, high fidelity and true amplitude estimates for the time-lapse images, ``\left\{\widehat{\mathbf{x}}_j\right\}_{j=1}^{n_v}``, from non-replicated sparsely sampled noisy time-lapse data.

We argue that our choice for linearized imaging is justified for four reasons. First, CO~2~-injection sites undergo extensive baseline studies, which means that accurate information on the background velocity model is generally available. Second, changes in the acoustic parameters induced by CO~2~ injection are typically small, so it suffices to work with one and the same background model for the baseline and monitor surveys. Third, when the background model is sufficiently close to the true model, linearized inversion, which corresponds to a single Gauss-Newton iteration of full-waveform inversion, converges quadratically. Fourth, because the forward model is linear, it is conducive to the use of the joint recovery model where inversions are carried out with respect to the common component, which is shared between all vintages, and innovations with respect to the common component. Because the common component represents an average, we expect this joint imaging method to be relatively robust with respect to kinematic changes induced by time-lapse effects or by lack of calibration of the acquisition [@oghenekohwo2017hrt]\.

By parameterizing time-lapse images, ``\left\{\mathbf{x}_j\right\}_{j=1}^{n_v}``, in terms of the common component, ``\mathbf{z}_0``, and innovations with respect to the common component, ``\left\{\mathbf{z}_j\right\}_{j=1}^{n_v}``, we arrive at the joint recovery model where representations for the images are given by

```math #eq-components
\mathbf{x}_j = \frac{1}{\gamma}\mathbf{z}_0 + \mathbf{z}_j \quad \text{for} \quad j=1,2,\cdots,n_v.
```

Here, the parameter, ``\gamma``, controls the balance between the common component, ``\mathbf{z}_0``, and innovation components, ``\left\{\mathbf{z}_j\right\}_{j=1}^{n_v}`` [@li2015weighted]. Compared to traditional time-lapse approaches, where data are imaged separately or where time-lapse surveys are subtracted, inversions for time-lapse images based on the above parameterization are carried out jointly and involve inverting the following matrix:

```math #eq-jrm
\mathbf{A} = \begin{bmatrix} 
\frac{1}{\gamma} \mathbf{A}_1 & \mathbf{A}_1 &  & \\ 
\vdots & & \ddots &   \\
\frac{1}{\gamma} \mathbf{A}_{n_v} & &  & \mathbf{A}_{n_v}
\end{bmatrix}.
```
While traditional time-lapse imaging approaches strive towards maximal replication between the surveys to suppress acquisition related artifacts, imaging with the joint recovery model---which entails inverting the underdetermined system in Equation #eq-jrm using structure-promotion techniques (e.g. via ``\ell_1``-norm minimization)---improves the image quality of the vintages themselves in situations where the surveys are not replicated. This occurs in cases where ``\mathbf{A}_i\neq\mathbf{A}_j`` for ``\forall i\neq j``, or in situations where there is significant noise. This remarkable result was shown to hold for sparsity-promoting denoising of time-lapse field data [@wei2018improve;@tian2018joint], for various wavefield reconstructions of randomized simultaneous-source dynamic (towed-array) and static (OBC/OBN) marine acquisitions [@oghenekohwo2017hrt;@kotsi2020time;@zhou2021non], and for wave-based inversion, including least-squares reverse-time migration and full-waveform inversion [@oghenekohwo2017THetl;@oghenekohwo2017EAGEitl]\. The observed quality gains in these applications can be explained by improvements in the common component resulting from complementary information residing in non-replicated time-lapse surveys. This enhanced recovery of the common component in turn improves the recovery of the innovations and therefore the vintages themselves. The time-lapse differences themselves also improve, or at the very least, remain relatively unaffected when the surveys are not replicated. Relaxing replication of surveys obviously leads to reduction in cost and environmental impact. Below, we show how GCS monitoring also benefits from this approach.

### Monitoring with curvelet-domain structure promotion

To obtain high resolution and high fidelity time-lapse images, we invert the system in Equation #eq-jrm [@yang2020tdsp;@witte2018cls;@yin2021SEGcts] with 

```math #eq-elastic
\underset{\mathbf{z}}{\operatorname{minimize}} \quad \lambda \|\mathbf{C}\mathbf{z}\|_1+\frac{1}{2}\|\mathbf{C}\mathbf{z}\|_2^2 \\
\text{subject to}\quad \|\mathbf{b}- \mathbf{A}\mathbf{z}\|_2^2 \leq \sigma,
```

where ``\mathbf{C}`` is the forward curvelet transform, ``\lambda`` the threshold parameter, and ``\sigma`` the magnitude of the noise. At iteration ``k`` and for ``\sigma=0``, solving Equation #eq-elastic corresponds to computing the following iterations:

```math #LBk
\begin{array}{lcl} 
  \mathbf{u}_{k+1} & = & \mathbf{u}_k-t_k \mathbf{A}_k^\top(\mathbf{A}_k\mathbf{z}_{k}-\mathbf{b}_k)\\
 \mathbf{z}_{k+1} & = & \mathbf{C}^\top S_{\lambda}(\mathbf{C}\mathbf{u}_{k+1}),
\end{array}
```

where ``\mathbf{A}_k``, with a slight abuse of notation, represents the matrix in Equation #eq-jrm for a subset of shots randomly selected from sources in each survey. The vector ``\mathbf{b}_k`` contains the extracted shot records from ``\mathbf{b}`` and the symbol ``^\top`` refers to the adjoint. Sparsity is promoted via curvelet-domain soft thresholding, ``S_{\lambda}(\cdot)=\max(|\cdot|-\lambda,0)\operatorname{sign}(\cdot)``, where, ``\lambda``, is the threshold. The vectors ``\mathbf{u}_k`` and ``\mathbf{z}_k`` contain the baseline and innovation components. 

<!-- When different vintages are imaged independently, time-lapse monitoring involves minimization of

```math #objtlimaging
\underset{\mathbf{x}_j}{\operatorname{minimize}} \quad \|\mathbf{b}_j-\mathbf{A}_j\mathbf{x}_j\|_2^2,\quad \text{for}\quad j=1,2,\cdots,n_v.
```

In this formulation, the ``\{\mathbf{b}_j\}_{j=1}^{n_v}`` and ``\{\mathbf{x}_j\}_{j=1}^{n_v}`` represent the observed time-lapse data and the unknown time-lapse images for the different vintages. The matrices ``\{\mathbf{A}_j\}_{j=1}^{n_v}`` represent the discretized linearized Born modeling operators, which relate the perturbations in the acoustic impedance to the linearized data for each vintage. Minimizing Equation #objtlimaging corresponds to carrying out independent least-squares migrations [@herrmann11GPelsqIm] for each vintage separately. When the same background model is used to parameterize the Born modeling operators for all vintages, the CO~2~-induced acoustic impedance change can be computed by subtraction of the recovered seismic images between different points in time. Since the difference of the seismic images due to CO~2~ plume is subtle, conventional time-lapse seismic methods usually call for dense replicated acquisitions between the baseline and monitor surveys [@zhang2013double;@yang2015double;@zhou2021central], which can be expensive and difficult to achieve in practice in the field.

In practice, the quality of the image from this subtraction is known to suffer from differences in seismic acquisition (source/receiver positions, source types), presence of noise in the data, and calibration errors [@oghenekohwo2017hrt;@kotsi2020time;@zhou2021non]. These complications might introduce artifacts into the recovered images, which can be wrongly attributed to time-lapse differences. It is also shown that these artifacts can be spurious given only a minor calibration error on the acquisition geometry of dense replicated surveys [@oghenekohwo2017hrt]. Instead of insisting on replicate the surveys, we take advantage of the joint recovery model (JRM) to focus on the commonalities between the different images and exploit the shared information between all surveys.

## Joint recovery model

Based on the fact that the Earth subsurface undergoes subtle change during the time range of the CCS project, the images for the different surveys are mostly the same, except for minor changes at the CO~2~ plume location. Therefore, we assume that there is a common component ``\mathbf{z}_0`` shared by all the images ``\{\mathbf{x}_j\}_{j=1}^{n_v}``, and there are innovation components with respect to the common component ``\{\mathbf{z}_j\}_{j=1}^{n_v}``, which are unique for the different images. The different components and the images are related as following:

```math #components
\mathbf{x}_j = \frac{1}{\gamma}\mathbf{z}_0 + \mathbf{z}_j \quad \text{for} \quad j=1,2,\cdots,n_v,
```

where the parameter ``\gamma`` controls the balance between common component and innovation components [@li2015weighted]. With this formulation, we make use of the joint recovery model (JRM) to solve for the common and innovation components instead of images themselves [@oghenekohwo2016GEOPctl;@wason2016GEOPctl]. This leads to inverting the following matrix:

```math #jrm
\widetilde{\mathbf{A}} = \begin{bmatrix} 
\frac{1}{\gamma} \mathbf{A}_1 & \mathbf{A}_1 &  & \\ 
\vdots & & \ddots &   \\
\frac{1}{\gamma} \mathbf{A}_{n_v} & &  & \mathbf{A}_{n_v}
\end{bmatrix}.
```

This matrix relates the common and innovation components collected in vector ``\mathbf{z}=\left[\mathbf{z}_0^{\top},\mathbf{z}_1^{\top},\cdots,\mathbf{z}_{n_v}^{\top}\right]^{\top}
``, to the linearized seismic dataset collected in vector ``\mathbf{b}=\left[\mathbf{b}_1^{\top},\cdots,\mathbf{b}_{n_v}^{\top}\right]^{\top}``. In JRM, the common component ``\mathbf{z}_0`` is related to all surveys through the first column of ``\widetilde{\mathbf{A}}``. Instead of replicating surveys to force the recoveries for each vintage to be more repeatable, JRM aims to find the shared information from different surveys. Therefore, when surveys are not replicated, the common component absorbs the complementary information from all surveys and can be built robustly with respect to noise and calibration errors [@oghenekohwo2017highly]. With a well recovered common component, the innovation components and the time-lapse differences are also expected to be better recovered. This has been confirmed by several synthetic experiments on wavefield recovery [@oghenekohwo2016GEOPctl;@wason2016GEOPctl], field data denoising [@wei2018improve;@tian2018joint], imaging and inversion [@oghenekohwo2017THetl], and recently on geological carbon storage monitoring [@yin2021SEGcts].

Apart from robustly monitoring the CO~2~ dynamics by time-lapse imaging, we further provide an explainable CO~2~ leakage detection network to identify anomalous behavior of the CO~2~ plume in the time-lapse images. -->

## Numerical case study: Blunt sandstone in the Southern North Sea 

Before discussing the impact of high resolution and high fidelity time-lapse imaging with the joint recovery model on the down-stream task of automatic leakage detection with a neural network classifier, we first detail the setup of our numerical experiments using techniques from simulation-based acquisition design as described by @yin2021SEGcts\. In order to generate realistic time-lapse data and training sets for the automatic leakage classifier, we follow the workflow summarized in Figure #fig:workflow\. In this approach, use is made of proxy models for seismic properties derived from real 3D imaged seismic and well data [@jones2012building]. With rock physics, these seismic models are converted to fluid-flow models that serve as input to two-phase flow simulations. The resulting datasets, which include pressure-induced leakage, will be used to create time-lapse data used to train our classifier. For more detail, refer to the caption of Figure #fig:workflow\.

<!-- We present a case study to validate the performance of the joint recovery model and demonstrate the efficacy of the deep neural network classifier for CO~2~ leakage detection. In this study, we conduct numerical experiments on the Compass model, a near-realistic model that is derived from the pre-stack seismic dataset and well logs of the North Sea area [@jones2012building]. This model is representative for the geology of saline aquifers in the Earth subsurface of the North Sea area, which is considered as a potential site for geological carbon storage. We first discuss the simulation-based monitoring design framework below (shown in Figure #fig:perm): -->

### Figure: {#fig:workflow}
![](figs/classification_workflow.png){width=98%}

:Simulation-based monitoring design framework. Starting with a proxy model for the wavespeed and density (a), the workflow proceeds by converting these seismic properties into permeability and porosity (b). These fluid flow properties are used to simulated CO~2~ plumes that behave regularly or exhibit leakage outside the storage complex (c). Induced changes by the CO~2~ plume for the wavespeed and density are depicted in (d) and serve as input to simulations of time-lapse seismic data (SNR ``8.0\,\mathrm{dB}``) and shot-domain time-lapse differences (SNR ``-31.4\,\mathrm{dB}``). Imaging results for regular and irregular plume developments are plotted in (f) and serve as input to the deep neural classifier (g), which determines whether the flow behaves regularly or leaks. Activation mappings in (h) show regions on which the network is basing its classification. As expected, the activation mapping is diffusive in case of regular CO~2~ plume development and focused on the leakage location when CO~2~ plume behaves irregularly.

### Proxy seismic and fluid-flow models 

Amongst the various CO~2~ injection projects, GCS in offshore saline aquifers has been most successful in reaching scale and in meeting injection targets. For that reason, we consider a proxy model constructed representative for CO~2~ injection in the South of the North Sea involving a saline aquifer made of the highly permeable Blunt sandstone. This area, which is actively being considered for GCS [@kolster2018impact], consists of the following three geologic sections (see Figure #fig:perm for the permeability and porosity distribution):

(i) the highly porous (average ``33\%``) and permeable (``>170\,\mathrm{mD}``) Blunt sandstone reservoir of about ``300-500\,\mathrm{m}`` thick. This section, denoted by red colors in Figure #fig:perm, corresponds to the saline aquifer and serves as the reservoir for CO~2~ injection;

(ii) the primary seal (permeability ``10^{-4}-10^{-2}\,\mathrm{mD}``) made of the Rot Halite Member, which is ``50\,\mathrm{m}`` thick and continuous (black layer in Figure #fig:perm);

(iii) the secondary seal made of the Haisborough group, which is ``>300\,\mathrm{m}`` thick and consists of low-permeable (permeability ``15-18\,\mathrm{mD}``) mudstone (purple section in Figure #fig:perm).

To arrive at the fluid-flow models, we consider 2D subsets of the 3D Compass model [@jones2012building] and convert these seismic models to fluid-flow properties (see Figure #fig:workflow (b)) by assuming a linear relationship between compressional wavespeed and permeability in each stratigraphic section. For further details on the conversion of compressional wavespeed and density to permeability and porosity, we refer to empirical relationships reported in @klimentos1991effects\. During conversion, an increase of ``1\mathrm{km/s}`` in compressional wavespeed is assumed to correspond to an increase of ``1.63\,\mathrm{mD}`` in permeability. From this, porosity is calculated with the Kozeny-Carman equation [@costa2006permeability] ``K = \mathbf{\phi}^3 \left(\frac{1.527}{0.0314*(1-\mathbf{\phi})}\right)^2``, where ``K`` and ``\phi`` denote permeability (``\mathrm{mD}``) and porosity (%) with constants taken from the Strategic UK CCS Storage Appraisal Project report.

### Figure: {#fig:perm}
![](figs/permeability.png){width=49% #fig:permeability}
![](figs/porosity.png){width=49% #fig:porosity}

:Permeability and porosity derived from a 2D slice of Compass model.


<!-- including injection (white ``\times``) and production well (yellow ``\bullet``). <!-- *(b)* The realistic geological stratigraphy reported in the Strategic UK CCS Storage Appraisal Project. -->

<!-- 
We start by converting the wavespeed and density of the rocks to permeability and porosity through empirical relationships [@klimentos1991effects]. There are 3 regions in the permeability: xxx, xxx and xxx. This process is shown in Figure #fig:perm from *(a)* to *(b)*\. -->

### Fluid-flow simulations

To model CO~2~ plumes that behave regularly and irregularly, the latter due to leakage, we solve the two-phase flow equations numerically[^1] for both pressure and concentration [@kailaix_2021_5528428;@li2020coupled]\. To mimic possible pressure-induced CO~2~ leakage, we increase the permeability at random distances away from the injection well within the seal from ``10^{-4}\,\mathrm{mD}`` to ``500\,\mathrm{mD}`` when the pressure exceeds ``\sim 15\,\mathrm{MPa}``. At that depth, the pressure is below the fracture gradient [@ringrose2020store]\. Since pressure-induced fractures come in different sizes, we also randomly vary the width of the pressure-induced fracture openings from ``12.5\,\mathrm{m}`` to ``62.5\,\mathrm{m}``. Examples of fluid-flow simulations without and with leakage are shown in Figure #fig:workflow (c)\.

[^1]: We used the open-source software [FwiFlow.jl](https://github.com/lidongzh/FwiFlow.jl) [@kailaix_2021_5528428;@li2020coupled] to solve the two-phase flow equations for both the pressure and concentration.

<!-- We inject CO~2~ in the speed of xxx at varying locations, and use a production well at the right-hand side of the model to take water out of this physical system. The fluid-flow simulator involves solving two-phase flow equations, which models the leading order behavior of supercritial CO~2~ replacing brine in the porous rocks. To generate realistic CO~2~ leakage examples, we select random locations of the seal as pressurized faults and open them by increasing the permeability by to ``100\,\mathrm{md}`` when the pressure at these locations increases by ``9\,\mathrm{MPa}``, lower than the fracturing pressure mentioned in @ringrose2020store\. The CO~2~ plume will move toward the fault and leaks through the seal over time. -->

### Rock-physics conversion

To monitor temporal variations in the plume's CO~2~ concentration seismically, we use the patchy saturation model [@avseth2010quantitative] to convert the CO~2~ concentration to decrease in compressional wavespeed and density. These changes are shown in Figure #fig:workflow from (d). The fact that these induced changes in the time-lapse differences in seismic properties are relatively small in spatial extent (``\sim 800\,\mathrm{m}`` for the plume and ``< 62.5 \,\mathrm{m}`` for the leakage) and amplitude (``18\%`` of the impedance perturbations with respect to the background model) calls for a time-lapse imaging modality with small normalized root-mean-square (NRMS) [@kragh2002seismic] values ``(\sim 2.5\%)``\. This NMRS value is based on the relative amplitude of impedance perturbations with respect to the background model ( ``14\%``)\.

### Time-lapse seismic simulations

To train and validate automatic detection of CO~2~ leakage from the storage complex requires the creation of realistic time-lapse datasets that contain the seismic imprint of regular as well as irregular (leakage) plume development. To this end, baseline surveys are simulated prior to CO~2~ injection for different subsets of the Compass model. Monitor surveys are simulated ``200`` days after leakage occurs to verify that potential leakage can be detected automatically early on. For regular plume development, we shoot monitor surveys for each subset at random times after CO~2~ injection. In order to strike a balance between acquisition productivity and time-lapse image quality, use is made of dense permanent acoustic monitoring at the seafloor with ``25\,\mathrm{m}`` receiver spacing. Time-lapse acquisition costs are reduced by non-replicated coarse shooting with the source towed at ``10\,\mathrm{m}`` depth below the ocean surface. Subsampling artifacts are reduced by using a randomized technique from compressive sensing where ``32`` sources are located at non-replicated jittered [@herrmann2008GJInps] source positions, yielding an average source sampling of ``125\,\mathrm{m}``. Given this acquisition geometry, linear data is generated[^2] with Equation #eq-lin-model for a ``25\,\mathrm{Hz}`` Ricker wavelet and with the band-limited noise term set so that the data's signal-to-noise ratio (SNR) is ``8.0\,\mathrm{dB}``\. This noise level leads to an extremely poor SNR of ``-31.4\,\mathrm{dB}`` for time-lapse differences in the shot domain. See Figure #fig:workflow (e).

[^2]: We used the open-source software [JUDI.jl](https://github.com/slimgroup/JUDI.jl) [@witte2018alf;@mathias_louboutin_2022_7086719] to model the wave propagation. This Julia package calls the highly optimized propagators of [Devito](https://www.devitoproject.org/) [@louboutin2018dae;@luporini2020architecture;@fabio_luporini_2022_6958070].

<!-- We collect a baseline survey before CO~2~ injection. To verify the detectability of leakage at early period, we shoot a monitor survey after two years of CO~2~ leakage over the seal. We use relatively dense ocean bottom nodes that permanently sit at the ocean bottom with ``25\,\mathrm{m}`` spacing. To keep the operational costs down, we shoot non-replicated source at ``10\,\mathrm{m}`` depth below the ocean surface. We jitter sample [@herrmann2008GJInps] ``16`` source locations for baseline survey and another ``16`` locations of monitor survey. We use the open-source software [JUDI.jl](https://github.com/slimgroup/JUDI.jl) [@witte2018alf;@@mathias_louboutin_2022_7023581] to model the wave propagation, which calls the highly optimized propagators of [Devito](https://www.devitoproject.org/) [@louboutin2018dae;@luporini2020architecture;@fabio_luporini_2022_6958070]. We generate linear data via born modeling on acoustic impedance perturbations of baseline and monitor surveys. In practice, we can obtain the linearized data for each vintage by subtracting wavefields modeled in the background velocity model from the observed data from the field. To mimic the challenging scenario in practice, we add non-replicated band-limited noise (xx dB) to the baseline and monitor data. This process is shown in Figure #fig:perm from *(d)* to *(e)*\. -->

### Imaging with joint recovery model versus reverse-time migration

Given the simulated time-lapse datasets with and without leakage, time-lapse difference images are created according to two different imaging scenarios, namely via independent reverse-time migration (RTM), conducted on the baseline and monitor surveys separately, and via inversion of the joint recovery model (cf. Equations #eq-jrm and #eq-elastic). To limit the computational cost of the Bregman iterations (Equation #LBk), four shot records are selected per iteration at random from each survey for imaging [@yin2008bregman;@witte2018cls;@yang2020tdsp;@yin2021SEGcts], limiting the cost of the joint inversion to three RTMs. The recovered baseline images are shown in Figures #fig:RTM for RTM and #fig:JRM for JRM. For the leakage scenario, the time-lapse differences are plotted in Figures #fig:diff_RTM and #fig:diff_JRM\, for RTM and JRM respectively. For the regular plume, the time-lapse differences are plotted in Figures #fig:diff_RTM_noleak and #fig:diff_JRM_noleak, for RTM and JRM respectively. From these images, it is clear that joint inversion leads to relatively artifact-free recovery of the vintages and time-lapse differences. This observation is reflected in the NRMS values, which improve considerably as shown by the histograms in Figure #fig-difference_hist_RTM_vs_JRM\ for ``1000`` imaging experiments. Not only do the NRMS values shift towards the left, their values are also more concentrated when inverting time-lapse data with the joint recovery model. Both features are beneficial to automatic leakage detection.

### Figure: {#fig-difference_RTM_vs_JRM}
![](figs/RTM.png){width=49% #fig:RTM}
![](figs/JRM.png){width=49% #fig:JRM}\
![](figs/diff_RTM.png){width=49% #fig:diff_RTM}
![](figs/diff_JRM.png){width=49% #fig:diff_JRM}\
![](figs/diff_RTM_noleak.png){width=49% #fig:diff_RTM_noleak}
![](figs/diff_JRM_noleak.png){width=49% #fig:diff_JRM_noleak}\

:Reverse-time migration (RTM) versus inversion joint recovery model (JRM). (a) RTM image of the baseline; (b) JRM image of the baseline; (c) time-lapse difference and CO~2~ plume for independent RTM images with leakage; (d) time-lapse obtained by inverting the time-lapse data jointly with leakage; (e) time-lapse difference and CO~2~ plume for independent RTM images without leakage; (f) time-lapse obtained by inverting the time-lapse data jointly without leakage. Notice improvement in the time-lapse image quality. This improvement in reflected in the NRMS values that decrease from ``8.48\,\%`` for RTM to ``3.20\, \%`` for JRM.

### Figure: {#fig-difference_hist_RTM_vs_JRM}
![](figs/nrms.png){width=98% #fig:nrms}

:NRMS values for ``1000`` time-lapse experiments.

<!-- from ``8.48\,\%`` for RTM to ``3.20\, \%`` for JRM. -->

<!-- We then image the baseline and monitor data with two schemes, one being reverse-time migration and the other one being joint recovery. The imaging process is shown in Figure #fig:perm from *(e)* to *(f)*\. To recover high-fidelity images and suppress the low-frequency artifacts, we use the inverse scattering imaging condition [@t2012linearized;@witte2017EAGEspl] in our time-lapse imaging scheme. We first show the reverse-time migration results for the baseline and monitor, and the time-lapse difference. We can see that there are strong coherent artifacts in the image due to the non-replicated acquisitions.

Next, we jointly recover the common and innovation components in Equation #jrm via sparsity-promoting linearized Bregman algorithm [@yin2008bregman;@witte2018cls;@yang2020tdsp;@yin2021SEGcts]. To reap the benefit of stochastic optimization, we randomly select two sources for each vintage during each Bregman iteration for imaging and run three datapass. The recovered baseline and monitor images are shown in xxx, with the time-lapse difference shown in xxx.

We use the normalized root mean square (NRMS) values [@kragh2002seismic] to measure the repeatability of the recovery results, where a lower NRMS value indicates a more repeatable recovery results between baseline and monitor surveys.  -->

## Deep neural network classifier for CO~2~ leakage detection

The injection of supercritical CO~2~ into the storage complex perturbs the physical, chemical and thermal environment of the reservoir [@newell2019overview]. Because CO~2~ injection increases the pressure, this process may trigger CO~2~ leakage across the seal when the pressure increase induces opening of pre-existing faults or fractures zones [@pruess2006co2;@ringrose2020store]. To ensure safe operations of CO~2~ storage, we develop a quantitative leakage detection tool based on a deep neural classifier. This classifier is trained on time-lapse images that contain the imprint of CO~2~ plumes that behave regularly and irregularly. In case of irregular flow, CO~2~ escapes the storage complex through a pressure induced opening in the seal, which causes a localized increase in permeability (shown in Figure #fig:diff_JRM).

Because time-lapse differences are small in amplitude, and strongly localized laterally when leakage occurs, highly sensitive learned classifiers are needed. For this purpose, we follow @erdinc2022AAAIdcc and adopt the Vision Transformer (ViT) [@dosovitskiy2020image]. This state-of-the-art classifier originated from the field of natural language processing (NLP) [@vaswani2017attention]. Thanks to their attention mechanism, ViTs have been shown to achieve superior performance on image classification tasks where image patches are considered as word tokens by the transformer network. As a result, ViTs have much less image-specific inductive bias compared to convolutional neural networks [@dosovitskiy2020image]\. 

To arrive at a practical and performant ViT classifier, we start from a ViT that is pre-trained on image tasks with ``16\times 16`` patches and apply transfer learning [@yosinski2014transferable] to fine-tune this network on ``1576`` labeled time-lapse images. Catastrophic forgetting is avoided by freezing the initial layers, which are responsible for feature extraction, during the initial training. After the initial training of the last dense layers, all network weights are updated for several epochs while keeping the learning rate small. The labeled (regular vs. irregular flow) training set itself consists of ``1576`` time-lapse datasets divided equally between regular and irregular flow. 

After the training is completed, baseline and monitor surveys are simulated for ``394`` unseen Earth models with regular and irregular plumes. These simulated time-lapse datasets are imaged with JRM by inverting the matrix in Equation #eq-jrm via Bregman iterations in Equation #LBk\. The resulting time-lapse difference images (see Figures #fig:diff_JRM and #fig:diff_JRM_noleak for two examples) serve as input to the ViT classifier. Refer to Figure #fig-confusion for performance, which corresponds to a two by two confusion matrix. The first row denotes the classification results for samples with regular plume (negative samples), where ``193`` (true negative) out of ``206`` samples are classified correctly. The second row denotes the classification results for samples with CO~2~ leakage over the seal (positive samples), where ``147`` (true positive) out of ``188`` samples are classified correctly. Due to the fact that JRM recovers relatively artifact-free time-lapse differences, the classifier does not pick up too many artifacts related to finite acquisition as CO~2~ leakage. This leads to much fewer false alarms for CO~2~ leakage.

### Figure: {#fig-confusion}
![](figs/confusion.png){width=98% #fig:confusion}

:Confusion matrix for classifier trained on recovery images from JRM.

## Class activation mapping based saliency map

While our ViT classifier is capable of achieving good performance (see Figure #fig-confusion), making intervention decisions during GCS projects calls for interpretability and trustworthiness of our classifier [@hooker2019benchmark;@zhang2021survey;@mackowiak2021generative]. To enhance these features, we take advantage of class activation mappings (CAM) [@zhou2016learning]. These saliency maps help us to identify the discriminative spatial regions in each image that support a particular class decision. In our application, these regions correspond to areas where the classifier deems the CO~2~ plume to behave irregularly (if the classification result is leakage). By overlaying time-lapse difference images with these maps, interpretation is facilitated, assisting practitioners to make decisions on how to proceed with GCS projects and take associated actions. Figure #fig-cam illustrates how the Score CAM approach [@wang2020score] serves this purpose[^3]. Figure #fig:cam-leak shows the CAM result for a time-lapse difference image classified as a CO~2~ leakage (in Figure #fig:diff_JRM). Despite few artifacts around the image, the CAM clearly focuses on the CO~2~ leakage over the seal, which could potentially alert the practitioners of GCS. When the plume is detected as growing regularly, the CAM result is diffusive (shown in Figure #fig:cam-noleak). This shows that the classification decision is based on the entire image and not only at the plume area. The scripts to reproduce the experiments are available on the SLIM GitHub page [https://github.com/slimgroup/GCS-CAM](https://github.com/slimgroup/GCS-CAM).

[^3]: We used the open-source software [PyTorch library for CAM methods](https://github.com/jacobgil/pytorch-grad-cam) [@jacobgilpytorchcam] to calculate the CAM images.

### Figure: {#fig-cam}
![](figs/cam-leak.png){width=49% #fig:cam-leak}
![](figs/cam-noleak.png){width=49% #fig:cam-noleak}
:CAM for time-lapse difference images with a leaking plume and with a regular plume.

<!-- There is an extensive literature on the interpretability and trustworthiness of deep neural networks [@hooker2019benchmark;@zhang2021survey;@mackowiak2021generative]. In our approach, we take advantage of the class activation mapping (CAM) [@zhou2016learning] to visualize the discriminative regions in the seismic images, which potentially indicate the CO~2~ leakage location when anomalous CO~2~ dynamics is detected in the images. The CAM, for each classification result, is a saliency map within the same size of the image. It is computed based on the result of last layer in a network before average pooling and can be acquired without training on any bounding box information. The discriminative regions of the images are highlighted in CAM, which can be helpful for practitioners to further investigate the CO~2~ anomalous behavior and take associated actions. -->

<!-- 

### Deep neural network classifier training

The aforementioned process is repeated for ``2000`` slices of the Compass model, with half of them being "leakage" samples (CO~2~ leaks over the seal) and the other half being "regular" samples (CO~2~ plume grows regularly in the Storage Complex). We then train the ViT with time-lapse difference images, labeled with "regular" or "leakage". This training process is shown in Figure #fig:perm from *(f)* to *(g)*\. need more detail here

### Classification metrics and saliency map

Finally, we follow the Score CAM approach [@wang2020score] in order for explainability of classification results from ViT. To visualize how this saliency map informs the CO~2~ leakage, we overlay this map on the time-lapse image difference. -->

## Conclusions & discussion

By means of carefully designed time-lapse seismic experiments, we have shown that highly repeatable, high resolution and high fidelity images are achievable without insisting on replication of the baseline and monitor surveys. Because our method relies on a joint inversion methodology, it also averts labor-intensive 4D processing. Aside from establishing our claim of relaxing the need for replication empirically, through hundreds of time-lapse experiments yielding significant improvements in NMRS values, we also showed that a deep neural classifier can be trained to detect CO~2~ leakage automatically. While the classification results are encouraging, there are still false negatives. We argue that this may be acceptable since decisions to stop injection of CO~2~ are also based on other sources of information such as pressure drop at the wellhead. In future work, we plan to extend this methodology to different leakage scenarios and quantification of uncertainty. 

## Acknowledgement

We would like to thank Charles Jones and Philipp A. Witte for the constructive discussion. The CCS project information is taken from the Strategic UK CCS Storage Appraisal Project, funded by DECC, commissioned by the ETI and delivered by Pale Blue Dot Energy, Axis Well Technology and Costain. The information contains copyright information licensed under ETI Open Licence. This research was carried out with the support of Georgia Research Alliance and partners of the ML4Seismic Center.