---
title: "Time-lapse full-waveform permeability inversion: a feasibility study"
author:
  - name:
      Ziyi Yin^1^, Mathias Louboutin^1,3^, Olav Møyner^2^, Felix J. Herrmann^1^ \
      
      ^1^ Georgia Institute of Technology, ^2^ SINTEF Digital, ^3^ Now at Devito Codes Ltd
bibliography: paper.bib
lightbox: True
filters:
  - SLIM
execute:
  freeze: auto
format:
  html:
    number-sections: true
toc: true
abstract: Time-lapse seismic monitoring necessitates integrated workflows that combine seismic and reservoir modeling to enhance reservoir property estimation. We present a feasibility study of an end-to-end inversion framework that directly inverts for permeability from multiple prestack time-lapse seismic datasets. To assess the method's robustness, we design experiments focusing on its sensitivity to initial models and potential errors in modeling. Our study leverages the publicly available Compass model to simulate CO~2~ storage in saline aquifers. This model is derived from well and seismic data from the North Sea in an area that is currently considered for geological carbon storage.
---

## Introduction

Despite significant advancements in reservoir monitoring over recent decades, time-lapse seismic technology continues to face challenges related to cost and efficiency [@lumley2001time; @chadwick2009; @chadwick2010quantitative; @FURRE20173916]. Employing 4D seismic workflows, including time-lapse full-waveform inversion [TL-FWI, @lumley20104d; @hicks2016], has become a common practice for estimating changes in the Earth's elastic properties, facilitating the quantitative interpretation of these changes as indicators of reservoir attributes like fluid content and pressure [@bosch2010; @wei2017]. Recent methodologies aim to leverage time-lapse seismic data for the joint estimation of both elastic and reservoir properties, with a focus on parameters such as saturation and porosity [@bosch2007; @hu2022]. However, the integration of seismic imaging workflows with reservoir simulation tools remains limited, constraining the direct application of time-lapse seismic data for permeability estimation directly from multiple time-lapse seismic surveys. A few exceptions exist. For example, @eikrem2016 uses ensemble Kalman filtering to refine permeability and porosity estimates. @vasco2004 and @vasco2008 have explored using time-lapse seismic data for linearized inversion to update permeability. Despite these initial attemps, a more systematic and integrated approach for reservoir characterization and monitoring deserves further investigation.

This paper introduces a novel 4D processing framework for estimating permeability directly from prestack time-lapse seismic data, offering a streamlined, geophysics based inversion process. Unlike traditional methods, this framework, tested on various synthetic case studies [@li2020coupled; @yin2022SEGlci; @louboutin2023; @yin2023solving], updates permeability models by exclusively matching against multiple observed time-lapse seismic surveys. Despite the potential for rapid model updates, initial results have not yet demonstrated significant alterations in fluid saturation predictions, and the resulting permeability models often lack the heterogeneity necessary for detailed analysis. To address these limitations and to assess the framework's real-world applicability for 4D monitoring, we undertake a feasibility study using a 2D slice of the Compass model shown in @fig-perms [@jones2012building]. The geological structures of this model were derived from well logs and imaged seismic from the South-West North Sea area --- a region under consideration for CO~2~ storage [@kolster2018impact; @meneguolo2024subsurface]. This region comprises a storage unit composed of Bunter sandstone $300-500 \mathrm{m}$ thick, depicted in orange in @fig-acquisition and characterized by high permeability values shown in @fig-true-perm. On top of the storage unit there is a primary seal about $50\,\mathrm{m}$ thick, made of the low-permeable Rot Halite Member, and a secondary seal, over $300\,\mathrm{m}$, made of low-permeability mudstone in the Haisborough Group, illustrated in blue and green in @fig-acquisition.

::: {#fig-perms layout="[[99],[99],[99]]"}
![](./figs/acquisition.png){#fig-acquisition width="60%"}

![](./figs/true-perm.png){#fig-true-perm width="60%"}

![](./figs/hist-perm.png){#fig-hist-perm width="60%"}

Experimental configuration. (a) Setup of seismic acquisition and well control. Dark blue $\star$ denotes the CO~2~ injection well. While $\times$ and black $\cdot$ represent source and receiver locations, respectively. The gray curve delineates the shape of the CO~2~ plume at the 25th year. (b) Unseen ground truth spatial distribution of horizontal permeability. (c) Histogram of the common logarithm of the permeability model with Kernel Density Estimation (KDE).
:::

This study evaluates the coupled fluid-flow, rock physics, seismic inversion framework's sensitivity to different starting models, forward modeling errors, and crosstalk during multiparameter inversion, omitting regularization techniques to focus on the impact of time-lapse seismic data on permeability updates. We explore the framework's ability to recover relatively fine-scale permeability structures, predict CO~2~ dynamics within the seismic monitoring period, and forecast CO~2~ dynamics in near future without any seismic observation. Recognizing the limitations of our simplifying assumptions, we conclude with suggestions for future research to advance this promising approach.

## Permeability inversion framework

Our feasibility study examines the time-lapse seismic monitoring of geological carbon storage (GCS), focusing on the integration of three fundamental physics disciplines: fluid-flow physics, rock physics, and wave physics, as illustrated in @fig-multiphysics. The dynamics of the CO~2~ plume during injection are modeled using multiphase flow equations [@pruess2011numerical], processed through a reservoir simulator [@krogstad2015mrst; @settgast2018geosx; @rasmussen2021open; @stacey2017validation]. While these simulations require detailed inputs, including well operation parameters and the spatial distribution of porosity and permeability, in this exposition we focus on the permeability, $\mathbf{K}$, particularly, as the parameter of interest. The output from the reservoir simulator, $\mathcal{S}$, primarily the time-varying CO~2~ saturation snapshots, compiled in $\mathbf{c}$, serves as the input to the rock physics model, $\mathcal{R}$. Based on the porosity and the brine-filled baseline velocity before CO~2~ injection, this model translates each CO~2~ saturation snapshot into altered seismic velocity models, compiled in the vector, $\mathbf{v}$, using the patchy saturation model proposed by @avseth2010quantitative. Lastly, based on the velocity models for each snapshot, the wave modeling operator [@tarantola1984inversion], $\mathcal{F}$, is used to generate the time-lapse seismic dataset, $\mathbf{d}$, which collects the seismic data from each vintage.

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

Multiphysics forward model. The reservoir simulator, $\mathcal{S}$, produces $n_k$ snapshots of time-varying CO~2~ saturation, compiled in $\mathbf{c}$, from the permeability model, $\mathbf{K}$. The rock physics model, $\mathcal{R}$, based on the porosity and the brine-filled velocity model, converts each CO~2~ saturation snapshot, $\mathbf{c}_k$, to the altered velocity model, $\mathbf{v}_k$ . The shaded area highlights the CO~2~-induced changes in velocity. Finally, the wave modeling, $\mathcal{F}$, generates a time-lapse seismic dataset, $\mathbf{d}_k$, for each velocity model, $\mathbf{v}_k$. These datasets are collected in the vector, $\mathbf{d}$.
:::

In practice, the prestack time-lapse seismic dataset, $\mathbf{d}$, is observed from the field, with the objective to estimate the past, current, and future dynamics of the CO~2~ plume in GCS projects. Our methodology diverges from traditional workflows that typically proceed from seismic inversion to quantitative interpretation and subsequent reservoir parameter updates based on the derived wave properties. Instead, we propose an integrated, end-to-end approach that directly inverts the time-lapse surveys collected in $\mathbf{d}$ for permeability, $\mathbf{K}$, by reducing the time-lapse seismic data misfit objective through an automatic optimization procedure. The core of our method is the composition of three physics-based modeling operators, formulated to minimize the following objective function:

$$\underset{\mathbf{K}}{\operatorname{minimize}} \quad \|\mathcal{F}\circ\mathcal{R}\circ\mathcal{S}(\mathbf{K})-\mathbf{d}\|_2^2.$${#eq-inv}

This optimization problem is reached with the assumption that the permeability model, $\mathbf{K}$, is the only unknown parameter during inversion. Especially, the porosity and the brine-filled velocity model before CO~2~ injection (as inputs to the patchy saturation model, $\mathcal{R}$) are assumed known and fixed during the inversion. This objective is minimized via an iterative procedure that includes:

-   Generating synthetic time-lapse seismic data using an initial guess for the permeability model;
-   Calculating the gradient of the permeability by backpropagating the residuals of the time-lapse seismic datasets;
-   Updating the permeability model to reduce the misfit between the synthetic and observed time-lapse seismic datasets.

The advantage of this end-to-end inversion framework lies in its ability to break down silos through multiphysics integration. Specifically, it eliminates the need for intermediate processing steps to update the saturation and velocity models. As we demonstrate in the subsequent feasibility study, the inverted permeability can produce accurate fluid saturation and velocity models, even for the near future without any additional seismic observations.

## Feasibility study on the Compass model

We evaluate the performance of this inversion framework through a synthetic case study on the Compass model [@jones2012building]. This model has a grid spacing of $6\mathrm{m}$ in both the horizontal and vertical directions. Compared to conventional reservoir models that often have nearly homogeneous layers and a coarse discretization in the horizontal direction (e.g., $100\,\mathrm{m}$), the high-resolution and spatially heterogeneous Compass model can help reveal the potential of inverting fine-scale geological structures in the permeability model.

Using five vintages of prestack time-lapse seismic surveys, we aim to invert for the spatial distribution of permeability. To this end, we utilize the aforementioned 2D slice of the velocity model, included in @fig-acquisition, where the orange region signifies the storage unit. Since the Compass model only includes velocities and densities and no permeability or porosity values, we build a fully heterogeneous ground truth permeability model, displayed in @fig-true-perm, by assuming the elementwise relationship in @eq-vtoK, between the entries of the brine-filled velocity model, $\mathbf{c}_p$, in $\mathrm{km/s}$, and the horizontal permeability, $\mathbf{K}$, in millidarcies (md)[^1].
$$
\mathbf{K} = 
\begin{cases} 
3000 \exp(\mathbf{c}_p-4) \quad \text{if}\quad \mathbf{c}_p\geq 4 \\
0.01\exp(25.22(\mathbf{c}_p-3.5)) \quad \text{if}\quad \mathbf{c}_p\geq 3.5 \\
0.01\exp(\mathbf{c}_p-3.5) \quad \text{else}
\end{cases} 
$$ {#eq-vtoK}

Because the baseline seismic Compass model is derived from well and imaged seismic, it contains subtle changes in the seismic properties related to sub-wavelength interference. The above elementwise relationship between the seismic velocity and permeability ensures that the derived permeability model shares the same spatial heterogeneity as exhibited by the seismic baseline. By construction, it also features significant permeability contrasts within different layers within the storage unit. The low-permeability layers range from approximately $10^{-3}$ to $1$ md, while the high-permeability layers vary between $600$ to $6000$ md. The histogram of the common logarithm of the permeability model is shown in @fig-hist-perm, which demonstrates that the permeability values do not follow log-normal distribution. A CO~2~ injection well, marked with a dark blue $\star$, is placed centrally to inject supercritical CO~2~ for 25 years at a constant rate of two million metric tons per year. We assume the porosity and the kv/kh ratio to be constant and given by 25% and 10%, respectively. The simulation of compressible and immiscible two-phase flow, where CO~2~ displaces brine in porous rocks, is performed using a fully implicit method implemented in [JutulDarcy.jl](https://github.com/sintefmath/JutulDarcy.jl) [@møyner2023; @JutulDarcyRules]. The boundary of the CO~2~ plume at the 25th year is depicted in grey in @fig-acquisition. After converting the CO~2~ saturation into seismic velocity models, $\mathbf{v}$, via the patchy saturation model, acoustic time-lapse seismic data is generated with constant density for five vintages at years 5, 10, 15, 20, and 25 using [Devito](https://github.com/devitocodes/devito) [@louboutin2018dae; @luporini2020architecture] and [JUDI.jl](https://github.com/slimgroup/JUDI.jl) [@witte2018alf; @louboutin2023a], employing a Ricker wavelet with a central frequency of 20 Hz. The well-bore source and receiver geometries are shown in @fig-acquisition.

[^1]: This relationship is only used to create the ground truth permeability model. When solving the optimization problem in @eq-inv, no assumptions are made on relationships between the permeability, porosity, and velocity. 

### Sensitivity with respect to starting models

To evaluate the efficacy of our end-to-end inversion framework, particularly its sensitivity to initial permeability models, we examine two distinct initial permeability models. In case 1, the initial permeability model, shown in @fig-init-perm-1, features homogeneous permeability values ($100\mathrm{md}$) across the entire reservoir. This model allows us to explore the extent of permeability updates achievable from a non-informative permeability model. In case 2, we apply a spatial distortion [@bloice2017] to the unseen ground truth permeability in @fig-true-perm to obtain the initial permeability model, shown in @fig-init-perm-2. The values of different permeability layers are near accurate, but the positions are misplaced.

In both cases, we employ a methodological shortcut often referred to as committing an "inversion crime", where the data generation and inversion processes share the same computational kernel. This ideal setup is used here to show what is ideally achievable by this inversion framework. To add a layer of realism, we incorporate 8 dB of incoherent band-limited Gaussian noise into the observed time-lapse datasets, which severely contaminates the seismic signal in the time-lapse difference.

To invert for the permeability model, we run 100 iterations of stochastic gradient descent (SGD), starting with @fig-init-perm-1 and @fig-init-perm-2. During each iteration, we randomly draw four sources out of the total of 32 sources to calculate the misfit and the gradient with respect to permeability model [@herrmann2013frugal]. This amounts to 12.5 datapasses through the entire time-lapse seismic dataset. We display permeability updates in logarithmic scale for both cases, in @fig-case1-update-perm and @fig-case2-update-perm, respectively. Additionally, @fig-case1-ideal-update and @fig-case2-ideal-update offer a visualization of "ideal" updates by showing the logarithmic differences between the ground truth and the initial permeability models.

The following observations can be made: First, the permeability updates are primarily confined to areas directly influenced by the CO~2~ plume's flow, as delineated by the gray curves. This outcome is expected since the time-lapse variations in wave properties are attributed to changes in fluid saturation exclusively. Consequently, without additional information, this inversion method does not alter permeability values outside the CO~2~ plume's extent where the flow of CO~2~ has not occurred. Second, the inverted permeability within the CO~2~ plume largely reflects the trend of the ground truth permeability model. In case 1, the framework successfully identifies major permeability layers --- both high and low (depicted in red and blue, respectively) at approximately $1600\,\mathrm{m}$ depth --- accurately capturing their depth and lateral distribution in alignment with the actual layers. In case 2, the inversion process introduces high-resolution details to the layers affected by the plume, aligning well with the ideal updates shown in @fig-case2-ideal-update. Despite these successes, the full magnitude of permeability contrasts is not entirely captured, pointing to the inherently ill-posed nature of permeability inversion [@zhang2014], necessitating workflows that include uncertainty quantification for future investigations [@gahlot2023inference].

::: {#fig-inv-perm layout="[[49,49],[49,49],[49,49]]"}
![](./figs/init-perm-1.png){#fig-init-perm-1}

![](./figs/init-perm-2.png){#fig-init-perm-2}

![](./figs/case1-ideal-update-perm.png){#fig-case1-ideal-update}

![](./figs/case2-ideal-update-perm.png){#fig-case2-ideal-update}

![](./figs/case1-update-perm.png){#fig-case1-update-perm}

![](./figs/case2-update-perm.png){#fig-case2-update-perm}

Permeability inversion results for cases 1 and 2. (a)(c)(e) display the initial permeability model in case 1, the logarithmic ratio of the ground truth permeability (@fig-true-perm) to the initial one, and the logarithmic ratio of the inverted permeability to the initial one. (b)(d)(f) display the same but for case 2 with a distorted initial permeability model. Gray curve indicates the boundary of the CO~2~ plume at the 25th year. "$\times$ initial" on the caption of the colorbar represents the factor by which the initial permeability is updated.
:::

### Sensitivity with respect to forward modeling errors

To extend our investigation beyond overly idealized scenarios, we examine the framework's robustness next, during scenarios that avoid "committing the inversion crime". A critical area of focus in time-lapse seismic is the error in the brine-filled baseline velocity model before CO~2~ injection. Errors in this baseline model, which feeds into the rock physics model, can produce inaccurate velocity models of the CO~2~-filled reservoir, leading to inaccuracies in the simulated time-lapse seismic datasets.

To construct an inaccurate but realistic baseline velocity model, we use a Ricker wavelet with central frequency of $20\,\mathrm{Hz}$ to generate cross-well and surface seismic data before CO~2~ injection, and employ the SGD method to run 10 datapass of FWI with a kinematically correct but smooth initial velocity model to obtain the inverted velocity model depicted in @fig-fwi-v. This imperfect brine-filled velocity model is subsequently fed into the rock physics model, $\mathcal{R}$, for permeability inversion (case 3).

We employ the initial permeability model from @fig-init-perm-1 to assess the impact of modeling errors on the inversion results. The update to the permeability, shown in logarithmic scale in @fig-non-inv-crime-perm-update, reveals some artifacts outside the CO~2~ plume area due to the modeling inaccuracies. Additionally, a high permeability zone within the plume is slightly misplaced when compared to the updates in @fig-case1-update-perm, yet the overall trend of permeability changes is correctly captured.

Following this permeability update, we proceed with reservoir simulations using the updated permeability models to assess the corrections made to the CO~2~ plume predictions. This step is crucial for validating the practical utility of the inversion framework for real-world seismic monitoring scenarios.

::: {#fig-non-inv-crime layout="[[49,49]]"}
![](./figs/fwi-v.png){#fig-fwi-v}

![](./figs/case3-update-perm.png){#fig-non-inv-crime-perm-update}

Permeability inversion results for case 3. (a) Inverted brine-filled baseline velocity used in permeability inversion. (b) The logarithmic ratio of the inverted permeability to the initial one. Gray curve indicates the boundary of the CO~2~ plume at the 25th year.
:::

### CO~2~ plume estimation and forecast

The primary objective of our end-to-end inversion framework is to accurately estimate reservoir permeability, a crucial step towards the ultimate goal of predicting CO~2~ saturation both historically and in the near future. Based on initial, inverted, and ground truth permeability models, we conduct a quality control involving CO~2~ saturation simulations, as depicted in @fig-co2. Across all simulations, we note substantial improvements in predictions of the CO~2~ plume shape, closely aligning with the boundaries of the ground truth CO~2~ plume. Notably, the initial simulations significantly misjudged the lateral spread of the CO~2~ plume. The corrections applied through the updated permeability models, however, yield accurate representations of the plume's lateral extent.

::: {#fig-co2 layout="[[22,22,22],[22,22,22],[22,22,22],[22,22,22],[22,22,22],[22,22,22]]"}
![](./figs/true-co2-5yr.png){width=90%}

![](./figs/true-co2-15yr.png){width=90%}

![](./figs/true-co2-25yr.png){width=90%}

![](./figs/case1-init-co2-5yr.png){width=90%}

![](./figs/case1-init-co2-15yr.png){width=90%}

![](./figs/case1-init-co2-25yr.png){width=90%}

![](./figs/case1-co2-5yr.png){width=90%}

![](./figs/case1-co2-15yr.png){width=90%}

![](./figs/case1-co2-25yr.png){width=90%}

![](./figs/case2-init-co2-5yr.png){width=90%}

![](./figs/case2-init-co2-15yr.png){width=90%}

![](./figs/case2-init-co2-25yr.png){width=90%}

![](./figs/case2-co2-5yr.png){width=90%}

![](./figs/case2-co2-15yr.png){width=90%}

![](./figs/case2-co2-25yr.png){width=90%}

![](./figs/case3-co2-5yr.png){width=90%}

![](./figs/case3-co2-15yr.png){width=90%}

![](./figs/case3-co2-25yr.png){width=90%}

Predicted CO~2~ saturation for 5th, 15th, and 25th years, shown in first, second, and third columns, respectively. The first row shows the (unseen) ground truth CO~2~ saturation. The second and fourth rows show the saturation predicted with initial permeability models in @fig-init-perm-1 and @fig-init-perm-2, respectively. The third, fifth, and sixth rows show the updated saturation after updating the initial permeability models by @fig-case1-update-perm, @fig-case2-update-perm, and @fig-non-inv-crime-perm-update, respectively. The boundaries of the (unseen) ground truth CO~2~ saturation are shown in white curves.
:::

Expanding our analysis to future forecasting, @fig-co2-forecast illustrates the predicted movement of the CO~2~ plume over a 40-year period, following a 25-year injection phase, without further CO~2~ injection or seismic observations. During this forecasted period, the CO~2~ plume primarily ascends due to buoyancy, while a portion (approximately 10%) remains trapped in the pore spaces, indicated in purple. This phenomenon, known as residual trapping [@rahman2016], is a critical factor in assessing the long-term storage capabilities of GCS projects. Initial forecasts tend to underestimate the extent of CO~2~ sequestration through residual trapping. In contrast, simulations driven by the updated permeability models not only provide a more accurate estimation of the permanently stored CO~2~ volume but also closely match the ground truth CO~2~ plume's boundaries, even without collecting further monitoring data.

::: {#fig-co2-forecast layout="[[49, 49],[49, 49],[49, 49],[49, 49],[49, 49],[49, 49]]"}
![](./figs/true-co2-45yr.png){width=90%}

![](./figs/true-co2-65yr.png){width=90%}

![](./figs/case1-init-co2-45yr.png){width=90%}

![](./figs/case1-init-co2-65yr.png){width=90%}

![](./figs/case1-co2-45yr.png){width=90%}

![](./figs/case1-co2-65yr.png){width=90%}

![](./figs/case2-init-co2-45yr.png){width=90%}

![](./figs/case2-init-co2-65yr.png){width=90%}

![](./figs/case2-co2-45yr.png){width=90%}

![](./figs/case2-co2-65yr.png){width=90%}

![](./figs/case3-co2-45yr.png){width=90%}

![](./figs/case3-co2-65yr.png){width=90%}

CO~2~ plume forecasts for 45th and 65th years, shown in first and second columns, respectively. The ordering of the rows remains the same as @fig-co2. Purple regions display the CO~2~ plume permanently stored via the residual trapping mechanism.
:::

### Multiparameter inversion

While case 1-3 demonstrate the performance of the inversion framework for permeability estimation, further scrutiny is in order to investigate its performance when multiple parameters are unknown and need to be jointly estimated. To this end, we design case 4 where spatial distributions of both porosity and permeability are unknown. While permeability only appears in the two-phase flow equations, porosity appears as an input not only in the two-phase flow reservoir simulator, $\mathcal{S}$, but also in the patchy saturation model, $\mathcal{R}$. This results in a more challenging inverse problem because porosity affects more than one physics-based forward models in @fig-multiphysics, and because there can be crosstalk in the gradient calculations during the inversion.

As a proof of concept, we simplify this multiparameter inversion problem by assuming that the porosity, $\boldsymbol{\phi}$, and the horizontal permeability, $\mathbf{K}$, are related by the following elementwise Kozeny-Carman relationship [@costa2006permeability]:

$$
\mathbf{K} = \mathcal{T}(\boldsymbol{\phi})=3.65 \times 10^4 \frac{\boldsymbol{\phi}^3}{(1-\boldsymbol{\phi})^2}.
$$ {#eq-phitoK}

Following this relationship, we artificially create a ground truth porosity model, shown in @fig-true-poro, according to the permeability values in @fig-true-perm, and then use the ground truth permeability and porosity models to simulate CO~2~ saturation, velocity models, and the five time-lapse seismic datasets. During inversion, we parameterize permeability by porosity according to @eq-phitoK, and minimize the following objective function to invert for the porosity:

$$
\underset{\boldsymbol{\phi}}{\operatorname{minimize}} \quad \|\mathcal{F}\left(\mathcal{R}\left(\boldsymbol{\phi}, \mathcal{S}\left(\boldsymbol{\phi}, \mathcal{T}\left(\boldsymbol{\phi}\right)\right)\right)\right)-\mathbf{d}\|_2^2.
$$ {#eq-multiparam}

The scalable and differentiable programming framework, proposed by @louboutin2023, allows for effortless and accurate gradient calculation with respect to porosity, $\boldsymbol{\phi}$, which otherwise requires labor-intensive and error-prone derivation of cross-gradient terms by hand. We initialize the reservoir with homogeneous porosity values of $12\%$, shown in @fig-init-poro. After 100 iterations of SGD, the inverted porosity is shown in @fig-inv-poro. While some layers in the inverted porosity are slightly misplaced in this preliminary study, the overall trend of porosity is adequately estimated in the center of the model.

::: {#fig-poro layout="[[33, 33, 33]]"}
![](./figs/init_poro.png){#fig-init-poro}

![](./figs/inv_poro.png){#fig-inv-poro}

![](./figs/true_poro.png){#fig-true-poro}

Porosity inversion in case 4. (a) Initial porosity. (b) Inverted porosity. (c) Unseen ground truth porosity. The gray curve delineates the shape of the CO~2~ plume at the 25th year.
:::

## Limitations

While our case studies offer promising insights, it is crucial to acknowledge the assumptions underpinning our approach and recognize the inherent limitations that merit further investigation. Additionally, we explore the potential for integrating this 4D processing workflow with other reservoir characterization and management strategies.

### Reservoir simulation

Our study assumes known values for all multiphase flow model parameters, including relative permeability functions, residual water saturation, temperature and capillary pressure. These parameters were kept constant in the simulations to isolate the impact of permeability on seismic data, but there can be significant rock-dependent variations in practice. A multiparameter inversion, indicated by the preliminary case study in case 4, is worthwhile for future investigation to extend this inversion framework through joint estimation of these parameters. Further exploration is also required to understand the crosstalk between these parameters. In addition, our assumption that supercritical CO~2~ miscibility in the resident brine is low could be removed by considering a compositional flow model that introduces additional uncertain parameters. The feasibility of such approaches hinges on the availability of a differentiable reservoir simulator, like [JutulDarcy.jl](https://github.com/sintefmath/JutulDarcy.jl), or the use of deep neural networks to approximate the physics of multiphase flow [@grady2023model; @louboutin2023; @yin2023solving] and serve as a surrogate during inversion. Moreover, multiphase flow equations may not hold in scenarios involving CO~2~ leakage, necessitating robust leakage detection methodologies [@erdinc2022; @yin2023derisking].

### Rock physics

The case studies currently ignore the pressure effect on the wave properties [@macbeth2004; @macbeth2006]. While this can be justified for some GCS projects where the pressure change is relatively subtle, the inversion framework can be extended to honor the relationship between pressure and wave properties, and include geomechanical effects [@chen2020coupled; @nagao2023efficient]. The patchy saturation model may also not fully capture the complexities of real-world reservoirs [@allo2024overcoming], indicating a need for calibration of the rock physics model against actual reservoir and seismic data.

### Wave physics

The omission of updates to the brine-filled baseline velocity model represents a simplification that warrants further exploration. Future research could extend the framework to jointly update this baseline alongside permeability, incorporating additional parameters like shear velocity and density, which are currently ignored in the modeling and inversion. Quantifying uncertainties in velocity [@yin2024wise; @yin2024wiser; @orozco2024aspire] and permeability models remains a critical challenge for enhancing the reliability of inversion results.

## Discussion and conclusion

Our feasibility studies demonstrate the performance of this inversion framework in estimating permeability models directly from multiple prestack time-lapse seismic datasets in a cross-well setting. The recovered horizontal details in the permeability, especially shown in @fig-case2-update-perm using a $6\,\mathrm{m}$ grid spacing, highlight the capability of this end-to-end inversion framework to provide high-resolution spatial information of the permeability model. This framework also has the potential to significantly reduce cycle time in 4D processing workflows by avoiding labor-intensive, step-by-step inversion schemes where seismic velocity and CO~2~ saturation are subsequently inverted from right to left in @fig-multiphysics. Furthermore, the proposed framework differs from conventional workflows [@hatab2021assessing; @hatab2021assessment] by coupling fluid-flow, rock, and wave physics, leveraging the sensitivities of the reservoir simulator through a differential programming software framework [@louboutin2023]. By utilizing the rock physics model that links changes in CO~2~ saturation to changes in seismic properties, we gain access to these sensitivities from the time-lapse seismic data, allowing us to invert for the permeability directly.

Opportunities for future research remain. This inversion framework can be further enhanced to incorporate multimodal observations, such as a combination of well measurements and seismic data. This enhancement can be readily achieved by integrating additional misfit terms into the objective function, as detailed by @yin2023solving. Moreover, the resolution of permeability inversion deserves more exploration. Currently, the inversion is limited to the resolution achievable by seismic methods. Adopting a multimodal data assimilation approach necessitates further studies into the upscaling of the permeability model [@wu2002analysis], which we leave for future work. Additionally, the permeability inversion framework is well-suited for integration with the digital twin framework, as reported by @herrmann2023president and in other ongoing projects at our research group. When time-lapse seismic and well measurements are collected from the field, performing permeability inversion facilitates more accurate estimations of reservoir properties. These estimations can then be used to forecast the CO~2~ plume's behavior and optimize well injectivity in GCS projects, thereby maximizing injection volumes while minimizing fracturing risks [@gahlot2024digital].

## Data availability

Software with this research can be accessed at [http://github.com/slimgroup/TL-FWPI.jl/](http://github.com/slimgroup/TL-FWPI.jl/), with the DOI link [https://doi.org/10.5281/zenodo.10910283](https://doi.org/10.5281/zenodo.10910283). The full 3D Compass model is open access with CC BY license, available at [ftp://slim.gatech.edu/data/synth/Compass/](ftp://slim.gatech.edu/data/synth/Compass/).

## Acknowledgements

The authors would like to thank Zhishuai Zhang and the other anonymous reviewer for the thorough reviews and valuable suggestions to our original manuscript. This research was carried out with the support of Georgia Research Alliance, partners of the ML4Seismic Center, and the US National Science Foundation grant OAC 2203821.

The North Sea CO~2~ storage project information is taken from the Strategic UK CCS Storage Appraisal Project, funded by the UK Department of Energy and Climate Change, commissioned by the Energy Technologies Institute (ETI), and delivered by Pale Blue Dot Energy, Axis Well Technology, and Costain. The information contains copyright information licensed under ETI Open License.

During the preparation of this work, the authors used ChatGPT to refine sentence structures and improve the readability of the manuscript. After using this service, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication.