---
title: "Advancing Geological Carbon Storage Monitoring with 3D Digital Shadow Technology"
author:
  - name: Abhinav Prakash Gahlot
    affiliations:
      name: Georgia Institute of Technology
  - name: Rafael Orozco
    affiliations:
      name: Georgia Institute of Technology
  - name: Felix J. Herrmann
    affiliations:
      name: Georgia Institute of Technology
pdf-engine: pdflatex
bibliography: abstract.bib
crossref:
  fig-prefix: figure
  eq-prefix: equation
filters:
    - SLIM
    - pseudocode
format:
  html:
    page-layout: full
    sidebar: false
    lightbox: true
    crossrefs-hover: true
  arxiv-pdf:
    keep-tex: true  
    linenumbers: false
    doublespacing: false
  pdf:
    template: EAGEabstractTemplate.latex
    csl: apa.csl
    keep-tex: true
# abstract:
---

::: {.hidden}

\newcommand{\argmin}{\mathop{\mathrm{argmin}\,}\limits}
\newcommand{\argmax}{\mathop{\mathrm{argmax}\,}\limits}
$$
\def\textsc#1{\dosc#1\csod} 
\def\dosc#1#2\csod{{\rm #1{\small #2}}} 
$$

:::

# Introduction

To limit global warming and achieve climate goals, the IPCC highlights the need for technologies capable of removing gigatonnes of CO~2~ annually (IPCC, @ipcc2018global), with Geological Carbon Storage (GCS) playing a key role in this effort. GCS captures CO~2~ emissions and stores them in deep geological formations. As a scalable solution for net-negative emissions [@ringrose2020store; @ringrose2023storage], GCS is considered essential for reducing anthropogenic CO~2~. Its safety depends on accurate monitoring of subsurface CO~2~ flow and its long-term behavior. Advanced time-lapse seismic imaging is particularly crucial for tracking CO~2~ migration, detecting leaks, optimizing operations, and ensuring the integrity of storage sites.

A Digital Shadow [@gahlot2024uads] is a framework for monitoring CO~2~ storage projects. It integrates field data from sources such as seismic and borehole well measurements to track changes in CO~2~ saturation over time. This framework updates a digital model of the CO~2~ using machine learning-assisted data assimilation techniques [@spantini2022coupling; @gahlot2024uads], including generative AI and nonlinear ensemble Bayesian filtering. By employing a Bayesian approach, it addresses uncertainties in reservoir properties, such as permeability, and incorporates this lack of information into uncertainties of plume predictions.

Compared to 2D monitoring, 3D monitoring is critical for accurately assessing and controlling GCS projects, as it provides a detailed, spatially accurate representation of induced changes in the subsurface. In GCS projects, CO~2~ injection and plume migration occur in three dimensions, making 3D seismic imaging and reservoir modeling essential for capturing the full extent of CO~2~ behavior. This abstract extends the uncertainty-aware 2D Digital Shadow framework [@gahlot2024uads] by incorporating 3D seismic imaging and 3D reservoir modeling to improve decision-making and mitigate risks associated with subsurface storage.

# Methodology
To develop an uncertainty-aware 3D digital shadow, we follow the formulation in @gahlot2024uads for the plume dynamics and observation models,

$$
\begin{aligned}
\mathbf{x}_{k} & = \mathcal{M}_k\bigl(\mathbf{x}_{k-1}, \boldsymbol{\kappa}_k; t_{k-1}, t_k\bigr)\\
               & = \mathcal{M}_k\bigl(\mathbf{x}_{k-1}, \boldsymbol{\kappa}_k\bigr),\ \boldsymbol{\kappa}_k\sim p(\boldsymbol{\kappa}) \quad \text{for}\quad k=1\cdots K.
\end{aligned}
$$ {#eq-dynamics}

where $\mathbf{x}_k$ is the time-varying spatial distribution for the CO~2~ saturation and pressure perturbation measured at time instants $k=1\cdots K$, $\mathcal{M}_k$ represents multi-phase fluid-flow simulations that transition the state from the previous timestep, $t_{k-1}$, to the next timestep, $t_k$, for all $K$ time-lapse timesteps. The symbol $\boldsymbol{\kappa}$ represents the spatial distribution of the permeability, which is highly heterogeneous [@ringrose2020store; @ringrose2023storage], and is assumed to be known only in terms of its statistical distribution $p(\boldsymbol{\kappa})$. To incorporate this lack of knowledge, every time the dynamics operator is applied new samples of the subsurface permeability distribution, $\boldsymbol{\kappa}\sim p(\boldsymbol{\kappa})$ are sampled while the CO~2~ saturation and reservoir pressure perturbations, collected at the previous timestep $k-1$ in the vector, $\mathbf{x}_{k-1}$, are advanced to the next timestep, $\mathbf{x}_{k}$. 

While the dynamics can be accurately modeled through reservoir simulations, relying solely on such simulations is impractical due to the inherent stochasticity of the permeability, which remains poorly constrained. To enhance estimates of CO~2~ plume behavior, real-time monitoring is essential. Time-lapse data from 3D surface seismic measurements [@lumley20104d] provide critical information, which we incorporate into our characterization of the CO~2~ plumes. These time-lapse datasets are integrated within a framework designed to better capture the state of CO~2~ plume, accounting for the uncertainties in the subsurface reservoir properties. The seismic datasets are modeled as

<!-- $$
\mathbf{y}_k = \mathcal{H}_k(\mathbf{x}_k,\boldsymbol{\epsilon}_k, \mathbf{\bar{m}}_k),\ \boldsymbol{\epsilon}_k \sim p(\boldsymbol{\epsilon}) \quad\text{for}\quad k=1,2,\cdots,K
$$ {#eq-obs}

where $\mathbf{y}_k$ represents the observed $3$D seismic data at timestep $k$, $\boldsymbol{\epsilon}_k$ is the colored gaussian noise added to the seismic shot records. The vector, $\mathbf{\bar{m}_k}=(\mathbf{\bar{v}}, \mathbf{\bar{\rho}})$ with $\mathbf{\bar{v}}$ and $\mathbf{\bar{\rho}}$ the reference velocity and density, represents a time-dependent seismic reference model which is used to derive the migration-velocity model for imaging. -->

$$
\mathbf{y}_k = \mathcal{H}_k(\mathbf{x}_k) + \boldsymbol{\epsilon}_k , \ \boldsymbol{\epsilon}_k \sim p(\boldsymbol{\epsilon}) \quad\text{for}\quad k=1,2,\cdots,K
$$ {#eq-obs}

where $\mathbf{y}_k$ represents the 3D seismic data at timestep $k$ and $\boldsymbol{\epsilon}_k$ is colored Gaussian noise added to the seismic shot records prior to reverse-time migration. To assimilate the fluid-flow simulations and the time-lapse seismic images, we use a nonlinear Bayesian filtering framework, using neural posterior density estimation to approximate the posterior distribution of the CO~2~ plume state conditioned on the time-lapse seismic data. This approach leverages advanced techniques from sequential Bayesian inference and amortized neural posterior density estimation [@nf] in a principled manner. To estimate the posterior distribution of the CO~2~ plume given seismic observations, we use conditional Normalizing Flows (CNFs) [@gahlot2024uads]. In this simulation-based framework, we carry out simulations of the state dynamics and observations (cf. equations [-@eq-dynamics] and [-@eq-obs]) and train the CNFs on these pairs. 

# Synthetic case study
To test our method, we use a synthetic 3D Earth model from the Compass model [@BG] representative for the North Sea. We select a subset of the Compass model that includes subsurface structures suitable for CO~2~ injection. The studied subset spans a physical size of 3.2 km³, discretized into a grid of dimensions $128 \times 128 \times 128$. To initialize the probabilistic digital shadow, we require an initial ensemble of potential plume scenarios. Since the reservoir properties of the CO~2~ storage sites are unknown, we establish a probabilistic baseline for the 3D permeability properties of the site and use this distribution to create the plume initialization. This permeability distribution is derived by first performing probabilistic 3D full-waveform inference [@orozco2024aspire] on the site, assuming a baseline seismic survey was conducted prior to injection. The resulting distribution of the unknown acoustic velocity models is then converted into $128$ permeability samples using the empirical relationship described in @gahlot2024uads.

## Multi-phase flow simulations 
The flow simulations are carried out using an open-source tool, JutulDarcy [JutulDarcy.jl](https://github.com/sintefmath/JutulDarcy.jl) [@jutuldarcy]. Initially, the reservoir is filled with brine, and then supercritical CO~2~ is injected at the rate chosen such that the storage capacity of the reservoir rocks does not exceed $3\%$. The injection depth is approximately $1130 m$. The simulation is performed over a 4-year interval for each time-step, $t_k$. It outputs predicted CO~2~ saturation for each of the $N=128$ ensemble members. One of the training (predicted) saturations is shown in Figure 1(c). 

## Seismic simulations
The output from the 128 flow simulations are converted to changes in the subsurface acoustic properties using the Patchy saturation model [@gahlot2024uads]. To monitor these acoustic changes, we shoot 4D seismic surveys using 16000 receivers and 16 sources with a dominant frequency of 24 Hz and a recording length of 3 sec. We added 18 dB SNR of colored Gaussian noise to the shot records. The nonlinear wave simulations and imaging are performed using the open-source software package [JUDI.jl](https://github.com/slimgroup/JUDI.jl) [@witte2018alf; @JUDI]. 

## CNFs training
The $128$ members of the ensemble serve as training pairs. These pairs are created from the above simulations and consist of a forecasted plume and the associated 4D seismic observation. For the CNF training, we implemented a first-of-its-kind 3D conditional normalizing flow in the open-source package [InvertibleNetworks.jl](https://github.com/slimgroup/InvertibleNetworks.jl) [@orozco2023invertiblenetworks]. Crucially, the scalability of this generative network to 3D volumes was enabled by the memory-frugal invertible layers of the normalizing flow.  The network weights $\boldsymbol{\phi}$ are trained by minimizing the objective @eq-loss-CNF with $\textsc{ADAM}$ [@Kingma2014AdamAM] optimizer and for a total of $120$ epochs,

$$
\widehat{\boldsymbol{\phi}} = \argmin_{\boldsymbol{\phi}} \frac{1}{M}\sum_{m=1}^M \Biggl(\frac{\Big\|f_{\boldsymbol{\phi}}(\mathbf{x}^{(m)};\mathbf{y}^{(m)})\Big\|_2^2}{2} - \log\Bigl |\det\Bigl(\mathbf{J}^{(m)}_{f_{\boldsymbol{\phi}}}\Bigr)\Bigr |\Biggr).
$$ {#eq-loss-CNF}

where $\mathbf{J}$ represents the Jacobian of the network $f_{\theta}$ with respect to its input and $\text M$ is the number of training samples. For details on the implemenation we refer to [@gahlot2024uads].

# Results
We illustrate the performance of our method on an unseen (not part of the training) plume shown in @fig-gt. To calculate posterior statistics, we use 128 posterior samples to calculate the mean and standard deviation. @fig-posterior shows the inferred plume conditioned on the observed seismic data in @fig-obs.  @fig-error shows the error between the mean of the posterior and the ground-truth plume in @fig-gt. @fig-uq shows the uncertainty of the inferred plume. We highlight that there is a good correlation between the error and uncertainty, which is an indication that the uncertainty is well calibrated. 

::: {#fig-sup layout-nrow=3}

![4D seismic observation](./figs/rtm_plume.png){#fig-obs width=85%}

![Unknown ground truth plume](./figs/gt_plume.png){#fig-gt width=85%}

![Training plume](./figs/train_plume.png){#fig-train width=85%}

![Assimilated plume](./figs/post_plume.png){#fig-posterior width=85%}

![Estimated plume error](./figs/error_plume.png){#fig-error width=85%}

![Estimated plume uncertainty](./figs/std_plume.png){#fig-uq width=85%}

*(a)* Time-lapse seismic measurements used as conditioning input to the trained conditional Normalizing Flow. *(b)* The ground-truth unknown plume that is the target to monitor *(c) An ensemble member for the state of the plume used for training. *(d)* A posterior sample conditioned on the seismic observation *(e)* Error between the inferred posterior mean and the ground truth *(f)* The variance between posterior samples. 

:::

# Conclusions
We present a 3D Digital Shadow to monitor a 3D time-lapse CO~2~ flow using 4D seismic data. Our method introduces a novel implementation of scalable normalizing flows, which facilitates the first known application of generative modeling to realistically sized 3D volumes. The results show significant improvements in CO~2~ plume estimates when seismic 4D data was incorporated at the monitoring step, particularly in bringing the predicted plume to the correct shape and size shown in the seismic and ground truth plume. The framework also demonstrates well-calibrated uncertainty quantification, with strong correlations between inferred uncertainty and ground-truth errors. This study confirms that the Digital Shadow, integrating 3D seismic data, effectively tracks CO~2~ plume dynamics and lays the groundwork for a 3D Digital Twin to optimize CO~2~ storage operations.

# Acknowledgement
This research was carried out with the support of Georgia Research Alliance, partners of the ML4Seismic Center and in part by the US National Science Foundation grant OAC 2203821.

# References {.unnumbered}

::: {#refs}
:::
