---
title: "An uncertainty-aware Digital Shadow for underground multimodal CO~2~ storage monitoring"
author:
  - name: Abhinav Prakash Gahlot
    orcid: 0009-0008-7580-1220
    email: agahlot8@gatech.edu
    affiliations:
      name: Georgia Institute of Technology
      url: https://slim.gatech.edu/
  - name: Rafael Orozco
    orcid: 0000-0003-0917-2442
    email: rorozco@gatech.edu
    affiliations:
      name: Georgia Institute of Technology
      url: https://slim.gatech.edu/
  - name: Ziyi Yin
    orcid: "0000-0002-5024-8771"
    email: ziyi.yin@gatech.edu
    affiliations:
      name: Georgia Institute of Technology
      url: https://slim.gatech.edu/
  - name: Felix J. Herrmann
    orcid: "0000-0003-1180-2167"
    email: felix.herrmann@gatech.edu
    affiliations:
      - name: Georgia Institute of Technology
        url: https://slim.gatech.edu/
date: last-modified
pdf-engine: pdflatex
toc: true
number-sections: true
number-depth: 3
format:
  html:
    page-layout: full
    sidebar: false
    lightbox: true
    crossrefs-hover: true
  arxiv-pdf:
    keep-tex: true  
    linenumbers: false
    doublespacing: false
    runninghead: "A Preprint"
  pdf:
    documentclass: article
    geometry:
      - top=1in
      - bottom=1in
      - left=1in
      - right=1in
crossref:
  fig-prefix: figure
  eq-prefix: equation
header-includes:
  - \usepackage{gensymb}
filters:
    - SLIM
    - critic-markup
    - pseudocode
bibliography: paper.bib
appendix-style: plain
#csl: geophysical-journal-international.csl
abstract: |
  As a society, we are faced with important challenges to combat climate change. Geological Carbon Storage (GCS), during which gigatonnes of super-critical CO~2~ are stored underground, is arguably the only scalable net-negative CO~2~-emission technology that is available. While promising, subsurface complexities and heterogeneity of reservoir properties demand a systematic approach to quantify uncertainty when optimizing production and mitigating storage risks, which include assurances of Containment and Conformance of injected supercritical CO~2~. As a first step towards the design and implementation of a Digital Twin for monitoring and control of underground storage operations, a machine-learning-based data-assimilation framework is introduced and validated on carefully designed realistic numerical simulations. Because our implementation is based on Bayesian inference, but does not yet support control and decision-making, we coin our approach an uncertainty-aware Digital Shadow. To characterize the posterior distribution for the state of CO~2~ plumes (its CO~2~ concentration and pressure), conditioned on multi-modal time-lapse data, the envisioned Shadow combines techniques from Simulation-Based Inference (SBI) and Ensemble Bayesian Filtering to establish probabilistic baselines and assimilate multi-modal data for GCS problems that are challenged by large degrees of freedom, nonlinear multiphysics, non-Gaussianity, and computationally expensive to evaluate fluid-flow and seismic simulations. To enable SBI for dynamic systems, a recursive scheme is proposed where the Digital Shadow's neural networks are trained on simulated ensembles for their state and observed data (well and/or seismic). Once training is completed, the system’s state is inferred when time-lapse field data becomes available. Contrary to ensemble Kalman filtering, corrections to the predicted simulated states are not based on linear updates, but instead follow during the Analysis step of Bayesian filtering from a prior-to-posterior mapping through the latent space of a nonlinear transform. Starting from a probabilistic model for the permeability field, derived from a baseline surface-seismic survey, the proposed Digital Shadow is validated on unseen simulated ground-truth time-lapse data. In this computational study, we observe that a lack of knowledge on the permeability field can be factored into the Digital Shadow's uncertainty quantification. Our results also indicate that the highest reconstruction quality is achieved when the state of the CO~2~ plume is conditioned on both time-lapse seismic data and wellbore measurements. Despite the incomplete knowledge of the permeability field, the proposed Digital Shadow was able to accurately track the unseen physical state of the subsurface throughout the duration of a realistic CO~2~ injection project. To the best of our knowledge, this work represents the first proof-of-concept of an uncertainty-aware, in-principle scalable, Digital Shadow that captures the uncertainty arising from unknown reservoir properties and noisy observations. This framework provides a foundation for the development of a Digital Twin aimed at mitigating risks and optimizing the management of underground storage projects.
---

\newpage

::: {.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}}} 
$$
\newcommand{\To}{\textbf{to}~}
\newcommand{\Output}{\textbf{Output:}~}
\renewcommand{\Return}{\State \textbf{return}~}
\renewcommand{\Require}{\textbf{Require:}~}

:::

# Introduction

## The need for Geological Carbon Storage

Geological carbon storage (GCS) constitutes capturing carbon dioxide (CO~2~) emissions from industrial processes followed by storage in deep geological formations, such as depleted oil and gas reservoirs, deep saline aquifers, or coal seams [@page2020global]. Because of its ability to scale [@ringrose2020store; @ringrose2023storage], this net-negative CO~2~-emission technology is widely considered as one of the most promising technologies to help mitigate the effects of anthropogenic CO~2~ greenhouse gas emissions, which are major contributors to global warming [@bui2018carbon; @metz2005ipcc]. According to the Intergovernmental Panel on Climate Change (IPCC, @ipcc2018global), it is imperative to contain atmospheric warming to $<2\degree$C by the year 2100 to avoid serious consequences of global warming. To achieve this goal, scalable net-negative emission technologies that can remove tens of gigatonnes of CO~2~ per year [@masson2018global] are needed in addition to reducing emissions. These needs call for a transition towards increased use of renewable energy in combination with the global deployment of GCS technology to store significant amounts of supercritical CO~2~ into geological formations to mitigate the effects of climate change [@metz2005carbon; @orr2009onshore]. When supercritical CO~2~---captured at anthropogenic sources (fossil-fuel power plants, steel, cement, and chemical industries), is stored securely and permanently in underground geological storage reservoirs, it can be scaled to mitigate upto $20\%$ of the total CO~2~ emission problem [@Ketzer2012]. This work is a first step towards an uncertainty-aware Digital Twin aimed at optimizing underground operations while mitigating risks.

## The challenge of monitoring subsurfcae CO~2~ injection projects

Our main technical objective is directed towards monitoring CO~2~-storage projects, so that assurances can be made whether storage proceeds as expected and whether CO~2~ plumes remain within the storage complex. These assurances are referred to as Conformance and Containment by the regulators and operators of GCS [@ringrose2020store; @ringrose2023storage]. Even though the proposed methodology could be applied to other settings, our focus is on CO~2~ storage in offshore saline aquifers in a marine setup as illustrated in @fig-acquisition. Compared to storage projects on land, this choice for marine settings circumvents potential complications and risks associated with future liability and geomechanically-induced pressure changes that may trigger earthquakes [@zoback2012earthquake; @he2011risk; @raza2016integrity]. While this approach simplifies the complexity of underground CO~2~ storage projects to some degree, simulations for the state of CO~2~ plumes are on balance more challenging than conventional reservoir simulations because the fluid-flow properties are less well-known. Besides, the process of displacing brine with supercritical CO~2~ is intricate, taxing multi-phase flow simulators [@ringrose2020store; @ringrose2023storage]. To illustrate the possible impact of lack of knowledge on the reservoir properties (the permeability), @fig-acquisition also includes a juxtaposition between multi-phase flow simulations for two different realizations of permeability fields drawn from a stochastic baseline established from active-source surface seismic data [@yin2024wise]. While both realizations for the permeability are consistent with this baseline, they yield vastly different outcomes where in one case the CO~2~ plume remains well within (@fig-acquisition bottom-left) the storage complex (denoted by the white dashed line in @fig-acquisition) whereas the other simulation obviously breaches Containment (@fig-acquisition bottom-right), due to the presence of an "unknown" high-permeability streak. By maximally leveraging information derived from multimodal time-lapse monitoring data, we plan to guide uncertain reservoir simulations, so outlying "unpredicted" geological events, such as the one caused by the high-permeability streak, can be accounted for.

## Proposed approach with Digital Shadows

As stated eloquently by @ringrose2020store; @ringrose2023storage, CO~2~-storage projects can not rely on static reservoir models alone [@ringrose2016reservoir]. Instead, reservoir models and simulations need to be updated continuously from multimodal sources (e.g. active-source surface seismic and/or pressure/saturation measurements at wells) of time-lapse data. By bringing this insight together with the latest developments in generative Artificial AI, which includes the sub-discipline of Simulation-Based Inference (SBI, @cranmer2020frontier), where conditional neural networks are trained to act as posterior density estimators, and nonlinear ensemble filtering [@spantini2022coupling], a Digital Shadow [@asch2022toolbox] is coined as a framework for monitoring underground CO~2~-storage projects. Digital Shadows differ from Digital Twins [@asch2022toolbox] in the sense that they integrate information on the state of a physical object or system --- e.g. the temporal behavior of CO~2~-injection induced changes in pressure/saturation for GCS --- by capturing information unidirectionally from real-time, possibly multimodal, observed time-lapse data without feedback loops. While our framework is designed to allow for future optimization, control, and decision-making to mitigate risks of GCS that would turn our Shadow into a Twin, our main focus in this work is to introduce a scalable machine-learning assisted data-assimilation framework [@asch2022toolbox; @spantini2022coupling]. To be relevant to realistic CO~2~-storage projects, use is made of scalable state-of-the-art reservoir simulation, wave-equation based seismic inversion [@virieux2009overview; @baysal1983reverse], ensemble Bayesian filtering [@asch2022toolbox], and the advanced machine-learning technique of neural posterior density estimation [@cranmer2020frontier; @radev2020bayesflow], to capture uncertainties due to our lack of knowledge on the reservoir properties, the permeability field, in a principled way. Because the envisioned Digital Shadow is aimed at monitoring CO~2~-injection projects over many decades, its performance is validated *in silico*---i.e. via a well-designed computer-simulation study. As with the development of climate models, this approach is taken by necessity because suitable time-lapse datasets spanning over long periods of time are not yet available. 

## Relation to exisiting work

The central approach taken in this work finds its basis in ensemble Kalman Filtering (EnKF, @evensen1994sequential; @asch2016data; @carrassi2018data), a widely employed and scalable Bayesian data-assimilation (DA) technique. Here, DA refers to a statistical approach where indirect and/or sparse real-world time-lapse measurements are integrated into a model that describes the dynamics of the hidden state of CO~2~ plumes. While other time-lapse data-collection modalities exist, including controlled-source electromagnetics [@AYANI2020103098; @9186368], gravity [@doi:https://doi.org/10.1002/9781119156871.ch14], history matching from data collected at wells [@ZHANG20143000; @STRANDLI20144473], and even tomography with muons [@gluyas2019passive], we employ active-source time-lapse surface-seismic surveys [@doi:10.1190/1.1444921] in combination with data collected at wells. Seismic surveys are known for their relatively high attainable resolution, sensitivity to small CO~2~-induced changes in seismic rock properties, and spatial coverage in both the vertical and horizontal coordinate directions [@ringrose2020store; @ringrose2023storage]. 

By including crucial multi-phase fluid-flow simulations [@nordbotten2011geological] for the state, and by exchanging the extended Kalman filter [@asch2016data; @asch2022toolbox] for better scalable EnKF [@bruer2024seismicmonitoringco2plume] type of approach, our work extends early seismic monitoring work by @huang2020towards, where full-waveform inversion (FWI, @virieux2009overview) is combined with DA. Our choice for a scalable ensemble filtering type of approach [@herrmann2023president; @gahlot2023NIPSWSifp; @spantini2022coupling] is supported by recent work of @bruer2024seismicmonitoringco2plume, who conducted a careful study on the computational complexity of realistic CO~2~ projects. This work is aimed at further improving several aspects of ensemble-filter approaches including *(i)* incorporation of stochasticity in the state dynamics by treating the permeability distribution as a random variable newly drawn after each time-lapse timestep; *(ii)* removal of reliance on limiting Gaussian assumptions and linearizations of the observation operator; *(iii)* replacement of the crucial linear Kalman correction step by a learned nonlinear mapping similar to the coupling approach recently proposed by @spantini2022coupling; and *(iv)* extending the approach to multimodal time-lapse data.

To meet these additional complexities and evident nonlinearities, we develop a new neural sequential Bayesian inference framework based on Variational Inference [@jordan1998introduction] with Conditional Normalizing Flows (CNFs, @rezende2016variational; @nf; @radev2020bayesflow). As in the work by @spantini2022coupling, nonlinearities in both the state and observation operators are accounted for by a prior-to-posterior mapping through the latent space of a nonlinear transform, given by a CNF in this work. However, by allowing the permeability field, which emerges in the state for the fluid flow, to appear as a stochastic nuisance parameter, our proposed formulation sets itself apart from existing work. By using the fact that CNFs implicitly marginalize over nuisance parameters during training [@cranmer2020frontier; @radev2020bayesflow], stochasticity due to the lack of knowledge on the permeability field is factored into the uncertainty captured by our neural approximation for the posterior distribution of the CO~2~ plume's state.

Finally, the presented approach provides an alternative to "pure" machine-learning approaches, where neural networks are either trained on pairs saturation/seismic shot records [@um2023real], on reservoir models and saturation plumes [@TANG2021103488], or on "time-lapse movies" [@huang2024diffusionbasedsubsurfacemultiphysicsmonitoring], to produce estimates for the dynamic state with uncertainties based on Mont-Carlo dropout. Aside from the argument that this type of uncertainty quantification may not be Bayesian [@folgoc2021mc], our proposed approach is solidly anchored by the fields of DA, where high-fidelity simulations for the state and observations are used to recursively create forecast ensembles, and by SBI [@cranmer2020frontier; @radev2020bayesflow], where advantage is taken of the important concept of *summary statistics* [@orozco2023adjoint; @orozco2024aspireiterativeamortizedposterior] to simplify the complex relationship between the state (saturation/pressure) and seismic observations (shot records). Instead of relying on the network to provide the uncertainty, our deployment of CNFs aims to learn the Bayesian posterior distribution directly that includes both aleatoric, by adding additive noise to the observed time-lapse data, and epistemic uncertainty, by treating the permeability field as a stochastic nuisance parameter. While our approach can be extended to include network-related epistemic uncertainty [@qu2024uncertaintyquantificationseismicinversion], we leave this to future work to keep our approach relatively simple.

## Outline

The Digital Shadow for underground CO~2~-storage monitoring is aimed to demonstrate how learned nonlinear Bayesian ensemble filtering techniques can be used to characterize the state from multimodal time-lapse monitoring data. Our work is organized as follows. In @sec-problem, we start by introducing the main ingredients of nonlinear ensemble filtering, followed in @sec-method by an exposition on how techniques from amortized Bayesian inference [@cranmer2020frontier; @radev2020bayesflow] can be used to derive our nonlinear extension of ensemble filtering based on CNFs. After specifying algorithms for Sequential Simulation-based Bayesian Inference and $\textsc{Forecast}$, $\textsc{Training}$, and $\textsc{Analysis}$, which includes an introduction of the important aspect of summary statistics, we proceed by illustrating the "lifecycle" of the proposed Digital Shadow, which includes the important prior-to-posterior mapping that corrects the predictions for the state from time-lapse observations. Next, we outline a detailed procedure to validate the proposed Digital Shadow prototype in a realistic off-shore setting representative of the South-Western North Sea. In addition to detailing this setting, which includes descriptions of the probabilistic baseline, ground-truth simulations, and monitoring with the Digital Shadow, @sec-shadow also introduces certain performance metrics. Guided by these metrics, a case study is presented in @sec-case with an emphasis on examining the advantages of working with multimodal data. This case study is followed by @sec-discussion in which the Digital Shadow is put in a broader context, followed by conclusions in @sec-conclusion.

::: {#fig-acquisition layout-nrow="2"}
![](./figs/Acquisition_setup_wells.png)

![](./figs/Vastly_diff_plumes.png)

Example of a setup of Geological Carbon Storage in a marine setting. Top: Schematic diagram showing the velocity model, the seismic acquisition setup with sources (denoted by red $X$ symbol), and receivers (denoted by the yellow $\nabla$ symbol). The setup also includes injection (left) and monitoring (right) wells. Bottom: the juxtaposition of multi-phase flow simulations for two different realizations of a permeability field drawn from a stochastic baseline established from active-source surface seismic data [@yin2024wise]. Both realizations for the permeability are consistent with this baseline seismic survey but yield vastly different outcomes. Bottom-left CO~2~ plume remains well within the storage complex (denoted by the white dashed box). Bottom-right: plume obviously breaches Containment, due to the presence of an "unknown" high-permeability streak.
:::

# Problem formulation {#sec-problem}

Monitoring CO~2~ plumes from multimodal time-lapse measurements from well and/or surface-seismic data calls for an approach where nonlinear causes of epistemic uncertainty, reflecting lack of knowledge on the reservoir properties determining the dynamics, and aleatoric uncertainty, due to noise in the data, are accounted for in a sequential Bayesian inference framework [@tatsis2022sequential; @rubio2023transport]. As a lead into the definition of our Digital Shadow based on amortized neural posterior estimation [@papamakarios2016fast; @lueckmann2017flexible; @cranmer2020frontier], we begin by defining general models for nonlinear state dynamics and observations, followed by formulating a general approach to nonlinear Bayesian filtering [@asch2022toolbox].

## Uncertain state dynamics and observation models

While the presented framework is designed to be general, the focus will be on geophysical monitoring of Geological Carbon Storage (GCS) aimed at keeping track of the dynamics of CO~2~ plumes via a combination of multi-phase flow simulations and time-lapse multimodal geophysical measurements, consisting of sparse localized direct measurements of the saturation and pressure perturbations (shown in @fig-ground-truth-multimodal) collected within boreholes of (monitoring) wells and indirect active-source seismic measurements, collected at the surface. In this scenario, the time-varying state-dynamics are measured at discrete time instants, $k=1\cdots K$, and consists of the time-varying spatial distribution for the CO~2~ saturation, $\mathbf{x}_k[S_{\mathrm{CO_2}}]$, and induced pressure perturbation, $\mathbf{x}_k[\delta p]$, which can be expressed as follows: 

$$
\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}

In this expression for the uncertain nonlinear dynamics, the possibly time-dependent *dynamics operator*, $\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 dynamics operator, $\mathcal{M}_k$, which implements the underlying continuous dynamics between times $t_{k-1}$ and $t_k$, may itself be time-dependent because of induced chemical and/or rock mechanical changes. Because the spatial distribution of the permeability field, collected in the vector $\boldsymbol{\kappa}$, is highly heterogeneous [@ringrose2020store; @ringrose2023storage], unknown, and possibly time-dependent, it is treated as a stochastic nuisance parameter [@radev2020bayesflow]. Without loss of generality, we assume samples for the spatial gridded permeability field to be drawn from the time-invariant probability distribution---i.e, $\boldsymbol{\kappa}\sim p(\boldsymbol{\kappa})$ where $p(\boldsymbol{\kappa})$ is assumed to be independent of the state and time-lapse observations. The dynamics operator takes samples from the spatial distribution for the subsurface permeability, $\boldsymbol{\kappa}$, and the CO~2~ saturation and reservoir pressure perturbations, collected at the previous timestep $k-1$ in the vector, $\mathbf{x}_{k-1}$, as input and produces samples for the CO~2~ saturation and reservoir pressure perturbations at the next timestep, $\mathbf{x}_{k}$, as output. While this procedure offers access to samples from the transition density, the transition density itself can not be evaluated because of the presence of the stochastic nuisance parameter, $\boldsymbol{\kappa}$. To leading order, the nonlinear dynamics consists of solutions of Darcy's equation and mass conservation laws [@pruess2011; @https://doi.org/10.1029/2019WR027032]. Despite the fact that the dynamics can be modeled accurately with reservoir simulations [@krogstad2015mrst; @settgast2018geosx; @rasmussen2021open; @stacey2017validation], relying on these simulations alone is unfeasible because of a lack of knowledge on the permeability.

To improve estimates for the state, CO~2~ plumes are monitored by collecting discrete time-lapse well [@freifeld2009recent] and surface-seismic data [@lumley20104d]. We model these multimodal time-lapse datasets 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 the vector, $\mathbf{y}_k$, represents the observed multimodal data at timestep $k$. To add realism, the seismic data is corrupted by the noise term, $\boldsymbol{\epsilon}_k$, sampled by filtering realizations of standard zero-centered Gaussian noise with a bandwidth-limited seismic source-time function. This filtered noise is normalized to attain a preset signal-to-noise ratio before it is added to the seismic shot records prior to reverse-time migration. The shot records themselves are generated with 2D acoustic simulations with the seismic properties (velocity and density) linked to the CO~2~ saturation through the patchy saturation model[^1] [@avseth2010quantitative]. 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 from which the migration-velocity model for imaging will be derived and with respect to which the time-lapse differences in the seismic shot records will be calculated. According to these definitions, the *observation operator*, $\mathcal{H}_k$, models multimodal, direct and indirect, time-lapse monitoring data, collected at the (monitoring) wells and migrated from surface-seismic shot records. While this observation operator may also depend on time, e.g. because the acquisition geometry changes, we assume it, like the dynamics operator, to be time-invariant, and drop the time index $k$ from the notation.

[^1]: To arrive at this acquisition geometry, use was made of reciprocity---i.e., the computational experiments are calculated for $16$ sources, reducing the computational costs.

Due to sparsely sampled receiver and well locations, and the non-trivial null-space of the modeling operators, the above measurement operator is ill-posed, which renders monitoring the CO~2~ dynamics challenging, a situation that is worsened by the intrinsic stochasticity of the dynamics. To overcome these challenges, we introduce a nonlinear Bayesian filtering framework where the posterior distribution for the state, conditioned on the collected multimodal time-lapse data, is approximated by Variational Inference [@jordan1999introduction; @doi:10.1080/01621459.2017.1285773] with machine learning [@rezende2016variational; @siahkoohi2022ravi; @orozco2024aspireiterativeamortizedposterior], neural posterior density estimation to be precise. With this approach, which combines the latest techniques from sequential Bayesian inference and amortized neural posterior density estimation [@nf], we aim to quantify the uncertainty of the CO~2~ plumes in a principled manner.

## Nonlinear Bayesian filtering

Because both the dynamics and observation operators are nonlinear, we begin by defining a general probabilistic model for the state and observations in terms of the following conditional distributions:

$$
\mathbf{x}_k \sim p(\mathbf{x}_k\mid\mathbf{x}_{k-1})
$$ {#eq-trans-density} 

and

$$
\mathbf{y}_k \sim p(\mathbf{y}_k\mid\mathbf{x}_{k}),
$$ {#eq-likelihood}

for $k=0\cdots K$, representing, respectively, the *transition density*, defined by @eq-dynamics, and the data *likelihood*, given by @eq-obs. These conditional distributions hold because the above monitoring problem (cf. equations [-@eq-dynamics] and [-@eq-obs]) is Markovian, which implies that the likelihood does not depend on the time-lapse history---i.e., $\left(\mathbf{y}_k\perp\mathbf{y}_{1:k-1}\right)\mid\mathbf{x}_k$ with the shorthand, $\mathbf{y}_{1:k}= \left\{\mathbf{y}_1,\ldots,\mathbf{y}_k\right\}$, so that $p(\mathbf{y}_k\mid \mathbf{x}_k, \mathbf{y}_{1:k-1})=p(\mathbf{y}_k\mid \mathbf{x}_k)$. This orthogonality assumption holds as long as the nuisance parameters are independent with respect to the state [@rubio2023transport]. The aim of Bayesian filtering is to recursively compute the posterior distribution of the state at time, $k$, conditioned on the complete history of collected time-lapse data---i.e., compute $p(\mathbf{x}_k\mid \mathbf{y}_{1:k})$. There exists a rich literature on statistical inference for dynamical systems [@mcgoff2015statistical] aimed to accomplish this goal, which entails the following three steps [@asch2022toolbox]:

**Initialization step:** Choose the prior distribution for the state $p(\mathbf{x}_0)$.

**Prediction step:** Compute the predictive distribution for the state derived from the Chapman-Kolmogorov equation yielding [@tatsis2022sequential]

$$
\begin{aligned}
p(\mathbf{x}_k\mid \mathbf{y}_{1:k-1})&=\int p(\mathbf{x}_k\mid \mathbf{x}_{k-1})p(\mathbf{x}_{k-1}\mid \mathbf{y}_{1:k-1})\mathrm{d}\mathbf{x}_{k-1}\\
& = \mathbb{E}_{\mathbf{x}_{k-1}\sim p(\mathbf{x}_{k-1}\mid\mathbf{y}_{1:k-1})}\bigl[p(\mathbf{x}_{k}\mid\mathbf{x}_{k-1})\bigr].
\end{aligned}
$$ {#eq-prediction}

This step uses @eq-dynamics to advance samples of the distribution $p(\mathbf{x}_{k-1}\mid \mathbf{y}_{1:k-1})$ to samples of $p(\mathbf{x}_k\mid \mathbf{y}_{1:k-1})$.

**Analysis step:** Apply Bayes' rule

$$
p(\mathbf{x}_{k}|\mathbf{y}_{1:k})  \propto  p(\mathbf{y}_{k}|\mathbf{x}_{k})
p(\mathbf{x}_{k}|\mathbf{y}_{1:k-1})
$$ {#eq-Bayes}

with $p(\mathbf{y}_{k}|\mathbf{x}_{k})$ and $p(\mathbf{x}_k\mid \mathbf{y}_{1:k-1})$ acting as the "likelihood" and "prior". During this step, the posterior distribution is calculated by correcting the "prior" with the "likelihood". To this end, @eq-likelihood will be used to draw samples from the likelihood.

These three steps are general and make no assumptions on the linearity and Gaussianity of the distributions involved. The general solution of the above Bayesian filtering problem is far from trivial for the following reasons: *(i)* equations [-@eq-prediction] and [-@eq-Bayes] only permit analytical solutions in cases where the state-dynamics and observation operators are linear and Gaussian distributed; *(ii)* the presence of the nuisance permeability parameter in @eq-dynamics requires access to the marginal likelihood, which is computationally infeasible [@cranmer2020frontier; @kruse2021hint]; *(iii)* likewise, the marginalization in @eq-prediction is also not viable computationally; and finally *(iv)* computing the above conditional probabilities becomes intractable in situations where the state is high dimensional and the state-dynamics and observational operators are computationally expensive to evaluate.

To address these challenges, numerous data-assimilation techniques have been developed over the years [@asch2016data]. The ensemble Kalman filter (EnKF) [@evensen1994sequential; @carrassi2018data] features amongst the most advanced techniques and is widely used for large nonlinear state-space and observation models. During the *prediction step* of the data-assimilation recursions, EnKF propagates a set of $M$ particles from the previous timestep, $\left\{\mathbf{x}^{(1)}_{k-1},\ldots,\mathbf{x}^{(M)}_{k-1}\right\}$, to the next timestep by applying the state-space dynamics (@eq-dynamics) to each particle, followed by generating simulated observation samples, $\left\{\mathbf{y}^{(1)}_{k},\ldots,\mathbf{y}^{(M)}_{k}\right\}$, based on the observation model of @eq-obs. During the subsequent *analysis step*, EnKF uses the pairs, $\{\mathbf{x}_k^{(m)}, \mathbf{y}_k^{(m)}\}_{m=1}^M$, to calculate empirical Monte-Carlo based linear updates to the forecast ensemble. While EnKF avoids the need for linearization, calculation of gradients, and is relatively straightforward to scale, it relies on linear updates that prevent it from solving the above nonlinear and non-Gaussian Bayesian filtering problem accurately, even when increasing the ensemble size [@spantini2022coupling].

Instead of relying on linear mappings, we consider the analysis step in [@eq-Bayes] as a trainable prior-to-posterior transformation and employ techniques from amortized neural posterior density estimation with CNFs [@kruse2021hintb; @gahlot2023NIPSWSifp; @herrmann2023president] to carry out the analysis step, an approach also taken by @spantini2022coupling, who proposes a closely related but different coupling technique. Our proposed approach derives from simulation-based inference [@cranmer2020frontier; @radev2020bayesflow] with invertible CNFs, which extends EnKF by allowing for learned nonlinear prior-to-posterior density transformations. CNFs represent a class of invertible neural networks that, dependent on a condition term, convert samples from a complicated distribution, through a series of nonlinear differentiable transformations, into samples from a simpler distribution, e.g. the standard Gaussian distribution. Instead of explicit evaluations of the corresponding densities, and similar to EnKF workflows, access is only needed to samples from the "likelihood" and "prior", available through simulations with equations [-@eq-dynamics] and [-@eq-obs].

# Methodology {#sec-method}

To free ourselves from limiting assumptions, which either rely on simplifying linearity and Gaussianity assumptions or on computationally unfeasible Markov chain Monte-Carlo (McMC) based approaches, extensive use will be made of Variational Inference (VI) [@jordan1999introduction] with CNFs [@rezende2016variational; @ardizzone2019guided; @kruse2021hint]. In these approaches, scalable optimization techniques are used to match parametric representations of distributions with CNFs to unknown target distributions. During simulation-based inference (SBI) [@cranmer2020frontier; @radev2020bayesflow], the forward Kullback-Leibler (KL) divergence [@kullback1951information; @siahkoohi2023reliable] is minimized on samples from the joint distribution for the state and observations, which involves a relatively easy to evaluate quadratic objective function. In this approach, the CNFs correspond to generative neural networks that serve as conditional density estimators [@radev2020bayesflow] for the posterior distribution. After explaining the training and inference phases of SBI, the technique of sequential simulation-based inference will be introduced. This technique is designed to serve as a nonlinear extension of EnKF. To this end, algorithms for the recurrent $\textsc{Forecast}$, $\textsc{Training}$, and $\textsc{Analysis}$ stages of the proposed neural Bayesian filtering approach will be presented.

## Amortized Bayesian Inference with CNFs {#sec-amortized}

Recent developments in generative Artificial Intelligence (genAI), where neural networks are trained as surrogates for high-dimensional (conditional) distributions, are key to the development and success of our uncertainty-aware Digital Shadow, which relies on capturing the posterior distribution from simulated samples for the state and observations. While any amortized neural network implementation capable of producing samples conditioned on samples for different observations will work in our setting, and this includes conditional neural networks based on conditional score-based diffusion [@batzolis2021conditional], CNFs have our preference because they allow us to quickly sample from the posterior, $\mathbf{x}\sim p(\mathbf{x}\mid\mathbf{y})$, and conduct density calculations---i.e., evaluate the posterior density $p(\mathbf{x}\mid\mathbf{y})$ and its derivatives with respect to the states, $\mathbf{x}$, and observations, $\mathbf{y}$. This latter feature is necessary when designing Digital Twins capable of uncertainty-informed decisions [@alsing2023optimal], control [@gahlot2024digital], and Bayesian experimental design [@orozco2024probabilistic]. While non-amortized approaches may be feasible, amortization of the neural estimate of the posterior distribution is essential because it allows us to marginalize over stochastic nuisance parameters (the permeability field, $\boldsymbol{\kappa}$, in our case) while capturing its uncertainty from the ensemble of dynamic states and associated time-lapse observations.

Specifically, our goal is to capture the posterior distribution from pairs, $\{\mathbf{x}^{(m)}, \mathbf{y}^{(m)}\}_{m=1}^M$ (time-dependence was dropped for simplicity). To achieve this goal, we make use of Bayesian Variational Inference [@rezende2016variational] in combination with the ability to carry out simulations of the state dynamics and observations (cf. equations [-@eq-dynamics] and [-@eq-obs]), which are of the type [@radev2020bayesflow]

$$
\mathbf{y}=\mathcal{G}(\mathbf{x},\boldsymbol{\zeta})\quad\text{with}\quad \boldsymbol{\zeta}\sim p(\boldsymbol{\zeta}\mid\mathbf{x}),\, \mathbf{x}\sim p(\mathbf{x})
$$ {#eq-simulations}

where $\mathcal{G}$ is a generic simulator that depends on samples of the prior, $\mathbf{x}\sim p(\mathbf{x})$, and nuisance parameters, $\boldsymbol{\zeta}\sim p(\boldsymbol{\zeta}\mid\mathbf{x})$, possibly conditioned on samples of the prior. These simulated samples give us implicit access to the likelihood. Given $M$ sampled pairs, $\{\mathbf{x}^{(m)}, \mathbf{y}^{(m)}\}_{m=1}^M$, simulation-based inference with CNFs proceeds by training an invertible neural network to approximate the posterior distribution, $p(\mathbf{x}\mid\mathbf{y})$, by minimizing the forward Kullback–Leibler (KL) divergence between the true and the approximate posterior for all possible data, $\mathbf{y}$. The training objective reads [@radev2020bayesflow; @siahkoohi2023reliable]

$$
\begin{aligned}
\widehat{\boldsymbol{\phi}}&= \argmin_{\boldsymbol{\phi}}\mathbb{E}_{p(\mathbf{y})}\bigl[\mathbb{KL}\bigl(p(\mathbf{x}\mid \mathbf{y})\mid\mid p_{\boldsymbol{\phi}}(\mathbf{x}\mid \mathbf{y})\bigr)\bigr]\\
& =\argmin_{\boldsymbol{\phi}}\mathbb{E}_{p(\mathbf{y})}\Bigl[\mathbb{E}_{p(\mathbf{x}\mid\mathbf{y})}\bigl[-\log p_{\boldsymbol{\phi}}(\mathbf{x}\mid\mathbf{y})+p(\mathbf{x}\mid\mathbf{y})\bigr]\Bigr] \\
& = \argmax_{\boldsymbol{\phi}}\mathbb{E}_{p(\mathbf{x},\mathbf{y})}\bigl[\log p_{\boldsymbol{\phi}}(\mathbf{x}\mid\mathbf{y})\bigr].
\end{aligned}
$$ {#eq-KL}

To arrive at this objective involving samples from the joint distribution for $\mathbf{x}$ and $\mathbf{y}$, use was made of the law of total probability that allows us to rewrite the marginalization over observations, $\mathbf{y}$, in terms of the expectation over samples from the joint distribution. The term $p(\mathbf{x}\mid\mathbf{y})$ can be dropped because it does not depend on the parameters, $\boldsymbol{\phi}$. For more details, see @siahkoohi2023reliable. The above objective is used to approximate the true posterior distribution as accurately as possible, $p_{\boldsymbol{\phi}}(\mathbf{x}\mid\mathbf{y})\approx p(\mathbf{x}\mid\mathbf{y})$ for all possible $\mathbf{x}$ and $\mathbf{y}$, drawn from the joint distribution, $p(\mathbf{x},\mathbf{y})$. Following @radev2020bayesflow, the approximate posterior, $p_{\boldsymbol{\phi}}(\mathbf{x}\mid\mathbf{y})$ is parameterized with a conditional Invertible Neural Network (cINN, @ardizzone2019guided), $f_{\boldsymbol{\phi}}$, which conditioned on, $\mathbf{y}$, defines an invertible mapping.

With the change-of-variable rule of probability

$$
p_{\boldsymbol{\phi}}(\mathbf{x}\mid\mathbf{y})=p\bigl(\mathbf{z}=f_{\boldsymbol{\phi}}(\mathbf{x};\mathbf{y})\bigr)\Bigl |\det\Bigl(\mathbf{J}_{f_{\boldsymbol{\phi}}}\Bigr)\Bigr |
$$

with $\mathbf{J}_{f_{\boldsymbol{\phi}}}$ the Jacobian of $f_{\boldsymbol{\phi}}$ evaluated at $\mathbf{x}$ and $\mathbf{y}$, the objective of maximizing the likelihood in @eq-KL becomes

$$
\begin{aligned}
\widehat{\boldsymbol{\phi}}&=\argmax_{\boldsymbol{\phi}}\mathbb{E}_{p(\mathbf{x},\mathbf{y})}\bigl[\log p_{\boldsymbol{\phi}}(\mathbf{x}\mid\mathbf{y})\bigr]\\
&= \argmax_{\boldsymbol{\phi}}\mathbb{E}_{p(\mathbf{x},\mathbf{y})}\Bigl[\log p\bigl(\mathbf{z}=f_{\boldsymbol{\phi}}(\mathbf{x};\mathbf{y})\bigr)+\log\Bigl |\det\Bigl(\mathbf{J}_{f_{\boldsymbol{\phi}}}\Bigr)\Bigr |\Bigr].
\end{aligned}
$$ {#eq-max-likelihood}

Given the simulated training pairs, $\{\mathbf{x}^{(m)}, \mathbf{y}^{(m)}\}_{m=1}^M$, maximization of the expectations in @eq-max-likelihood can be approximated by minimization of the sum:

$$
\begin{aligned}
\widehat{\boldsymbol{\phi}}&= \argmin_{\boldsymbol{\phi}}\frac{1}{M}\sum_{m=1}^M \Bigl(-\log p\bigl(f_{\boldsymbol{\phi}}(\mathbf{x}^{(m)};\mathbf{y}^{(m)})\bigr)-\log\Bigl |\det\Bigl(\mathbf{J}^{(m)}_{f_{\boldsymbol{\phi}}}\Bigr)\Bigr |\Bigr)\\
& =\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).
\end{aligned}
$$ {#eq-loss-CNF}

During derivation of this training loss for CNFs, use was made of the property: $\mathcal{N}(0,\mathbf{I})\propto \exp(-(1/2)\|\mathbf{z}\|_2^2)$. In the ideal situation where the estimated network weights, $\widehat{\boldsymbol{\phi}}$, attain a global minimum of the objective in @eq-KL, the latent distribution will be statistically independent of the conditioning data [@radev2020bayesflow]. In that case, samples from the posterior can be drawn based on the following equivalence:

$$
\mathbf{x}\sim p_{\widehat{\boldsymbol{\phi}}}(\mathbf{x}\mid\mathbf{y}) \,\Longleftrightarrow\, f^{-1}_{\widehat{\boldsymbol{\phi}}}(\mathbf{z};\mathbf{y})\sim p(\mathbf{x}\mid\mathbf{y})\quad\text{with}\quad \mathbf{z}\sim \mathcal{N}(0,\mathbf{I}). 
$$ {#eq-ideal-sampling}

According to this equivalence, sampling from the posterior corresponds to generating samples from the multidimensional standard Gaussian distribution, $\mathbf{z}\sim \mathcal{N}(0,\mathbf{I})$, followed by computing $f^{-1}_{\widehat{\boldsymbol{\phi}}}(\mathbf{z};\mathbf{y})$. In the next section, we will employ and adapt these training and inference phases for the less-than-ideal recursive training regime of data-assimilation problems where the ensemble sizes are relatively small for which the global minimum of @eq-KL may not be attained.

::: {.content-visible when-format="html" #alg-sSBI}

``` pseudocode
#| label: alg-sSBI
#| html-indent-size: "1.2em"
#| html-comment-delimiter: "//"
#| html-line-number: true
#| html-line-number-punc: ":"
#| html-no-end: false
#| pdf-placement: "htb!"
#| pdf-line-number: true

\begin{algorithm}
\caption{Sequential Simulation-based Bayesian Inference}
\begin{algorithmic}
  \Require{Datasets $\mathbf{y}_k^{\mathrm{obs}}$, collected sequentially in the field at times $k=1,\ldots,K$}
  \Output{Estimate ensemble for the states $\{\widehat{\mathbf{x}}^{(m)}_k\}_{m=1}^M$, for timesteps $k=1,\ldots,K$}
  \For{$m = 1$ \To $M$}
    \State Sample from initial state: $\widehat{\mathbf{x}}^{(m)}_0\sim p(\mathbf{x}_0)$ 
  \EndFor
  \For{$k = 1$ \To $K$}
    \State $\{(\mathbf{x}_k^{(m)},\mathbf{y}_k^{(m)}) \}_{m=1}^{M}= $\Call{Forcast}{$ p(\boldsymbol{\kappa}), p(\boldsymbol{\epsilon}_k),\{\widehat{\mathbf{x}}_{k-1}^{(m)}\}_{m=1}^M$}
    \State $\widehat{\boldsymbol{\phi}}_k$=\Call{Training}{$\{(\mathbf{x}_k^{(m)},\mathbf{y}_k^{(m)}) \}_{m=1}^{M}$}
    \State Collect data in the field $\mathbf{y}_k^{\mathrm{obs}}$
    \State $\{\widehat{\mathbf{x}}_k^{(m)}\}_{m=1}^M= $\Call{Analysis}{$\mathbf{y}^{\mathrm{obs}}_{k}, \widehat{\boldsymbol{\phi}}_k,\{(\mathbf{x}_k^{(m)},\mathbf{y}_k^{(m)}) \}_{m=1}^{M}$}
  \EndFor
\end{algorithmic}
\end{algorithm}
```

:::

::: {.content-hidden unless-format="pdf"}

``` pseudocode
#| label: alg-sSBI
#| html-indent-size: "1.2em"
#| html-comment-delimiter: "//"
#| html-line-number: true
#| html-line-number-punc: ":"
#| html-no-end: false
#| pdf-placement: "htb!"
#| pdf-line-number: true

\begin{algorithm}
\caption{Sequential Simulation-based Bayesian Inference}
\textbf{Require:} Datasets $\mathbf{y}_k^{\mathrm{obs}}$, collected sequentially in the field at times $k=1,\ldots,K$\\
\textbf{Output:} Estimate ensemble for the states $\{\widehat{\mathbf{x}}^{(m)}_k\}_{m=1}^M$, for timesteps $k=1,\ldots,K$
\begin{algorithmic}
  \For{$m = 1$ \To $M$}
    \State Sample from initial state: $\widehat{\mathbf{x}}^{(m)}_0\sim p(\mathbf{x}_0)$ 
  \EndFor
  \For{$k = 1$ \To $K$}
    \State $\{(\mathbf{x}_k^{(m)},\mathbf{y}_k^{(m)}) \}_{m=1}^{M}= $\Call{Forcast}{$ p(\boldsymbol{\kappa}), p(\boldsymbol{\epsilon}_k),\{\widehat{\mathbf{x}}_{k-1}^{(m)}\}_{m=1}^M$}
    \State $\widehat{\boldsymbol{\phi}}_k$=\Call{Training}{$\{(\mathbf{x}_k^{(m)},\mathbf{y}_k^{(m)}) \}_{m=1}^{M}$}
    \State Collect data in the field $\mathbf{y}_k^{\mathrm{obs}}$
    \State $\{\widehat{\mathbf{x}}_k^{(m)}\}_{m=1}^M= $\Call{Analysis}{$\mathbf{y}^{\mathrm{obs}}_{k}, \widehat{\boldsymbol{\phi}}_k,\{(\mathbf{x}_k^{(m)},\mathbf{y}_k^{(m)}) \}_{m=1}^{M}$}
  \EndFor
\end{algorithmic}
\end{algorithm}
```

:::

## Neural sequential simulation-based inference {#sec-sSBI}

After randomized initialization of $M$ particles of the CO~2~ plume's initial state, $\widehat{\mathbf{x}}^{(m)}_0\sim p(\mathbf{x}_0)$ (line 2 of @alg-sSBI), training pairs for the state (@eq-dynamics) and observations (@eq-obs), $\{\mathbf{x}_k^{(m)},\mathbf{y}_k^{(m)} \}_{m=1}^{M}$, for the next time step, are predicted by simulations with the state dynamics and detailed in the $\textsc{Forecast}$ step of @alg-training). For notational convenience, we already introduced the symbol, $\widehat{\quad}$, to distinguish between predicted "digital states" and analyzed predictions (denoted by the $\widehat{\quad}$) for the state conditioned on the field observed data. During the $\textsc{Forecast}$ step of @alg-training, random samples of the stochastic permeability field and noise are drawn and used to draw samples from the transition density (line 3 of @alg-training), conditioned on analyzed samples, $\{\widehat{\mathbf{x}}_k^{(m)}\}_{m=1}^M$, from the previous time step. These predictions for the advanced state are subsequently employed to sample observations (line 4 of @alg-training), producing the simulated training ensemble, $\{\mathbf{x}_k^{(m)},\mathbf{y}_k^{(m)} \}_{m=1}^{M}$. These simulated pairs are utilized during the $\textsc{Training}$ step (line 6 of @alg-sSBI) where @eq-loss-CNF is minimized with respect to the CNF's network weights, yielding $\widehat{\boldsymbol{\phi}}_k$. Given these network weights, $M$ latent-space representations for each pair are calculated via

$$
\mathbf{z}_k^{(m)}= f_{\widehat{\boldsymbol{\phi}}_k}(\mathbf{x}_k^{(m)};\mathbf{y}_k^{(m)}),
$$ {#eq-latent}

for $m=1\cdots M$ [@spantini2022coupling]. Given observed field data, $\mathbf{y}_k^{\mathrm{obs}}$, collected at time $k$ (line 7 of @alg-sSBI), these latent variables are used to correct the state for $m=1\cdots M$ via

$$
\widehat{\mathbf{x}}_k^{(m)}= f^{-1}_{\widehat{\boldsymbol{\phi}}_k}(\mathbf{z}_k^{(m)};\mathbf{y}_k^{\mathrm{obs}}).
$$ {#eq-analyses-samples}

In this expression, $\{\widehat{\mathbf{x}}_k^{(m)}\}_{m=1}^M$, corresponds for $k\geq 1$ to analyzed states, produced by the prior-to-posterior mapping (Algorithm 2 in @spantini2022coupling). To offset the potential effects of working with a limited ensemble size, $M$, we choose to adapt this approach of producing samples from the latent space calculated from the forecast ensemble. According to @spantini2022coupling, this process mitigates the effects when a global minimum of @eq-loss-CNF is not attained, implying that $p_{\widehat{\boldsymbol{\phi}}}(\mathbf{z})\perp \mathbf{x}$ may no longer hold. In that case, the latent output distribution is no longer guaranteed to be the same for any fixed, $\mathbf{y}$ [@radev2020bayesflow]. By producing latent samples that derive directly from the predicted state-observation pairs, this issue is resolved. When conditioned on the observed field data, $\mathbf{y}_k^{\mathrm{obs}}$, these latent samples improve the inference after running the CNF in the reverse direction (cf. @eq-analyses-samples). @alg-training summarizes the three key phases, namely the $\textsc{Forecast}$, $\textsc{Training}$, and $\textsc{Analysis}$ steps, which are used during the recursions in @alg-sSBI that approximate equations [-@eq-prediction] and [-@eq-Bayes]. During these recursions, for timesteps $k>1$, samples of the posterior distribution from the previous timestep, $\{\widehat{\mathbf{x}}^{(m)}_{k-1}\}_{m=1}^M$, serve as "priors" for the next timestep, producing, $\{\widehat{\mathbf{x}}^{(m)}_k\}_{m=1}^M$, which are conditioned on the observations, $\mathbf{y}_k^{\mathrm{obs}}$. This procedure corresponds to the nonlinear prior-to-posterior of @spantini2022coupling.

``` pseudocode
#| label: alg-training
#| html-indent-size: "1.2em"
#| html-comment-delimiter: "//"
#| html-line-number: true
#| html-line-number-punc: ":"
#| html-no-end: false
#| pdf-placement: "htb!"
#| pdf-line-number: true

\begin{algorithm}
\caption{Forecast, Training, and Analysis}
\begin{algorithmic}
\Procedure{Forecast}{$ p(\boldsymbol{\kappa}), p(\boldsymbol{\epsilon}_k),\{\widehat{\mathbf{x}}_{k-1}^{(m)}\}_{m=1}^M$}
  \For{$m = 1$ \To $M$}
    \State Sample from transition density: $\mathbf{x}^{(m)}_k\sim p(\cdot\mid \widehat{\mathbf{x}}^{(m)}_{k-1})$ by running dynamics simulations (equation 1)
    \State Sample from the likelihood: $\mathbf{y}^{(m)}_k\sim p(\cdot\mid \mathbf{x}^{(m)}_{k})$ by simulating observations (equation 2)
  \EndFor 
  \Return $\{(\mathbf{x}_k^{(m)},\mathbf{y}_k^{(m)}) \}_{m=1}^{M}$
\EndProcedure 
\Procedure{Training}{$\{(\mathbf{x}_k^{(m)},\mathbf{y}_k^{(m)}) \}_{m=1}^{M}$}
  \State $\widehat{\boldsymbol{\boldsymbol{\phi}}}_k=\mathop{\mathrm{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)$ 
  \Return $\widehat{\boldsymbol{\phi}}_k$
\EndProcedure
\Procedure{Analysis}{$\mathbf{y}^{\mathrm{obs}}_{k}, \widehat{\boldsymbol{\phi}}_k,\{(\mathbf{x}_k^{(m)},\mathbf{y}_k^{(m)}) \}_{m=1}^{M}$}
  \For{$m = 1$ \To $M$}
    \State Compute latent variables: $\mathbf{z}^{(m)}_k=f_{\widehat{\boldsymbol{\phi}}_k}(\mathbf{x}_k^{(m)};\mathbf{y}_k^{(m)})$
    \State Apply prior-to-posterior map: $\widehat{\mathbf{x}}_k^{(m)}=f^{-1}_{\widehat{\boldsymbol{\phi}}_k}(\mathbf{z}_k^{(m)};\mathbf{y}_k^{\mathrm{obs}})$
  \EndFor 
  \Return $\{\widehat{\mathbf{x}}^{(m)}_k\}_{m=1}^M$
\EndProcedure
\end{algorithmic}
\end{algorithm}
```

To illustrate the different stages of the Digital Shadow's recursions, @fig-life-cycle is included to depict its "lifecycle" at the first timestep, $k=1$. <!--This figure is a snapshot of a movie included in the ancilary material, which shows the different samples.--> The lifecycle entails the following stages. Samples for the initial state of the CO~2~ saturation, denoted with a slight abuse of notation by $\widehat{\mathbf{x}}^{(m)}_0\sim p(\mathbf{x}_0)$, are drawn first from the uniform distribution, $\widehat{\mathbf{x}}^{(m)}_0\sim \mathcal{U}(0.1,\, 0.6)$ (@fig-life-cycle (top-left)). By simulating the state-dynamics @eq-dynamics (@fig-life-cycle top-row middle), predictions for the state at the next timestep are sampled, $\mathbf{x}_{k}\sim p(\mathbf{x}_{k}| \widehat{\mathbf{x}}_{k-1})$. Given these predicted samples for the state, corresponding samples for the observations are simulated with @eq-obs (@fig-life-cycle top-row right). These two types of simulations make up the $\textsc{Forecast}$ step, which produces the predicted training ensemble, $\bigl\{\mathbf{x}_{k}, \mathbf{y}_{k}\bigr\}_{m=1}^M$, consisting of $M$ samples drawn from the joint distribution, $\mathbf{y}_{k}\sim p(\mathbf{x}_{k}, \mathbf{y}_{k})$. This predicted ensemble serves as input to the $\textsc{Training}$ step (@fig-life-cycle bottom right) during which the weights of the Digital Shadow's neural networks are trained by minimizing @eq-max-likelihood. This produces the optimized network weights, $\widehat{\boldsymbol{\phi}}_k$, at timestep, $k$. After collecting field observations, $\mathbf{y}_k^{\mathrm{obs}}$, posterior inference is conducted with the neural posterior density estimator, $p_{\widehat{\boldsymbol{\phi}}_k}(\widehat{\mathbf{x}}_{k} \mid\mathbf{y}_{k}^{\mathrm{obs}})$, conditioned on the observed field data, $\mathbf{y}_k^{\mathrm{obs}}$. To avoid finite sample size and training effects during the $\textsc{Analysis}$ step (@fig-life-cycle lower left), a sample from the posterior is drawn by evaluating the CNF in the forward latent-space direction on the predicted pair, $(\mathbf{x}_{k}, \mathbf{y}_{k})$, yielding $\mathbf{z}_k = f_{\widehat{\boldsymbol{\phi}}_k}(\mathbf{x}_k;\mathbf{y}_k)$. This operation is followed by running the CNF, conditioned on the observed field data, $\mathbf{y}_{k}^{\mathrm{obs}}$, in the reverse direction, yielding the sample $\widehat{\mathbf{x}}_k= f^{-1}_{\widehat{\boldsymbol{\phi}}_k}(\mathbf{z}_k;\mathbf{y}_k^{\mathrm{obs}})$. By design, this sample approximates a sample from the true posterior, $\mathbf{x}_k\sim p(\mathbf{x}_k\mid \mathbf{y}_k^{\mathrm{obs}})$. At the next timestep, this sample becomes a sample for the "prior"---i.e., for $k\rightarrow k+1$, and the recursions proceed by repeating themselves for all, $K$, timesteps. By comparing @fig-life-cycle's sample for the plume prediction, $\mathbf{x}_{k}\sim p(\mathbf{x}_{k}| \widehat{\mathbf{x}}_{k-1})$, with the corresponding analyzed sample conditioned on the observed field data, $\widehat{\mathbf{x}}_{k} \sim p_{\widehat{\boldsymbol{\phi}}_k}(\widehat{\mathbf{x}}_{k} \mid\mathbf{y}_{k}^{\mathrm{obs}})$, it becomes apparent that the predicted state produced by the $\textsc{Forecast}$ step undergoes a significant correction during the $\textsc{Analysis}$ step. This mapping during the $\textsc{Analysis}$ step between the predicted and inferred states---i.e., $\mathbf{x}_k\mapsto \widehat{\mathbf{x}}_k$, derives from the invertibility of CNFs and the implicit conditioning of the predictions on $\mathbf{y}^{\mathrm{obs}}_{k-1}$ (via $\widehat{\mathbf{x}}_{k-1}$) and explicit conditioning on $\mathbf{y}^{\mathrm{obs}}_k$. Despite the fact that the state at the previous timestep is conditioned on $\mathbf{y}^{\mathrm{obs}}_{k-1}$, the prediction varies considerably. This variability is caused by the unknown permeability field, which appears as a stochastic nuisance parameter in the dynamics (cf. [@eq-dynamics]), and explains the significant time-lapse data-informed correction performed by the $\textsc{Analysis}$ step. The importance of the correction is nicely illustrated by the example in @fig-life-cycle where thanks to the correction the plume is extended significantly in the leftward direction, a feature completely missed by the predicted state. 

By running through the Digital Shadow's lifecycle $k$ times recursively for an ensemble consisting of $M$ particles, the $\textsc{Prediction}$ and $\textsc{Analysis}$ steps approximate samplings from $p(\mathbf{x}_k\mid \mathbf{y}^{\mathrm{obs}}_{1:k-1})$ and $p(\mathbf{x}_k\mid \mathbf{y}^{\mathrm{obs}}_{1:k})$ of equations [-@eq-prediction] and [-@eq-Bayes]. Because the proposed posterior density estimation scheme recurrently trains networks, it can be considered as an instance of neural sequential Bayesian inference, which forms the basis of our Digital Shadow. Before going into further details on how the simulations are performed, we introduce the important concept of *summary statistics*, which simplifies the observable-to-state mapping, renders typically non-uniform seismic data space uniform, and also learns to detect features that inform the posterior.

::: {#fig-life-cycle}
![](./figs/DT-life.png){width="100%"}

The life cycle of the recurrent Digital Shadow at timestep $k=1$. Starting from the initial randomly chosen CO~2~ saturation, $\widehat{\mathbf{x}}_0\sim p(\mathbf{x}_0)$, samples from the joint distribution,$\mathbf{y}_{k}\sim p(\mathbf{x}_{k}, \mathbf{y}_{k})$, are simulated during the $\textsc{Forecast}$ step, consisting of fluid-flow and observation simulations. Samples from the joint distribution form a simulated training ensemble, $\{\mathbf{x}_k^{(m)},\mathbf{y}_k^{(m)} \}_{m=1}^{M}$, consisting of $M$ training pairs that are used to train the CNF, $p_{\widehat{\boldsymbol{\phi}}_k}(\mathbf{x}_{k} |\mathbf{y}_{k})$,  which approximates the posterior. After training, the predicted plume is, during the $\textsc{Analysis}$ step, conditioned on the observed field data, $\mathbf{y}^{\mathrm{obs}}_k$, and samples from the CO~2~ plume posterior distribution, $p_{\widehat{\boldsymbol{\phi}}_k}(\widehat{\mathbf{x}}_{k} \mid\mathbf{y}_{k}^{\mathrm{obs}})$, are produced. These samples for the state are used as "priors" for the next time step. The symbol $\widehat{\quad}$ is used to distinguish between predicted "digital states" and analyzed states, conditioned on observed field data, $\mathbf{y}^{\mathrm{obs}}_k$. During the $\textsc{Analysis}$ step, the predicted states are corrected by the mapping $\mathbf{x}_k\mapsto \widehat{\mathbf{x}}_k$. 
:::

## Physics-based summary statistics

Despite the adequate expressiveness of CNF's neural architecture based on cINNs [@cnf], the mapping of CO~2~ plume-induced acoustic changes to seismic shot records remains extremely challenging because of the large degrees of freedom of Earth models and seismic shot records [@cranmer2020frontier], which live in three dimensions and five dimensions for 2- and 3D- Earth models, respectively. This challenge is compounded by the relatively weak signal strength of time-lapse changes, the complexity of the time-lapse changes to seismic observations map, and the fact that raw seismic shot records are collected non-uniformly---i.e. different seismic surveys may have different numbers of shot records or shot records may vary in the number of traces (receivers), time samples, etc. To overcome these issues, we introduce a hybrid approach, consisting of physics-based summary statistics [@radev2020bayesflow; @orozco2023adjoint], and a trainable U-net [@ronneberger2015u] that forms the conditional arm of the cINN [@ardizzone2019guided]. By mapping the shot records to the image domain, uniformity is assured because the discretization can be controlled in that domain. Moreover, the mapping from the image domain to the velocity model is much simpler, so the conditional arm of the cINN can dedicate itself to learning image-domain features that maximally inform the posterior while the network's invertible arm can concentrate on learning how to produce samples of the posterior from realizations of standard Gaussian noise conditioned on these features. More importantly, @jalsing2018massive also showed that under certain conditions the gradient of the score function (the $\log$-likelihood) compresses physical data optimally. Motivated by this work, @orozco2023adjoint demonstrated that information in the original shot records can be preserved during pre-stack reverse-time migration (RTM) as long as the migration-velocity model is sufficiently accurate. @yin2024wise showed that this condition can be relaxed in the WISE and WISER [@yin2024wiser] frameworks. The network weights for the U-net in the conditional arm are absorbed in the network weights, $\boldsymbol{\phi}$, of the CNFs, and trained during the minimization of the objective in @eq-loss-CNF.

To meet the challenges of imaging subtle time-lapse changes induced by the replacement of brine by supercritical CO~2~, we calculate the residual with respect to a reference model, $\mathbf{\bar{m}}_k$. To avoid tomographic artifacts in the gradient, the residual is migrated with the inverse-scattering imaging condition [@stolk2009inverse; @douma2010connection; @whitmore2012applications; @wangconsis2024]. This corresponds to computing sensitivities with respect to the acoustic impedance for a constant density. Given this choice, the observation operator linking the state of the CO~2~ plume to migrated time-lapse seismic images is given by

$$
\begin{aligned}
\mathbf{y}_k&=\mathcal{H}_k(\mathbf{x}_k,\boldsymbol{\epsilon}_k, \mathbf{\bar{m}}_k)\\
&:=\Bigl({\mathbf{J}_{\mathcal{F}[\mathbf{v}]}}\Bigr |_{{\mathbf{\widetilde{v}}_k}}\Bigr)_{\mathsf{ic}}^{\top} \Bigl(\mathcal{F}[\mathbf{\bar{m}}_k]-\mathbf{d}_k\Bigr)
\end{aligned}
$$ {#eq-image-observed}

In this expression, the linear operator $\Bigl({\mathbf{J}_{\mathcal{F}[\mathbf{v}]}}\Bigr |_{{\mathbf{\widetilde{v}}_k}}\Bigr)_{\mathsf{ic}}$ corresponds the Jacobian of the forward wave operator, $\mathcal{F}[\mathbf{m}]$, calculated at $\mathbf{\widetilde{v}}_k$, which is a smoothed version of the reference velocity model, $\mathbf{\bar{v}_k}$. The symbol $^\top$ denotes the adjoint. Together with the reference densities, the reference wavespeeds are collected in the seismic reference model, $\bar{\mathbf{m}}_k$. The nonlinear forward model, $\mathcal{F}$, represents the active-source experiments part of the survey. Use of the inverse-scattering imaging condition is denoted by the subscript, $\mathsf{ic}$ in the definition of the reverse-time migration operator, $\Bigl({\mathbf{J}_{\mathcal{F}[\mathbf{v}]}}\Bigr |_{{\mathbf{\widetilde{v}}_k}}\Bigr)^\top_{\mathsf{ic}}$. Finally, the vector, $\mathbf{d}_k$, collects the simulated shot records at time-lapse timestep, $k$, for all active-source experiments and is simulated with

$$
\mathbf{d}_k=\mathcal{F}[\mathbf{m}_k]+\boldsymbol{\epsilon}_k,\quad \text{where} \quad \mathbf{m}_k = \mathcal{R}({\mathbf{m}_0, \mathbf{x}_k}).
$$ {#eq-forward-seismic}

In this nonlinear forward model at timestep, $k$, the CO~2~saturation induced corrected time-lapse seismic medium properties---i.e., the wavespeed, $\mathbf{v}_k$, and the density, $\boldsymbol{\rho}_k$, are computed with respect to the baseline, $\mathbf{m}_0$, via $\mathbf{m}_k = \mathcal{R}({\mathbf{m}_0, \mathbf{x}_k})$. This nonlinear operator, $\mathcal{R}(\cdot,\cdot)$, captures the rock physics, which models how changes in the state---i.e., the CO~2~ saturation, map to changes in the seismic medium properties, $\mathbf{m}_k$. Because "observed field" seismic time-lapse data is assumed to come from the same forward model, @eq-image-observed will also be used to *in silico* simulate $\mathbf{y}_k^{\mathrm{obs}}$, simply by replacing the simulated time-lapse shot records, $\mathbf{d}_k$, by the simulated observed time-lapse seismic shot records, $\mathbf{d}_k^{\mathrm{obs}}$, that stand for data collected in the field.

Given the exposition of the amortized Bayesian inference, the neural sequential simulation-based Bayesian inference, and physics-based summary statistics, we are now in a position to setup our Digital Shadow for GCS prior to conducting our case study.

# Digital Shadow for CO~2~ storage monitoring {#sec-shadow}

To validate our Digital Shadow, an *in silico* prototype is developed and implemented for the off-shore setting of the South-Western North Sea, an area actively considered for GCS [@eti395ccs]. For this purpose, we consider the 2D cross-section of the Compass model [@BG] and the experimental setup of @fig-acquisition, which consists of sparse ocean-bottom nodes, dense sources, an injection, and a monitor well. Based on the velocity model in @fig-acquisition, a probabilistic model for reservoir properties (permeability) will be derived first, followed by the establishment of the ground truth for the *in sillico* simulated states and associated observed time-lapse seismic and well data, which serve as input to our Digital Shadow for CO~2~ storage monitoring.

## Probabilistic baseline model for the permeability field {#sec-baseline}

To mimic a realistic CO~2~ storage project where reservoir properties are not known precisely, a stochastic baseline for the permeability field is established. To this end, the velocity model depicted in @fig-acquisition is used to generate noisy baseline seismic data (signal-to-noise ratio S/N $12\mathrm{dB}$). This seismic baseline dataset serves as input to full-waveform variational Inference via Subsurface Extensions (WISE [@yin2024wise]). Given a poor 1D migration-velocity model and noisy shot records, WISE draws samples from the posterior distribution for the velocity. These samples for the wavespeed conditioned on shot data are plotted in @fig-samples-wise (a) and are converted to the permeability fields plotted in @fig-samples-wise (b). This elementwise velocity-to-permeability conversion is done with the following nonlinear elementwise mapping: 

$$
\boldsymbol{\kappa}=\begin{cases}
\boldsymbol{\zeta}\odot\exp_\circ(\mathbf{v}-4.5)\oslash \exp_\circ(4.0-4.3) \quad &\mathrm{if}\quad \mathbf{v}\succcurlyeq 4.0\\
0.01 \exp_\circ\bigl((\log_\circ (\boldsymbol{\zeta}/0.01)/0.5)\odot(\mathbf{v}-3.5)\bigr) \quad &\mathrm{if}\quad \mathbf{v}\succcurlyeq 3.5 \\
0.01\exp_\circ(\mathbf{v}-1.5)\oslash\exp_\circ(3.5-1.5) \quad &\mathrm{else.}
\end{cases}
$$ {#eq-conversion}

In this expression, the symbols $\odot$,$\oslash$, $\exp_\circ$, and $\log_\circ$, refer to elementwise multiplication, division, exponentiation, and logarithm. The conversion itself consists of elementwise exponential scaling and randomization with $\boldsymbol{\zeta}\sim\mathrm{smooth}\bigl(\mathcal{N}(1000, 1200)\bigr)$, where is $\boldsymbol{\zeta}$ a realization of the standard Normal distribution smoothed with Gaussian kernel of width $5\mathrm{m}$. The conversion is also controlled by elementwise (denoted by the elementwise greater or equal symbol $\succcurlyeq$) operations, which involve the wavespeeds. These inequalities are chosen such that the permeability values are high (between $1600 \mathrm{mD}$ and $3000 \mathrm{mD}$) for regions where the wavespeeds are above $4.0\mathrm{km/s}$. These regions correspond to clean sandstones. Areas with wavespeeds that vary within the interval $3.5-4.0\mathrm{km/s}$ stand for dirty sands with intermediate permeability values between $10 \mathrm{mD}$ and $1600 \mathrm{mD}$ while areas with lower wavespeeds correspond to shales with the lowest permeability values of $0.001 \mathrm{mD}$. 

Because of the noise and poor migration-velocity model, the samples for the baseline velocity models produced by WISE (see @fig-samples-wise (a) for a subset of these samples) exhibit strong variability. This variability amongst the samples of the baseline $\mathbf{v}\sim p(\mathbf{v}\mid \mathbf{y}_0^{\mathrm{obs}})$ produces, after conversion by @eq-conversion, highly variable samples for the permeability field, $\boldsymbol{\kappa}\sim p(\boldsymbol{\kappa})$, included in @fig-samples-wise (b). Because the permeability values are small to extremely small in the primary and secondary seals, they are denoted by solid blue-color overburden areas in @fig-samples-wise (b). These plots exhibit high variability, especially in areas where the permeable sandstones are interspersed with low-permeable shales and mudstones. This high degree of variability in the permeability field reflects our lack of knowledge of this crucial reservoir property, which complicates CO~2~ plume predictions based on reservoir simulations alone. Our Digital Shadow is aimed at demonstrating how time-lapse data can be used to condition the simulations, so the ground-truth CO~2~ plume can be estimated more accurately from monitoring data.

::: {#fig-samples-wise layout-ncol="2"}
![a.](./figs/wise-velocities.png){width="90%"}

![b.](./figs/wise-permeabilities.png){width="90%"}

Samples of the probabilistic baseline model for the permeabily field. *(a)* samples from the posterior for the wavespeeds, $\mathbf{v}_0\sim p(\mathbf{v}_0\mid \mathbf{y}_0^{\mathrm{obs}})$, produced by WISE [@yin2024wise]; *(b)* samples of the permeability field, $\boldsymbol{\kappa}\sim p(\boldsymbol{\kappa})$, derived from the velocities with @eq-conversion.
:::

Geologically, the considered CO~2~-injection site corresponds to a proxy for a storage complex in the South-Western corner of the North Sea with a reservoir made of a thick stack of relatively clean highly-permeable Bunter sandstone [@kolster2018impact]. The study area corresponds to a 2D slice selected from the synthetic Compass model [@BG], which itself was created from imaged seismic and well-log data collected in a region under active consideration for Geological Carbon Storage (GCS) by the North-Sea Transition Authority [@eti395ccs]. The seals of the storage complex comprise the Rot Halite Member, as a primary seal about $50 \mathrm{m}$ thick, situated directly on top of the Blunt sandstone, and the Haisborough Group, as the secondary seal over $300\mathrm{m}$ thick consisting of low-permeability mudstones.

The velocities of the 3D Compass model were also used to train the networks of WISE [@yin2024wise]. The reader is referred to @yin2024wise for further details on the training of these neural networks. To create a ground-truth dataset, our simulation-based workflow follows the numerical case study by @erdinc2022AAAIdcc and @yin2022TLEdgc, where realistic time-lapse seismic surveys were generated in response to CO~2~ injection in a saline aquifer.

## *In silico* ground-truth simulations {#sec-ground-truth}

The Digital Shadow will be validated against unseen ground-truth multimodal simulations. At this early stage of development, we opt for *in silico* simulations for the ground-truth CO~2~-plume dynamics and the seismic images of the induced time-lapse changes. This choice, which is not an uncommon practice in areas such as climate modeling where future data does not exist, also offers more control and the ability to validate with respect to a ground truth to which there is no access in practice. The fluid-flow simulations themselves are carried out over $24$ timesteps of $80$ days each. The time-lapse seismic images are derived from four sparsely sampled noisy ocean-bottom node surveys, collected over a period of $1920=4\times 480$ days.

### Multi-phase flow simulations

Figure [-@fig-Ground_truth] contains a schematic overview on how the ground-truth flow simulations are initialized by randomly sampling the ground-truth (denoted by the symbol $^\ast$) permeability field, $\boldsymbol{\kappa}^\ast$, and the initial state for the CO~2~ plume at timestep, $k=0$---i.e. the CO~2~ saturation is sampled, $\widehat{\mathbf{x}}^{\ast}_0\sim \mathcal{U}(0.2,\, 0.6)$. The ground-truth anisotropic permeability field, $\boldsymbol{\kappa}$, consists of $512\times 256$ gridpoints with a grid spacing of $6.25\mathrm{m}$ and has a vertical-to-horizontal permeability ratio of $\kappa_v/\kappa_h = 0.36$. The reservoir extends $100\mathrm{m}$ in the perpendicular $y$-direction. Supercritical CO~2~ is injected at a depth of $1200\mathrm{m}$ over an injection interval of $37.5 \mathrm{m}$ for the $1920$ day duration of the project with an injection rate of $0.0500\mathrm{m^3/s}$, pressurizing the reservoir from a maximum pressure of $9.7e6$ to $1.1e7\mathrm{Pa}$. The total amount of injected CO~2~ is $5.8 \mathrm{Mt}$ with $1.1. \mathrm{Mt}$ injected annually. The hydrostatic pressure at the injection depth of $1200\mathrm{m}$ is well below the depth at which CO~2~ becomes critical.

The resulting multi-phase fluid-flow simulations [@nordbotten2011geological] by the reservoir simulator [@krogstad2015mrst; @settgast2018geosx; @rasmussen2021open; @stacey2017validation; @moyner2024nonlinear] are collected in the vector, $\mathbf{x}^\ast_{1:24}$, see @fig-Ground_truth top-row for the time-lapse states, $\mathbf{x}^\ast_{1:6:24}$ ($1:6:24$ means taking each $6^{\mathrm{th}}$ sample). The simulations themselves are conducted with an implicit method implemented with the open-source software package [JutulDarcy.jl](https://github.com/sintefmath/JutulDarcy.jl) [@møyner2023]. This industry-scale reservoir modeling package simulates how the injected CO~2~ gradually displaces brine within the connected pore space of rocks in the storage complex. Initially, the reservoir rocks are filled with brine. Both brine and CO~2~ residual saturations are set to $10\%$ each. As can be seen from @fig-Ground_truth, the CO~2~ plume follows regions of high permeability, forming extensive fingers in the lateral direction. Because the density of critical CO~2~ is smaller than that of brine ($700\mathrm{kg/m^3}$ for CO~2~ versus $1000\mathrm{kg/m^3}$ for water) buoyancy effects drive the CO~2~ plume upwards.

::: {#fig-Ground_truth}
![](./figs/Ground_truth_new.png)

Schematic diagram showing ground-truth data generation using a fixed realization for the permeability field (CO~2~ saturation shown in bright overlay on top of plots for permeability field). 
:::

### Seismic simulations

To convert changes in the reservoir's state due to CO~2~ injection to modifications in the seismic properties, a nonlinear rock physics model, $\mathcal{R}$, is introduced as in @yin2022TLEdgc and @louboutin2023lmi. At each of the $K=4$ time-lapse time instances, changes in the CO~2~ saturation are related to differences in the seismic properties by this model. When ignoring velocity changes due to pressure changes and assuming a heterogeneous distribution of the CO~2~, the patchy saturation model [@avseth2010quantitative] can be used to relate increases in CO~2~ saturation to reductions in the acoustic wavespeed and density. Given baseline wavespeeds and densities gathered in the seismic medium properties vector, $\mathbf{m}_0=(\mathbf{v}_0,\boldsymbol{\rho}_0)$, the rock physics at each timestep, $k$, reads

$$
\mathbf{m}_k=\mathcal{R}(\mathbf{m}_0,\mathbf{x}_k),\, k=1\cdots K.
$$ {#eq-rock-physics}

@Fig-velrho includes an example of how changes in the CO~2~ saturation at timestep, $k=1$, are mapped to changes in the acoustic wavespeed and density. In turn, these time-lapse changes in the seismic properties lead to changes in the simulated shot records as shown in @fig-shotrecord, which depicts time-lapse seismic differences with respect to data generated with the reference model for $k=1\cdots 4$. Shot records are simulated by solving the acoustic wave equation for varying density and acoustic wavespeeds using the open-source software package [JUDI.jl](https://github.com/slimgroup/JUDI.jl), a Julia Devito Inversion framework [@witte2018alf; @JUDI], which uses highly optimized time-domain finite-difference propagators of [Devito](https://www.devitoproject.org/) [@louboutin2018dae; @luporini2020architecture].

The four ground-truth time-lapse seismic surveys are generated after $480, 960, 1440, \mathrm{and}\, 1920$ days since the start of the CO~2~ injection. There is no free surface. Throughout the fixed acquisition setup, $200$ shot records with airgun locations sampled at $16\mathrm{m}$ interval are collected sparsely by eight Ocean Bottom Nodes (OBNs) with a receiver interval of $400\mathrm{m}$, yielding a maximum offset of $3.2\mathrm{km}$. During the seismic simulations, a Ricker wavelet with a dominant frequency of $15\mathrm{Hz}$ is used, and band-limited noise with a signal-to-noise ratio (SNR) of $28.0 \mathrm{dB}$ is added to the $1800\mathrm{ms}$ recordings sampled at $4\mathrm{ms}$. Due to the complexity of the model and time-lapse changes, the time-lapse differences with respect to the reference model are noisy and have normalized root-mean-square (NRMS) values of $2.65,\, 4.60,\, 4.60$ and $4.61$ after $480, 960, 1440, \mathrm{and}\, 1920$ days.

Ground-truth time-lapse shot records, $\mathbf{d}^{\mathrm{obs}}_{1:K}$ serve as input to the time-lapse imaging scheme, detailed below, which produces the ground-truth time-lapse images, $\mathbf{y}^{\mathrm{obs}}_{1:K}$, depicted in the bottom row of @fig-Ground_truth. Despite the strong artifacts, consisting of epistemic linearization errors and aleatoric observation noise, these time-lapse seismic images display a distinct imprint of the growing CO~2~ plume. How these images are formed and used to condition estimates for the dynamic state, will be discussed next.

::: {#fig-velrho}
![](./figs/GT_Perm_TL_vel_density.png)

Seismic property changes induced by CO~2~ injection where brine is replaced by super-critial CO~2~. Left: CO~2~ saturation shown as bright overlay on top of the permeability field. Middle: time-lapse changes in compressional wavespeed and Right: in density due to CO~2~ injection.
:::


::: {#fig-shotrecord layout-ncol="4"}
![$k=1$ with $\mathrm{NRMS}=2.56$](./figs/Shot_difference1.png){width="90%"}

![$k=2$ with $\mathrm{NRMS}=4.60$](./figs/Shot_difference2.png){width="90%"}

![$k=3$ with $\mathrm{NRMS}=4.60$](./figs/Shot_difference3.png){width="90%"}

![$k=4$ with $\mathrm{NRMS}=4.61$](./figs/Shot_difference4.png){width="90%"}

Time-lapse difference seismic shot records computed with respect to data generated with the reference model. Except for the first timestep, the NRMS values remain relatively constant.
:::

## CO~2~ monitoring with the Digital Shadow

Amortized posterior inference stands at the basis of our Digital Shadow. As described in @sec-amortized, training CNFs requires the creation of predicted training ensembles during the $\textsc{Forecast}$ step, which consists of input-output pairs made out of samples from the "prior" and "likelihood". Forming these ensembles is followed by the $\textsc{Training}$ step during which the network weights are calculated. Neural inference happens at the final $\textsc{Analysis}$ step during which members of the predicted ensemble for the state are corrected by measurements from the field, $\mathbf{y}^{\mathrm{obs}}_{k}$. After exhaustive testing, the size of the ensemble was chosen to be $M=128$, a number that will be used throughout our experiments. Aside from determining the ensemble size, this number also determines how Monte-Carlo estimates the conditional mean

$$
\mathbf{\bar{x}}= \mathbb{E}_{ p(\mathbf{x}\mid\mathbf{y})} \left[ \mathbf{x}\right]\approx \frac{1}{M}\sum_{m=1}^M\mathbf{x}^{(m)}\quad 
$$ {#eq-cmean}

and conditional variance 

$$ 
\boldsymbol{\bar{\sigma}}^2 = \mathbb{E}_{ p(\mathbf{x}\mid\mathbf{y})} \Bigl[\bigl(\mathbf{x} - \mathbb{E}_{ p(\mathbf{x}\mid\mathbf{y})}   \bigr )^2\Bigr] \approx \frac{1}{M}\sum_{m=1}^M\bigl(\mathbf{x}^{(m)}-\mathbf{\bar{x}}\bigr)^2
$$ {#eq-cstdev}

are calculated.

### The Forecast step

The creation of realistic simulation-based ensembles forms a crucial part of our Digital Shadow. Details on how these training ensembles for the Digital Shadow are assembled are provided next.

#### Uncertain dynamics

The main purpose of the Digital Shadow for CO~2~ monitoring lies in its ability to handle uncertainties, in particular epistemic uncertainty due to the unknown permeability field, which throughout the simulations for the dynamics is considered as a random vector (cf. @eq-dynamics). This implies that during each evaluation of the dynamics, $\mathcal{M}$, a new sample for the permeability field is drawn, $\boldsymbol{\kappa}_k\sim p(\boldsymbol{\kappa})$, according the procedure describe in @sec-baseline.  Because the high-permeability zones differ between samples, the injection depth varies between $1200-1250\mathrm{m}$. As can be seen from @fig-k-1-K-x, the variability amongst the permeabilities, injection depth, and initial saturation, are all responsible for the variable CO~2~ plume predictions, even after a single timestep at time $k=1$. Both the CO~2~ saturation and pressure perturbations differ for each permeability realization. The predicted CO~2~ saturations have complex shapes that vary significantly. The pressure perturbations vary less, are smooth, and as expected cover a larger area.


::: {#fig-k-1-K-x}
![](./figs/Train_perm_sat_press_k1.png)

Examples of the predicted state of the CO~2~ plume at timestep $k=1$. Left column: random realizations of the permeability field. Middle column: predictions for the CO~2~ concentration. Right column: predictions for the pressure perturbations. The different realization for the permeability field changes the shape of both the predicted saturations and pressure perturbations.
:::

#### Seismic reference model {#sec-reference}

While the different realizations for the permeability field give rise to different CO~2~ saturations, revealing their time-lapse seismic signature calls for a dedicated RTM imaging procedure with migration-velocity models, $\widetilde{\mathbf{v}}_k$, that are corrected at each time-lapse timestep, and with residuals that are computed with respect to a carefully chosen reference model. To derive these migration-velocity and reference models, the forecast ensemble for the state, $\{\mathbf{x}^{(m)}_k\}_{m=1}^M$, is used to update the seismic properties with @eq-rock-physics. Given this predicted ensemble for the state, the seismic properties and simulated shot records are given by

$$
\mathbf{m}_k^{(m)} = \mathcal{R}({\mathbf{{m}}_0}, \mathbf{x}_k^{(m)})\quad \text{and}\quad \mathbf{d}_k^{(m)} = \mathcal{F}(\mathbf{m}_k^{(m)})+\boldsymbol{\epsilon}_k,\quad m=1\cdots M,
$$ {#eq-shot-data-esemble}

by virtue of @eq-forward-seismic. 

Since the time-lapse changes are small, and epistemic and aleatoric uncertainties significant, the challenge is to define for each time-lapse survey an appropriate reference velocity model. This reference model serves two purposes. First, it is used to generate reference shot records with respect to which the time-lapse differences are calculated. Second, after smoothing, it yields migration-velocity models, $\widetilde{\mathbf{v}}_k$, for each timestep.

To compute this reference model, the seismic properties for the baseline, $\mathbf{m}_0$, are updated with the rock physics model, $\mathbf{\bar{m}}_k = \mathcal{R}({\mathbf{{m}}_0}, \mathbf{\bar{x}}_k)$ where $\mathbf{\bar{x}}_k=\frac{1}{M}\sum_{m=1}^M\mathbf{x}_k^{(m)}$ is the average state at timestep $k$. From these reference models, $\mathbf{\bar{m}}_k$, the corresponding migration-velocity models, $\widetilde{\mathbf{v}}_k$, are obtained via smoothing. Using @eq-image-observed, the migrated images for the ensemble are now computed as follows:

$$
\begin{aligned}
\mathbf{y}_k^{(m)}&=\mathcal{H}(\mathbf{x}_k^{(m)},\boldsymbol{\epsilon}_k, \mathbf{\bar{m}}_k)\\
&:=\bigl({\mathbf{J}_{\mathcal{F}[\mathbf{v}]}}\Bigr |_{{\mathbf{\widetilde{v}}_k}}\bigr)_{\mathsf{ic}}^{\top} \Bigl(\mathcal{F}[\mathbf{\bar{m}}_k]-\mathbf{d}_k^{(m)}\Bigr), \quad m=1\cdots M.
\end{aligned}
$${#eq-image-observed-final}

In this expression, seismic images for each ensemble member are calculated by migrating differences between shot records computed in the reference model and shot records of the ensemble given by  @eq-shot-data-esemble with a migration-velocity model that is given by $\widetilde{\mathbf{v}}_k$. Because the surveys are replicated, the time-lapse time-index subscript in the observation operator is dropped.

With the uncertain dynamics and observation operators in place, we are after initialization (lines 1---3 of @alg-sSBI) ready to execute the $\textsc{Forecast}$ step (cf. @alg-training line 5). During this $\textsc{Forecast}$ step, tuples $\{\mathbf{x}_k^{(m)},\mathbf{y}_k^{(m)} \}_{m=1}^{M}$ are produced that serve as input to the CNF training. By replacing the simulated ensemble shot records with the ground-truth "field" shot records, $\mathbf{d}^{\mathrm{obs}}_k$, in @eq-image-observed-final, the observed images, $\mathbf{y}^{\mathrm{obs}}_k$, introduced in @sec-ground-truth are formed.

### The Training step

Given the $\textsc{Forecast}$ ensemble, $\{\mathbf{x}_k^{(m)},\mathbf{y}_k^{(m)}\}_{m=1}^{M}$, the subsequent step is to compute the optimized training weights, $\widehat{\boldsymbol{\phi}}_k$. After random initialization, the network weights are trained by minimizing the objective on line 9 of @alg-training with the optimizer $\textsc{ADAM}$ [@Kingma2014AdamAM] and for a total of $120$ epochs. As we can see from @fig-learning-curve, the plots for the $\ell_2$-norm training objective and SSIM (see Appendix B) both behave as expected during training and validation and there is no sign of overfitting. 

::: {#fig-learning-curve layout="[[40, 40]]"}
![a.](./figs/Train_val_loss_epochs.png)

![b.](./figs/SSIM_epochs.png)

Training logs. *(a)* $\ell_2$-norm loss during training and validation. *(b)* plot for SSIM. Both plots behave as expected.
:::

To measure the performance of our network, we also plotted the structural similarity index measure (SSIM, @wang2004image), which is a metric designed to mimic human perception of image quality. See Appendix B for further details.

### The Analysis step

With the networks for the amortized posterior inference trained, the predicted samples in the ensemble for the state, produced by the $\textsc{Forecast}$ step, are corrected by conditioning the posterior distribution on the imaged "field" data. In our case, the field data consists of the *in silico* ground-truth simulations outlined in @sec-ground-truth. Given these "time-lapse field datasets", $\mathbf{y}_k^{\textrm{obs}}, \, k=1\cdots K$, the predicted states of the ensemble, $\{\mathbf{x}_k^{(m)}\}_{m=1}^{M}$, are corrected during the $\textsc{Analysis}$ step, which corresponds to executing lines 14---15 of @alg-training for each timestep, $k=1\cdots K$. To illustrate the importance of the $\textsc{Analysis}$ step at timestep, $k=3$, the first two columns of @fig-analysis-samples-mean include corrections on a single instance of the predicted state, $\mathbf{x}_3\sim p(\mathbf{x}_3\mid \widehat{\mathbf{x}}_2)$ with $\widehat{\mathbf{x}}_2\sim p_{\widehat{\phi}_2}(\cdot\mid \mathbf{y}_2^{\text{obs}})$, and the corrected plume, $\widehat{\mathbf{x}}_3\sim p_{\widehat{\phi}_3}(\cdot\mid \mathbf{y}_3^{\text{obs}})$. The saturations are plotted in the top row and pressure differences in the bottom row. The correction for the CO~2~ saturation at timestep, $k=3$, is smaller than the correction at time step, $k=1$, as shown in @fig-life-cycle. The reason for this is that the early correction at the time step, $k=1$, establishes the leftward fingering feature that persists to the later timesteps. However, this does not apply to the protruding rightward CO~2~ plume features that come in at later times and which are partially missed by the saturation prediction. While less pronounced, the individual pressure perturbation prediction is also updated considerably during the $\textsc{Analysis}$ Step.

The impact of the correction also becomes clear from the right-most two columns of @fig-analysis-samples-mean where the conditional mean before, $\bar{\mathbf{x}}_3$, and after the corrections, $\bar{\widehat{\mathbf{x}}}_3$, are juxtaposed. After conditioning on the observed data, the corrected samples exhibit less variation, which results in a significantly better resolved conditional mean for the corrected CO~2~ saturation. Similarly, the pressure difference corrections yield a pressure field that is more consistent. These improvements during the $\textsc{Analysis}$ step are also apparent when comparing the errors w.r.t. the ground truth, $\text{Error}=|\mathbf{x}_3^\ast-\mathbf{\bar{x}}_3|$ vs. $\text{Error}=|\mathbf{x}_3^\ast-\mathbf{\bar{\widehat{x}}}_3|$, plotted before and after corrections during the $\textsc{Analysis}$ step in @fig-analysis-CM-errors for the CO~2~ saturations (top-row: first two columns). The effect of the correction for the pressure differences (bottom-row: first two columns) is more difficult to interpret. However, the uncertainty in terms of the conditional deviations before, $\boldsymbol{\bar{\sigma}}_3$, and after the correction, $\bar{\boldsymbol{\widehat{\sigma}}}_3$, tightens significantly for both the CO~2~ saturations (top-row: rightmost columns) and pressure differences (bottom row: rightmost columns). The plots for the pressure contain locations of the wells where pressure differences can be measured.


::: {#fig-analysis-samples-mean}
![](./figs/Analysis_step_sample_mean.png){width="95%"}

Corrections by the $\textsc{Analysis}$ step at timestep, $k=3$. Corrections on predicted samples for the saturation, $\mathbf{x}_k[S_{\mathrm{CO_2}}]$ (top-row), and pressure perturbation, $\mathbf{x}_k[\delta p]$ (bottom-row). Left two columns include individual predicted sample, $\mathbf{x}_3\sim p(\mathbf{x}_3\mid \widehat{\mathbf{x}}_2)$ with $\widehat{\mathbf{x}}_2\sim p_{\widehat{\phi}_2}(\cdot\mid \mathbf{y}_2^{\text{obs}})$ for the tuple $(\mathbf{x}_k[S_{\mathrm{CO_2}}]$, $\mathbf{x}_k[\delta p])$ and corrected sample, $\mathbf{\widehat{x}}_3=\mathrm{Analysis}\bigl(\mathbf{y}^{\mathrm{obs}}_{3}, \widehat{\boldsymbol{\phi}}_3,(\mathbf{x}_3,\mathbf{y}_3)\bigr)$. Conditional means before the $\textsc{Analysis}$$ step, $\mathbf{\bar{x}}_3$ and after the $\textsc{Analysis}$$ step, $\mathbf{\bar{\widehat{x}}}_3$, are included in the third and fourth column. Notice the impact of conditioning the state on seismic observations, $\mathbf{y}^{\mathrm{obs}}_{3}$. 
:::

::: {#fig-analysis-CM-errors}

![](./figs/Analysis_step_errvar.png){width="95%"}


Effect of the $\textsc{Analysis}$ step on predicted and corrected posterior errors w.r.t. to the ground-truth $\mathbf{x}^\ast_3$ and posterior standard deviations. Two leftmost colums contain errors, $\text{Error}=|\mathbf{x}_3^\ast-\mathbf{\bar{x}}_3|$, before, and after, $\text{Error}=|\mathbf{x}_3^\ast-\mathbf{\bar{\widehat{x}}}_3|$, correction by the $\textsc{Analysis}$ step. Two rightmost columns depict the posterior standard deviations before and after correction.
:::

While these results where the $\textsc{Analysis}$ step reduces the uncertainty are encouraging, we will introduce in the next section a set of quantitative performance metrics before considering a case study where the benefits of working with multimodal will be evaluated.

## Performance metrics {#sec-performance}

While the main goal of this paper is to present a principled workflow to assess uncertainties in the state of the CO~2~ plume, another goal of establishing a Digital Shadow for Geological CO~2~ Storage is to accelerate innovations in monitoring to ensure safe operations towards the regulators and the general public. For this reason, we follow @donoho2023data's recipe to set the stage for accelerated development in this field by

-   releasing all our data under the Creative Commons License. This includes the ground-truth seismic and reservoir models, the creation of the samples for permeability, the generation of the ground-truth state and, "observed" time-lapse data;

-   making all our code available under the open-source software MIT license. While still work in progress, the aim is to make use of our software as frictionless as possible;

-   introducing benchmark metrics by which the success of the Digital Shadow is measured quantitatively. The aim is that new methods and proposed improvements to the monitoring life-cycle can quickly be validated in a trusted testbed of benchmarks.

For now, the proposed benchmarks include the methodology presented in this paper. In a later companion paper, comparisons between results obtained by alternative data-assimulation methods will be reviewed.

### Benchmark I: reconstruction quality

Since the ground-truth state is available as part of our *in silico* simulations, quantitative assessments of the reconstruction quality can be made by computing the Structural Similarity Index Metric Error (SSIM-ERR) and Root-Mean-Square-Error (RMSE), both detailed in Appendix B.

### Benchmark II: uncertainty calibration

While the high-dimensionality and computational complexity of our inference problem prevents rigorous validation of the presented statistical analysis, quantitative establishment of how well the predicted uncertainties, $\widehat{\boldsymbol{\sigma}}_k$, and the errors, $\text{Error}_k=|\mathbf{x}_k^\ast-\mathbf{\bar{x}}_k|$, correlate is computationally feasible. To this end, we follow @orozco2024aspireiterativeamortizedposterior and introduce the calibration test developed by @guo2017calibration; @laves2020well. During this test, Euclidean distances are calculated between posterior mean estimates for the state and the ground-truth state. According to this test, uncertainty is well-calibrated when it is proportional to the error. To establish this test, deviations for each gridpoint in $\widehat{\boldsymbol{\sigma}}$ are designated into $L$ bins of equal width. The uncertainty for each bin, $B_l$, is then calculated with 

$$
\mathrm{UQ}(B_l) = \frac{1}{|B_l|} \sum_{i\in B_l}\bar{\boldsymbol{\widehat{\sigma}}}[i]
$$ 

where $\bar{\boldsymbol{\widehat{\sigma}}}[i]$ refers to the $i^{\text{th}}$ entry of the estimated deviation vector. For each bin, the mean-square error w.r.t. the ground truth is calculated with 

$$
\mathrm{MSE}(B_l) = \frac{1}{|B_l|} \sum_{i\in B_l}(\mathbf{x}^{\ast}[i] - \widehat{\mathbf{x}}[i])^2.
$$

Qualitatively, good uncertainty calibration implies that points on the crossplot between the predicted error and actual uncertainty fall along the line $y=x$. To quantify the quality of uncertainty calibration, we follow @orozco2024aspireiterativeamortizedposterior and calculate the Uncertainty Calibration Error (UCE), defined as the average absolute difference between the predicted error and actual uncertainty across all bins. A lower UCE indicates more precise calibration, which is approximately the case in @fig-Calibration for the posterior samples conditioned on both seismic and wells. We observe a trend towards poorer calibration as the monitoring time step increases. We hypothesize that as the plume grows larger, there are more areas where errors can occur, making the calibration test more difficult to pass. When comparing a new method under this benchmark, we would qualitatively evaluate the calibration plots on a per-monitoring time basis, and also in aggregate, to summarize performance across all monitoring times. As a quantitative metric, we would compare the UCE values between methods and select the one with the lowest UCE as having the best calibration.

::: {#fig-Calibration}
![](./figs/uce1.png){width="50%"}

Calibration plot of binned mean-square errors vs. posterior deviation for the plumes conditioned on both seismic and wells.
:::

# Multimodal case study {#sec-case}

The task of the proposed Digital Shadow is to produce high-fidelity estimates for the state of the CO~2~ plume using a combination of data-assimilation and machine learning techniques that use multimodal time-lapse data as input. Aside from conducting principled uncertainty quantification, the aim is also to maximally benefit from having access to two complementary time-lapse datasets, namely time-lapse pressure/saturation measurements at the two wells (injection and monitoring) and 2D active-source time-lapse seismic surveys. The aim of these *in silico* experiments is to demonstrate the possible benefits of multimodal data-assimilation over data comprised of well or seismic data only. To this end, inferences will be carried out on both time-lapse datasets separately and jointly. After briefly describing the relative ease by which the presented Digital Shadow can be extended to include data collected at (monitoring) well(s), systematic quantitative comparisons will be made between the different approaches. This includes quantitative analyses based on the metrics introduced in @sec-performance.

## Data collection, simulation, and training

While inverting multimodal data in classical (read without machine learning) inversion/inference problems has proven challenging, especially when the different datasets are scaled differently, or collected differently, extending the proposed neural data-assimilation framework is relatively straightforward as long as the time-lapse data are collected on the same grid and at the same discrete observation times. Under that assumption, observations of the saturation/pressure ($\mathbf{x}_k[S_{\mathrm{CO_2}}]$ and $\mathbf{x}_k[\delta p]$) data can be modeled by a restriction operator, collecting these quantities at gridpoints that correspond to the well locations. During training and inference, time-lapse saturation, pressure, and seismic images each serve as input to separate input channels of the conditional branch of the cINN. The saturation and pressure are both converted to all-zero images except for values of the saturation and pressure at the well locations. Based on conditioning, by seismic only, well only, or both, the cINN is trained to produce samples of the state, consisting of gridded images for saturation and pressure perturbation. These gridded quantities are produced by separate output channels of the network. The experiment setup for the ground truth and is the same as outlined in @sec-shadow.

## Experimental protocol

Experiments to investigate the possible benefits of multimodal acquisition proceed as follows:

**Initialization step:** Draw $M=128$ samples from the initial distribution for the state, $\widehat{\mathbf{x}}_0\sim p(\mathbf{x}_0)$. The initial saturation is selected from the uniform distribution on the interval $[0.1,\, 0.6]$ and chosen to occupy a cylinder with a diameter of four grid cells ($4\times 6.25\mathrm{m}$). The location of the cylinder coincides with the injection interval at the well, which extends $37.5\mathrm{m}$, vertically. The vertical location of the well's injection interval is positioned adaptively within the closest high-permeability layer in the randomly drawn permeability fields, $\boldsymbol{\kappa}\sim p(\boldsymbol{\kappa})$. These permeability fields are drawn with replacement from a set of $1000$ converted velocity models produced by WISE. The resulting injection depth varies between $1200-1250\mathrm{m}$. The initial pressure is set to the hydraulic pressure, which increases with depth. This initialization produces the initial ensemble for the state, $\{\widehat{\mathbf{x}}_0^{(m)}\}_{m=1}^{128}$.

**Forecast step:** For each timestep, the state and multimodal observations are predicted by time advancing the state ensemble, $\{\widehat{\mathbf{x}}_{k-1}^{(m)}\}_{m=1}^{128}$, to $\{\mathbf{x}_k^{(m)}\}_{m=1}^{128}$ by running the dynamics (cf. @eq-dynamics), followed by sampling from the "likelihood" (cf. equations [-@eq-obs], [-@eq-image-observed], [-@eq-forward-seismic]). This produces the predicted ensemble, $\{(\mathbf{x}_k^{(m)}, \mathbf{y}_k^{(m)})\}_{m=1}^{128}$. These predicted states consist of the saturation, $S_{CO_2}$, and pressure perturbation, $\delta\mathbf{p}$, with respect to the initial hydraulic pressure. Depending on the time-lapse data modality, the simulated observations, $\{\mathbf{y}_k^{(m)}\}_{m=1}^{128}$, consist of imaged seismic, saturation/pressure, or both.

**Training step:** The predicted ensembles, $\{\mathbf{x}_k^{(m)}\}_{m=1}^{128}$, are split into $120$ training samples and eight validation samples. After the network training is completed, the inference is carried out on one sample of the test set, which consists of the ground-truth seismic images, $\mathbf{y}_k^{\mathrm{obs}}$.

**Analyses step:** Conditioned on observed time-lapse "field" data, $\mathbf{y}^{\mathrm{obs}}_k$, the predicted states are corrected by computing the latent representation for each ensemble pair, followed by sampling from the posterior conditioned on the field observations by running the cINN in reverse on the latent variables (cf. equations [-@eq-latent] and [-@eq-analyses-samples]). This produces the corrected ensemble for the states, $\{\widehat{\mathbf{x}}_{k}^{(m)}\}_{m=1}^{128}$.

By applying the $\textsc{Forecast}$ and $\textsc{Analysis}$ steps recursively, the full history of the estimated state for each member of the ensemble is collected in the vectors $\{\widehat{\mathbf{x}}_{1:24}^{(m)}\}_{m=1}^{128}$. These complete state histories are formed by appending six dynamic simulation timesteps after each member of the ensemble is corrected during the $\textsc{Analysis}$ step. Given these complete histories for the state, we proceed by calculating point estimates, in terms of the conditional mean, $\bar{\widehat{\mathbf{x}}}_{1:24}$ and their uncertainty in terms of the conditional deviations, $\bar{\widehat{\boldsymbol{\sigma}}}_{1:24}$, for all $24$ timesteps. To quantify the data-assimilation accuracy of the inferred states, the errors with respect to the ground truth, $\mathbf{e}_{1:24}=|\mathbf{x}^\ast_{1:24}-\mathbf{\bar{\widehat{x}}}_{1:24}|$, will be calculated as well.

## Qualitative results

To establish the performance of the proposed Digital Shadow, juxtapositions are made between data-assimilations for three different data modalities, namely well only, seismic only, and well + seismic. The importance of data-assimilation compared to relying on reservoir simulations alone will also be demonstrated by comparing plots for the conditional mean ($\bar{\widehat{\mathbf{x}}}_{1:6:24}$), the errors with respect to the ground truth, $\mathbf{e}_{1:6:24}=|\mathbf{x}^\ast_{1:6:24}-\mathbf{\bar{\widehat{x}}}_{1:6:24}|$, and the conditional standard deviation, $\bar{\widehat{\boldsymbol{\sigma}}}_{1:6:24}$ for these three different modalities. These quantities are calculated from posterior samples, conditioned on the simulated ground-truth time-lapse data included in @fig-ground-truth-multimodal for $k=1\cdots 4$. This figure contains plots for the ground-truth saturation, $\mathbf{x}^\ast_k[S_{\mathrm{CO_2}}]$; pressure perturbations, $\mathbf{x}^\ast_k[\delta p]$, including locations of the wells, and plots for the corresponding noisy imaged seismic observations in the third column. Observations at wells and/or imaged surface seismic are collected in the vectors, $\mathbf{y}^{\mathrm{obs}}_k$, for $k=1\cdots 4$, and serve as input to our Digital Shadow whose task is to infer the state of the CO~2~ plume for each timestep. 

Estimates for the CO~2~ saturation states $\bar{\widehat{\mathbf{x}}}_{1:6:24}$, $\mathbf{e}_{1:6:24}=|\mathbf{x}^\ast_{1:6:24}-\mathbf{\bar{\widehat{x}}}_{1:6:24}|$, and $\bar{\widehat{\boldsymbol{\sigma}}}_{1:6:24}$ are, for each of the four scenarios (unconditioned, conditioned on well only, seismic only, well + seismic) included in the rows of figures [-@fig-inference-t1]---[-@fig-inference-t4] for $k=1\cdots 4$. The corresponding estimates for the pressure perturbations $\mathbf{x}_{1:6:24}[\delta p]$ are shown in Appendix B. The columns of these figures contain estimates for the conditional mean, errors, and uncertainty in terms of the standard deviation. From these plots, the following qualitative observations can be made. 

First, and foremost, a comparison between the unconditioned and conditioned CO~2~ plume estimates shows that conditioning on any type of time-lapse data nudges, for each timestep, the corrected fluid-flow predictions closer to the ground-truth CO~2~ saturations. As expected, the unconditioned CO~2~ plume diverges due to epistemic uncertainties induced by the unknown permeability field at each time-lapse timestep ($k=1\cdots 4$). Remember, every time the dynamic state equation (cf. @eq-dynamics) is run, a new sample for the permeability field is drawn, introducing new uncertainties into the system. This effect is also present for the CO~2~ plume estimates conditioned on time-lapse data but to a much lesser degree because the insertion of time-lapse data cancels some of the systematic uncertainty. However, including time-lapse data at the wells by itself does not move the overall shape of the inferred CO~2~ plume much closer to the ground truth. Nevertheless, when the conditioning involves seismic data, the corrections capture the overall shape of the CO~2~ saturations very well.

Second, the conditional mean estimates for the CO~2~ saturation appear smooth compared to the ground truth. This smoothing effect is due to the remaining uncertainty of the estimated plumes (see column three of figures [-@fig-inference-t1]---[-@fig-inference-t4]) and persists, and even increases somewhat, as time progresses. 

Third, the CO~2~ plume estimates at $k=1$ obtained for the wells only completely miss the important leftward protruding CO~2~-plume feature, which eventually leads to a violation of Containment. Conversely, this feature is well recovered by inferences based on time-lapse data that includes imaged seismic. Compared to seismic only, joint inference from well + seismic data shows further improvements near the wells and away from the wells, reducing the error with respect to the ground-truth CO~2~ saturation (see column two of figures [-@fig-inference-t1]---[-@fig-inference-t4]). As expected, including pressure and saturation data at the wells reduces uncertainty in the direct vicinity of the wells but there are also indications that it may have an effect away from the wells when combined with time-lapse seismic. Finally, the errors and uncertainty in terms of the conditional standard deviations correlate reasonably well when the CO~2~ plume is conditioned on both time-lapse seismic and well data.

Fourth, the overall trend continues at timestep, $k=2$, leading to reduced errors and uncertainty in the main body of the plume when the CO~2~-plume inference is conditioned on well + seismic data. However, because of a lack of lateral information the results based on well data alone more or less completely miss the important leftward-protruding streamline feature induced by the unknown high-permeability streak.

Fifth, the complexity of the CO~2~ plume increases further at later timesteps. This challenges the inference. However, we note again that the estimated conditional mean improves considerably when the inference is conditioned on time-lapse seismic data. As can be observed in @fig-inference-t3 and [-@fig-inference-t3] bottom first and third column, the inference improves due to the inclusion of information at the wells. As a result, uncertainty is reduced in regions away from the wells both horizontally, and vertically below the wellhead. Interestingly, while the well-only inference result lacks detail, its uncertainty at the horizontal extremities seems to suggest protruding fingers at the periphery of the CO~2~ plume. These features are indeed present in the ground-truth CO~2~ plume and in the inferences that include time-lapse seismic data.

Sixth, the conditional mean estimates for the CO~2~ plume are for most features best recovered when the inference is conditioned on time-lapse seismic data with visible improvements when pressure and saturation data collected at the wells are included. When comparing the uncertainty (column three of figures [-@fig-inference-t1]---[-@fig-inference-t4]), it appears that adding pressure/saturation data from the wells also reduces the uncertainty in areas relatively far away from the well-bore. Finally, the well-only inference seems to include increased CO~2~ saturation at the seal, a feature that is only partly present in the ground truth. The presence of this feature can attribute to the average behavior of the permeability, an observation shared by the unconditioned results. It is encouraging to see that this feature at the seal has large uncertainty judged by the large $\bar{\widehat{\boldsymbol{\sigma}}}_{24}$ in this region. Except for an area near the well-bore, CO~2~-plume saturations are not visible on the conditional mean estimates conditioned on the seismic. However, there is a small imprint on the plots for the uncertainty, which means that the posterior samples do include sporadic examples of detectable CO~2~ saturation in this narrow area near the seal.

While these inference results by the Digital Shadow are encouraging, they only provide snapshots at the four timesteps time-lapse data is collected. They also do not provide quantitative information on the performance of the proposed Digital Shadow.


::: {#fig-ground-truth-multimodal}

![](./figs/Ground-truth-multimodal.png)

Ground-truth $\mathbf{x}^\ast_k[S_{\mathrm{CO_2}}]$ (first column) and $\mathbf{x}^\ast_k[\delta p]$ (second column) for $k=1\cdots 4$, including locations for the wells. The corresponding seismic observations, $\mathbf{y}^{\mathrm{obs}}_k$, are included in the third column. Notice that the pressure difference field precedes the CO~2~ saturation.
:::

::: {#fig-inference-t1}
![](./figs/TL1_meanerrvar_s.png)

At k=1, (from left to right column) plots for the posterior mean, the error, and standard deviations of the inferred CO~2~ saturations. Unconditioned CO~2~ saturation (first row); CO~2~ saturation conditioned on pressure and saturation wells (second row); on seismic data (third row); and on both seismic and wells (fourth row).
:::

::: {#fig-inference-t2}
![](./figs/TL2_meanerrvar_s.png)

At k=2, (from left to right column) plots for the posterior mean, the error, and standard deviations of the inferred CO~2~ saturations. Unconditioned CO~2~ saturation (first row); CO~2~ saturation conditioned on pressure and saturation wells (second row); on seismic data (third row); and on both seismic and wells (fourth row).
:::

::: {#fig-inference-t3}
![](./figs/TL3_meanerrvar_s.png)

At k=3, (from left to right column) plots for the posterior mean, the error, and standard deviations of the inferred CO~2~ saturations. Unconditioned CO~2~ saturation (first row); CO~2~ saturation conditioned on pressure and saturation wells (second row); on seismic data (third row); and on both seismic and wells (fourth row).
:::

::: {#fig-inference-t4}
![](./figs/TL4_meanerrvar_s.png)

At k=4, (from left to right column) plots for the posterior mean, the error, and standard deviations of the inferred CO~2~ saturations. Unconditioned CO~2~ saturation (first row); CO~2~ saturation conditioned on pressure and saturation wells (second row); on seismic data (third row); and on both seismic and wells (fourth row).
:::

## Quantitative results

Since the inference results produced by the Digital Shadow will be used during decision making by a Digital Twin, a better understanding needs to be gained on the performance of the Digital Shadow by means of a quantitative analysis of the different inference strategies. To this end, we include in @fig-DT-Statistics calculations of the SSIM-ERR defined as $1-\mathrm{SSIM}(\cdot)$, the mean absolute error (MAE), and the mean standard deviation ($\bar{\sigma}_k$), for all $24$ timesteps spanning a total of $1920$ days. To compensate for biases induced by the growth of the CO~2~ plumes themselves, the following relative measures will be used for the errors and standard deviations:

$$
\Delta\mathrm{RMSE}_k(\widehat{\mathbf{x}}_k) = \frac{\mathrm{RMS}(\mathbf{e}_k)}{\mathrm{RMS}(\mathbf{x}^\ast_k)}\quad \text{and}\quad \Delta\mathrm{STD}_k=\frac{\mathrm{mean}(\bar{\mathbf{\widehat{\boldsymbol{\sigma}}}}_k)}{\mathrm{RMS}(\widehat{\mathbf{\bar{x}}}_k)}
$${#eq-relative-metrics}

where $\mathrm{RMS}$ stands for Root-Mean Square. The normalization of the $\mathrm{RMS}$ for the error and the mean standard deviation by the $\mathrm{RMS}$ of the ground truth and estimated conditional mean is meant to correct for increased variability exhibited by growing CO~2~ plumes.

Each subplot in @fig-DT-Statistics contains estimates for these quantities based on inferences conditioned on seismic only (dashed blue lines), on wells only (dotted orange lines), and seismic + wells (dashed-dotted green lines). To exemplify the importance of conditioning fluid-flow simulations on time-lapse observations, we also plot estimates obtained without conditioning (solid black lines) obtained by running @eq-dynamics recursively for $k=1\cdots 4$. From these three plots, the following observations can be made.

First, the curves for the SSIM-Error, $\Delta\mathrm{RMSE}_k$, and $\Delta\mathrm{STD}_k$ overall show improved quality in all three metrics when constrained by time-lapse observations (all other than black curves). Moreover, significant enhancements as a result of the $\textsc{Analysis}$ step can be observed thanks to time-lapse data-informed corrections on the predictions of the $\textsc{Forecast}$ step. These corrections start by being relatively small when conditioning on the well measurements alone (cf. the orange star (corrected) and pentagon symbols) but become more significant when conditioning on the seismic (cf. the blue cross and square symbols) and even more pronounced when conditioning on well + seismic (cf. green diamond and round symbols). As time evolves, the magnitude of the data-informed corrections increases for all metrics with the largest improvement in $\mathrm{SSIM-Err}$ at the last timestep, which can be explained because the bulk of the CO~2~ plume will have reached the monitor well at that time. Since the state dynamics and observation model are both nonlinear, this significant improvement can not be explained by the corrections from the well data alone. 

Second, for each timestep in between the time-lapse data collection, there is, as expected, a general downward trend in plume recovery quality from the time the corrections are made until the next timestep during which the next time-lapse data is collected and the subsequent corrections are made. This deterioration in recovery quality is observed for the $\mathrm{SSIM-Err}$ for all modalities and is to leading order the result of differences between the randomly drawn permeability field and the ground truth that lead during the intermediate timesteps to increasing deviations between the actual CO~2~ plume and the predicted ensemble produced by the $\textsc{Forecast}$ step. The normalization in the definition of $\Delta\mathrm{RMSE}_k$ in @eq-relative-metrics compensates partially for this effect, yielding more prominent improvements, especially between conditioning on seismic data alone or on wells + seismic. In addition, relative errors after correction remain reasonably flat and even decrease at the final timestep. This suggests that relative errors remain more or less constant thanks to the conditioning on multimodal time-lapse observations.

Third, the relative quantity, $\Delta\mathrm{STD}_k$, aims to answer the question of whether including time-lapse data, and multimodal time-lapse data in particular, reduces the uncertainty. Thanks to the normalization by the $\mathrm{RMS}$ of the (conditional) mean, the uncertainty plotted in @fig-DT-Statistics (c) remains relatively constant with slightly larger corrections produced by conditioning the CO~2~ plumes on multimodal time-lapse data. While this behavior is consistent with the behavior displayed by the relative errors in @fig-DT-Statistics (b), the uncertainty measured in terms of the variability amongst the intermediate plume predictions seems to dip further during the immediate timesteps after the corrections are made. This behavior is observed for the unconditioned and conditioned by well-only modalities and can be explained by the fact that after each correction new samples for the permeability distribution field are drawn at which point the flow will first explore regions of high permeability before growing. This explains the observed reduction in variability for the unconditioned and conditioned on well-data-only modalities that produce CO~2~ plumes that are more widely spread and therefore can explore more regions to fill in. After these regions are filled, the plumes grow again increasing the variability and therefore the uncertainty.

In summary, the overall trends observed in @fig-DT-Statistics, are consistent with a highly nonlinear data-assimilation problem where the unknown permeability field that controls the flow is modeled as a random nuisance parameter. Even though there is a slight deterioration in the quality of the inferred CO~2~ plumes, the observed relative errors and uncertainty remain relatively stable over time when the predicted simulations are conditioned on time-lapse data that includes seismic surveys. While including time-lapse well data alone does not lead to significant improvements in CO~2~ plume reconstruction quality and uncertainty, combining this data modality with time-lapse seismic data leads to significantly improved results.

::: {#fig-DT-Statistics layout-ncol="3"}
![a.](./figs/DT-SSIM.png)

![b.](./figs/DT-RMS.png)

![c.](./figs/DT-STD.png)

Quantitative performance metric of CO~2~ plume recovery with the proposed Digital Shadow. Before and after corrections by the $\textsc{Analysis}$ step are denoted by blue square/cross symbols for seismic/corrected seismic; the orange pentagon/star symbols for wells/corrected wells, and the green circles/diamodns for seismic + wells/corrected seismic + wells.  *(a)* $\mathrm{SIMMD-Err}$ as a function of time since injection; *(b)* $\Delta\mathrm{RMSE}_k$, and *(c)* $\Delta\mathrm{STD}_k$.

:::

# Discussion {#sec-discussion}

The main goal of this work is the design, implementation, and *in sillico* validation of a Digital Shadow prototype for Geological Carbon Storage. IBM's definition of a Digital Twin reads:

> A digital twin is a virtual representation of an object or system designed to *reflect a physical object accurately*. It spans the object's *lifecycle*, is *updated from real-time data* and uses *simulation, machine learning* and reasoning to help make decisions.

When considering this definition, it is clear that the presented neural data-assimilation framework does indeed capture several key aspects (denoted by italics) of a Digital Twin. Because our formulation does not yet support control, reasoning, and decision-making, we follow @asch2022toolbox and use the term Digital Shadow for our prototype monitoring framework instead. To help future decision-making and optimize underground CO~2~ storage operations, we leverage recent breakthroughs in neural posterior-density estimation and generative AI to make our Digital Shadow uncertainty aware, which will help the Digital Twin with future tasks including underground monitoring and control, vital for cost-effective risk mitigation and management.

By taking a neural approach, where the simulated state-observation predictions of the well-established ensemble Kalman filter (EnKF, @asch2022toolbox) are considered as training samples, we arrive at a formulation where the important $\textsc{Analysis}$ step no longer relies on empirical estimates for the first- and second-order moments (mean and covariance). Instead, our approach hinges on generative conditional neural networks' ability to characterize non-Gaussian posterior distributions with high degrees of freedom. These distributions result from nonlinear models for the state dynamics (multi-phase flow) and observations (e.g. seismic images), which contain probabilistic components that give rise to complex uncertainties. Compared to other ensemble-based data-assimilation methods [@frei2013bridging], which are typically contingent on Monte-Carlo approximations of the classical Kalman recursions together with limiting Gaussian assumptions, our nonlinear extension of the EnKF enjoys certain important advantages. Backed by the recent work of @spantini2022coupling, there are strong indications that our approach does not suffer from an unresolvable lack of accuracy of EnKF whose accuracy can not be improved by increasing the ensemble size. Similar to @spantini2022coupling, we identify corrections during the stochastic prior-to-posterior map as a Kalman-like $\textsc{Analysis}$ step, implemented in the Gaussianized latent space of a conditional Invertible Neural Network (cINN, @cnf). Compared to existing ensemble-based methods, the use of conditional neural networks has the advantage that their accuracy, and hence the latent-space correction, increases when adding computational training effort, e.g. by growing the ensemble size or by training with more epochs. When the network is adequately expressive, these accuracy improvements are likely when the set of simulated training pairs, which correspond to the ensemble, are sufficiently large, diverse, and close in distribution to the actual state and observations. 

Throughout each lifecycle, the Digital Shadow trains its networks with Simulation-Based Inference (SBI, @cranmer2020frontier) by using nonlinear simulations of the time-advanced state (calculated with @eq-dynamics) and associated, possibly multimodal, observations (calculated with @eq-obs). Because our approach is simulation-based, uncertainties associated with the previous state, the random permeability field, and observational noise are marginalized implicitly [@cranmer2020frontier]. As a result, conditional neural networks act as amortized surrogates for the posterior of the state given (multimodal) observations. Thanks to the marginalization, contributions due to the stochasticity of the permeability distribution and other noise terms are factored into the resulting posterior, which is no longer dependent on this nuisance parameter and on specifics of the noise. Moreover, theoretical connections exist between the EnKF's $\textsc{Analysis}$ step and the proposed stochastic latent space prior-to-posterior mapping. Similar to @spantini2022coupling, this stochastic mapping reads in the latent space of the CNF as follows:

$$
\mathbf{x}_k^{(m)}\mapsto \widehat{\mathbf{x}}_k^{(m)}= f^{-1}_{\widehat{\boldsymbol{\phi}}_k}(\mathbf{z}_k^{(m)};\mathbf{y}_k^{\mathrm{obs}})\quad\text{with}\quad \mathbf{z}_k^{(m)}= f_{\widehat{\boldsymbol{\phi}}_k}(\mathbf{x}_k^{(m)};\mathbf{y}_k^{(m)}). 
$$ {#eq-analysis-equivalence}

Unlike the EnKF, this correction is nonlinear and corresponds to a transport map where samples from the joint distribution for the state and observations are mapped to the latent space, followed by conditioning on the observed time-lapse data during inversion from the latent space back to the physical state. As reported by @ramgraber2023ensembleb, conditioning on the observed data in the latent space can have certain advantages thanks to the Gaussinization of the latent space, a topic which we will further explore in future work. Finally, the use of physics-based summary statistics in our approach to our knowledge also differs from current data-assimilation approaches where dimensionality reduction approaches are more common for the state and not for the observations.

Even though the current implementation of our Digital Shadow is specialized towards multimodal monitoring of Geological Carbon Storage, the presented neural data-assimilation approach is general and ready to be applied to other fields as long as implementations for the offline simulation-based $\textsc{Forecast}$ step, conditional neural networks, and network training procedure are available to approximate the prior-to-posterior mapping. However, specific choices for the implementation of our Digital Shadow prototype were made contingent on scalability, computational cost, memory footprint, and accuracy of the multi-phase flow and wave simulations. For this reason, our implementation is based on the low-memory footprint Julia implementation for CNFs, `InvertibleNetworks.jl` [@Orozco2024], in combination with the industry-strength parallel multi-phase fluid-flow simulator, [JutulDarcy.jl](https://github.com/sintefmath/JutulDarcy.jl) [@møyner2023], and parallel wave-based Julia Devito Inversion framework inversion package [JUDI.jl](https://github.com/slimgroup/JUDI.jl),[@witte2018alf; @JUDI]. While our framework allows for other choices, our prototype for the Digital Shadow benefits from the intrinsic parallel capabilities of the fluid-flow and wave simulators and from the fact that the ensemble calculations are embarrasingly parallel. In that respect, the presented scalable approach differs from recent score-based work by [@rozet2023score; @bao2024nonlinear; @huang2024diffusionbasedsubsurfacemultiphysicsmonitoring], which either rely on many, and therefore fast evaluations of the dynamics and observation operators, or on fast access to the complete history of the CO~2~ plume's dynamics, neither of which are conducive to realistic scale data-assimilation for GCS projects. 

Because of the relatively long lead time between the collection of time-lapse seismic surveys, we so far did not opt for the use of our parallel implementation of Fourier Neural Operators (FNO, @grady2023model), which in principle can serve as a suitable substitute to the expensive multi-phase flow simulations [@louboutin2023lmi; @yin2023solving]. We also reserve inclusion of geomechanical, chemical, and (poro-) elastic effects, to future research even though these effects can in principle be included in the Digital Shadow's simulations. These extensions would allow for the use of our Digital Shadow in more complex situations, including GCS on land, where geomechanics could play a more important role, or within reservoirs that undergo geochemical changes.

Compared to contemporary approaches to CO~2~ monitoring [@huang2020towards; @grana2020prediction], the proposed Digital Shadow's conditional neural networks are capable of capturing the nonlinear state dynamics, given by the multi-phase flow, and nonlinear seismic observation model, involving nonlinear solutions of the wave equation. Our current implementation achieves this at a scale where extended Kalman filters struggle to store explicit matrix representations for the Kalman filter that includes explicit inverses of large matrices. Like recent work by @bruer2024seismicmonitoringco2plume, the proposed approach reaps information from an ensemble of states, which is an approach that scales well. By treating the permeability field as a stochastic nuisance parameter, which is drawn independently at each time-lapse timestep for each ensemble member, our Digital Shadow distinguishes itself from this latest work on an EnKF for GCS where the randomly drawn permeability field remains fixed within each member of the ensemble. Aside from incorporating stochasticity induced by the unknown permeable field, which can be interpreted as a modeling error connecting this work to @si2024latentensflatentensemblescore, our approach is also based on a nonlinear seismic observation model. 

Finally, the presented approach made assumptions on the distribution of the permeability field as well as on the underlying rock and wave physics. In future work, we will relax these assumptions by *(i)* extending recent work by @yin2024time to carry out permeability inference from time-lapse seismic and/or well data, which would allow for narrowing of the baseline probability distribution for the permeability field, as time-lapse data becomes available; *(ii)* marginalizing over uncertainties in the rock-physics model, e.g. by augmenting the Forecast ensemble with samples generated from models other than the patchy saturation model [@avseth2010quantitative]; *(iii)* calibrating rock-physical relations at (monitoring) wells [@feng2022extremelyweaksupervisioninversion]; *(iv)* robustifying the inference with respect to distribution shifts either through physics-based latent-space corrections, as proposed by @yin2024wiser and @orozco2024aspireiterativeamortizedposterior, or by adapting score-based methods to solve the weak formulation proposed by @siahkoohi2020weak and @orozco2024aspireiterativeamortizedposterior; *(iv)* no longer insisting on costly seismic survey replication to achieve repeatability of time-lapse seismic surveys [@oghenekohwo2017low; @wason2017low; @lopez2023spectral;@zhang2023optimized] *(v)* increasing detectibility of time-lapse changes in the density by using SH-polarized shear waves [@clochard2018quadri].

# Conclusions {#sec-conclusion}

The inherent high degrees of freedom and complications arising from the intricate interplay between the nonlinear multi-phase flow, nonlinear rock physics, and possibly nonlinear multimodal time-lapse observations, challenge the development of monitoring systems for subsurface storage of supercritical CO~2~. In this work, we meet these challenges by introducing an uncertainty-aware Digital Shadow, which combines state-of-the-art high-fidelity simulation capabilities with data-assimilation and machine-learning techniques to arrive at a scalable formulation that through the deployment of a stochastic nuisance parameter incorporates nonlinearities and inherent uncertainties in the permeability field. Thanks to the Digital Shadow's use of simulation-based inference, its neural networks marginalize over this stochastic permeability field, factoring lack of knowledge on this important reservoir property into Digital Shadows' uncertainty quantification. As a result, dependence on the permeability field is eliminated and its uncertainty is absorbed into the recursively inferred posterior distributions for the CO~2~ plume's state (pressure/saturation), conditioned at each timestep on multimodal time-lapse measurements. The training itself takes place recurrently on $\textsc{Forecast}$ ensembles, consisting of predictions for the state and simulations of induced time-lapse data collected at the surface and/or in wells.

To validate the proposed Digital Shadow, a realistic simulation-based case study was conducted, starting from a probabilistic model for the permeability field derived from a baseline surface seismic survey. Next, the Digital Shadow's network was trained on $\textsc{Forecast}$ ensembles, consisting of simulations for the time-advanced CO~2~ plume's state and corresponding multimodal time-lapse data. After training was completed, the Digital Shadow's networks were applied to unseen simulated ground-truth time-lapse data. From this computational case study, the following observations can be made. First, CO~2~plume estimates based on averages over the $\textsc{Forecast}$ ensembles behave poorly compared to estimates corrected by observed time-lapse observations. These improvements are achieved irrespective of the permeability field. Second, the crucial imprint of an unknown high-permeability streak was missed when the CO~2~ plume's state was conditioned on pressure/saturation measurements at wells alone. Plume inferences conditioned on seismic or seismic and well measurements succeeded in detecting this event that ultimately leads to Containment failure. Third, the inferred pointwise standard deviation amongst the samples of the posterior, which serves as a measure of their uncertainty, correlates well with errors with respect to the ground-truth state. This property persists over the different timesteps and is an indication that the uncertainty quantification is reasonably well calibrated. Fourth, visual as well as quantitative measures of the reconstruction quality show across the board the best results when the CO~2~ plume's state is conditioned on both time-lapse seismic and well measurements. While conditioning the CO~2~ plume on wells only leads to marginal improvements compared to uncorrected forecasts, combination with seismic leads to disproportional improvements even in areas away from the wells. We attribute these global improvements in the inferred state of the CO~2~ plume to the strong nonlinearity of multi-phase flow. All in all, we argue from this case study that the proposed Digital Shadow tracks the unseen ground-truth physical state accurately for the duration of a realistic CO~2~injection project. To our knowledge, the proposed methodology constitutes the first proof-of-concept of an uncertainty-aware Digital Shadow that principally encapsulates uncertainty due to unknown reservoir properties and noise present in the observed data. This Digital Shadow will form the basis for a Digital Twin designed to mitigate risks and optimize operations of underground storage projects. 


# Data and Software Availability {.unnumbered}

The data underlying the findings of this study and the software developed for it will be made publicly available upon acceptance of the manuscript for publication. Also, full movies of the snapshots shown in the manuscripts will be included during final submission in the online supplementary material. 

# Acknowledgement {.unnumbered}

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

***

# References {.unnumbered}

::: {#refs}
:::

\setcounter{section}{0}
\renewcommand{\thesection}{\Alph{section}}

\setcounter{table}{0}
\renewcommand{\thetable}{A\arabic{table}}

\setcounter{figure}{0}
\renewcommand{\thefigure}{A\arabic{figure}}

# Appendices {.unnumbered}

## Appendix A {.appendix .unnumbered}

## Network architecture cINNs {.appendix .unnumbered #sec-network}

For CNFs, we followed conditional GLOW network architecture [@kingma2018glow] implemented in the open-source package [InvertibleNetworks.jl](https://github.com/slimgroup/InvertibleNetworks.jl) [@InvertibleNetworks]. The architecture of Conditional GLOW typically consists of invertible neural network layers, coupled with affine coupling layers, which enable efficient computation of the determinant of the Jacobian matrix. These layers are organized in a hierarchical manner, which offers better expressiveness when compared to the conventional invertible architectures [@dinh2016density]. This expressivity comes from applying a series of conventional invertible layers [@dinh2016density] to the input in a hierarchical manner which results in an invertible architecture exhibiting the potential to depict intricate bijective transformations [@siahkoohi2023reliable; @orozco2023invertiblenetworks] and hence the ability to learn complex distributions.

### Training Setting {.appendix .unnumbered}

| Hyperparameters           | Values                              |
|:--------------------------|:------------------------------------|
| Batch Size                | $8$                                 |
| Optimizer                 | $\mathrm{Adam}$ [@Kingma2014AdamAM] |
| Learning rate (LR)        | $5^{-4}$                            |
| No. of training epochs    | $120$                               |
| Fixed Noise Magnitude     | $0.005$                             |
| No. of training samples   | $120$                               |
| No. of validation samples | $8$                                 |
| No. of testing samples    | $1$                                 |

: Training Hyperparameters {#tbl-hyperparameters}

## Appendix B {.appendix .unnumbered #sec-appendix-B}

### Useful Definitions {.appendix .unnumbered}

After the completion of training, we use the following procedure to calculate the posterior mean: 

$$
\mathbf{x}_{\text{pm}} = \mathbb{E}_{\mathbf{x} \sim p(\mathbf{x}|\mathbf{y})} [ \mathbf{x} ]
\approx 
\frac{1}{M}  \sum^{M}_{n=1}\mathbf{x}_{\text{gen}}^{i} \text{ where } \mathbf{x}_{\text{gen}}^{i} = f_{\hat{\boldsymbol{\phi}}}^{-1}(\mathbf{z_i};\mathbf{y}) \text{ with } \mathbf{z_i} \sim \mathcal{N}(0,I),
$$ {#eq-posteriormean}

where $\hat{\boldsymbol{\phi}}$ is the minimizer of Equation @eq-loss-CNF.

<!-- Metric definitions used in @fig-inference-t1 through @fig-inference-t4 and text -->

$\textbf{SSIM}$ - Structural Similarity Index Metric quantifies the similarity between two images and is commonly used to assess how closely a generated image resembles a ground truth or reference image. It considers image quality aspects such as luminance, contrast, and structure. For the mathematical formulation of SSIM, please refer to the study by [@1284395].

$\textbf{RMSE}$ - Root mean squared error is used to represent the measure of the difference between ground truth CO$_2$ saturation image and the posterior mean of the samples generated by the trained network.

## Appendix C {.appendix .unnumbered #sec-appendix-C}

$\textbf{Pressure perturbation}$ - Total reservoir pressure, excluding hydrostatic pressure. 

Similar to CO~2~ saturation states, estimates for the CO~2~ pressure perturbation, $\bar{\widehat{\mathbf{x}}}_{1:6:24}[\delta p]$, $\mathbf{e}_{1:6:24}=|\mathbf{x}^\ast_{1:6:24}[\delta p]-\mathbf{\bar{\widehat{x}}}_{1:6:24}[\delta p]|$, and $\bar{\widehat{\boldsymbol{\sigma}}}_{1:6:24}$ are, for each of the four scenarios (unconditioned, conditioned on well only, seismic only, well + seismic) included in the rows of figures [-@fig-inference-t1P]---[-@fig-inference-t4P] for $k=1\cdots 4$. The columns of these figures contain estimates for the conditional mean, errors, and uncertainty in terms of the standard deviation. From these plots, the following qualitative observations can be made. 


Unlike the CO~2~ saturation states, the pressure perturbations estimates are relatively smooth similar to the ground truth (cf. @fig-ground-truth-multimodal, second column), so the conditional mean estimates (first column of figures [-@fig-inference-t1P]---[-@fig-inference-t4P]) do not appear to differ substantially amongst the four different scenarios during the four time-lapse time steps. As observed for infered saturations, variations between the posterior samples conditioned on the wells or wells + sesimic (third columns) show reduced uncertainty near the wells. Errors (second colums) with respect to the ground truth do not change much between the different modalities.

::: {#fig-inference-t1P}
![](./figs/TL1_meanerrvar_p.png)

At k=1, (from left to right column) plots for the posterior mean, the error, and standard deviations of the inferred CO~2~ pressure perturbation. Unconditioned pressure perturbation (first row); pressure perturbation conditioned on pressure and saturation wells (second row); on seismic data (third row); and on both seismic and wells (fourth row).
:::


::: {#fig-inference-t2P}
![](./figs/TL2_meanerrvar_p.png)

At k=2, (from left to right column) plots for the posterior mean, the error, and standard deviations of the inferred CO~2~ pressure perturbation. Unconditioned pressure perturbation (first row); pressure perturbation conditioned on pressure and saturation wells (second row); on seismic data (third row); and on both seismic and wells (fourth row).
:::


::: {#fig-inference-t3P}
![](./figs/TL3_meanerrvar_p.png)

At k=3, (from left to right column) plots for the posterior mean, the error, and standard deviations of the inferred CO~2~ pressure perturbation. Unconditioned pressure perturbation (first row); pressure perturbation conditioned on pressure and saturation wells (second row); on seismic data (third row); and on both seismic and wells (fourth row).
:::


::: {#fig-inference-t4P}
![](./figs/TL4_meanerrvar_p.png)

At k=4, (from left to right column) plots for the posterior mean, the error, and standard deviations of the inferred CO~2~ pressure perturbation. Unconditioned pressure perturbation (first row); pressure perturbation conditioned on pressure and saturation wells (second row); on seismic data (third row); and on both seismic and wells (fourth row).
:::



