# Monitoring fluid saturation in reservoirs using time-lapse full-waveform
inversion

black Amir Mardan
INRS-ETE
Québec, QC, Canada
<EMAIL_ADDRESS>
&black Bernard Giroux
INRS-ETE
Québec, QC, Canada
<EMAIL_ADDRESS>
&black Gabriel Fabien-Ouellet
Polytechnique Montréal
Montréal, QC, Canada
<EMAIL_ADDRESS>
&black Mohammad Reza Saberi
GeoSoftware
The Hague, Netherlands
<EMAIL_ADDRESS>
###### Abstract

Monitoring the rock-physics properties of the subsurface is of great
importance for reservoir management. For either oil and gas applications or
CO2 storage, seismic data are a valuable source of information for tracking
changes in elastic properties which can be related to fluids saturation and
pressure changes within the reservoir. Changes in elastic properties can be
estimated with time-lapse full-waveform inversion (TL-FWI). Monitoring rock-
physics properties, such as saturation, with TL-FWI is usually a two-step
process: first, elastic properties are estimated with full-waveform inversion
(FWI), then, the rock-physics properties are estimated with rock-physics
inversion. However, multiparameter TL-FWI is prone to crosstalk between
parameter classes across different vintages. This leads to leakage from one
parameter class to another, which, in turn, can introduce large errors in the
estimated rock-physics parameters. To avoid inaccuracies caused by crosstalk
and the two-step inversion strategy, we reformulate TL-FWI to estimate
directly the changes in the rock-physics properties. Using Gassmann’s model,
we adopt a new parameterization containing porosity, clay content, and water
saturation (PCS). In the context of reservoir monitoring, changes are assumed
to be induced by fluid substitution only. The porosity and clay content can
thus be kept constant during time-lapse inversion. We compare the PCS
parameterization with the usual density-velocity (DV) parameterization for
different benchmark models. Results indicate that the proposed
parameterization eliminates crosstalk between parameters of different
vintages, leading to more accurate estimation of saturation changes. We also
show that using the PCS parameterization, the elastic changes can be monitored
more accurately.

_Keywords_ full-waveform inversion $\cdot$ CO2 monitoring $\cdot$ time-lapse
inversion $\cdot$ reservoir monitoring $\cdot$ rock-physics monitoring

## Introduction

Either for oil and gas applications or CO2 sequestration, the ability to
monitor reservoirs is highly important as it allows for safe and profitable
operations. Time-lapse seismic monitoring is a powerful tool that allows to
investigate large-scale changes in a reservoir Lumley, (2001); Landrø et al.,
(2003); Dupuy et al., 2016a ; Maharramov et al., (2016); Fabien-Ouellet et
al., (2017); Lang and Grana, (2019); Zhou and Lumley, (2021); Mardan et al.,
2022a ; Mardan et al., 2022c . Seismic monitoring is sensitive to changes in
elastic parameters, which can be related to reservoir parameters, such as
saturation, through a rock-physics model.

Numerous studies have been conducted to monitor saturation changes in
reservoirs. Landrø et al., (2003) used $PP$ and $PS$ seismic data to better
estimate the time-lapse changes. Buland and Ouair, (2006) proposed a Bayesian
time-lapse inversion method and applied it to monitor the Norne Field offshore
Norway with seismic data acquired from 2001 to 2003. Vedanti and Sen, (2009)
used prestack seismic data inversion to monitor the in situ combustion in the
Balol oil field. Lang and Grana, (2019) studied the potential of the Johansen
field in Norway by simulating $\text{CO}_{2}$ injection for 10 years and
monitoring the reservoir. All of these studies are based on linearized
amplitude variation with offset (AVO) and the convolutional model whose
accuracy decreases by increasing the incidence angle Mallick, (2007). Although
AVO has been very successful and is commonly used in industry, it has some
other limitations. AVO inversion only uses the amplitude of the reflected
waves while for better imaging, all type of waves should be involved during
the inversion Virieux et al., (2017). In addition, AVO relies on the
assumption that reflectors are flat, which is a restrictive hypothesis in many
cases Queißer and Singh, (2013). Last but not least, AVO suffers from the
uncertainties in data preprocessing such as velocity model errors Naeini et
al., (2017); Hu et al., (2021). These assumptions and the incomplete use of
waveforms can lead to suboptimal use of the resources invested in monitoring.

A variety of studies have used time-lapse full-waveform inversion (TL-FWI) to
monitor reservoirs. For example, Watanabe et al., (2005) used time-lapse
crosswell seismic data to monitor the changes during a gas hydrate thermal
production test at the Mallik site. Raknes et al., (2013) monitored the
leakage of one producing well over a field in North Sea by studying the
seismic data of two vintages with two years time difference. Hicks et al.,
(2016) and Lescoffit et al., (2016) showed the interest of TL-FWI when
studying seismic data and changes in $P$-wave velocity ($V_{P}$) at Grane
field due to gas-oil replacement. Maharramov et al., (2016) showed the
potential of TL-FWI to image the changes in the Genesis field in the Green
Canyon area of the central Gulf of Mexico. In these studies, TL-FWI is carried
out to map the acoustic changes in the subsurface using a monoparameter
acoustic formulation (analyzing only $V_{P}$). To obtain the changes in terms
of rock-physics properties, full-waveform inversion (FWI) is employed to
estimate $V_{P}$ and in the next step, the estimated $V_{P}$ is inverted using
rock-physics inversion. For example, this two-step procedure has been used to
estimate the changes in the saturation of CO2 at Sleipner field in North Sea
Queißer and Singh, (2013); Dupuy et al., 2016a ; Dupuy et al., 2016b ; Yan et
al., (2019); Dupuy et al., (2021); Romdhane et al., (2022).

To approach the true properties of the subsurface, an elastic formulation is
preferred over acoustic as it allows taking into account the different seismic
phases in the measured data. An elastic medium requires multiparameter TL-FWI.
In this case, model parameters can have coupled effects on the seismic data,
which is referred to as crosstalk between parameters Operto et al., (2013);
Yang et al., (2018). Ma et al., (2016) studied the Marmousi model for
multiparameter acoustic TL-FWI and showed that the crosstalk problem is more
severe in TL-FWI. Although Ma et al., (2016) recovered good images of the
subsurface using FWI, the estimated time-lapse images were inaccurate due to
the crosstalk problem (Figure 2 and 3 of Ma et al., (2016)). This problem can
be alleviated by using different parameterizations Tarantola, (1986); Operto
et al., (2013); Ma et al., (2016); Hu et al., (2021). A comprehensive
discussion is provide by Pan et al., (2019) where the authors analyzed the
efficiency of different parameterizations such as velocity-density (DV)
($P$-wave velocity, $S$-wave velocity, and density), modulus-density (DM)
(bulk modulus, shear modulus, and density), impedance-density (D-IP) ($P$-wave
impedance, $S$-wave impedance, and density), velocity-impedance-$I$ (V-IP-$I$)
($P$-wave velocity, $S$-wave velocity, and $P$-wave impedance), and velocity-
impedance-$II$ (V-IP-$II$) ($P$-wave velocity, $S$-wave velocity, and $S$-wave
impedance) for performing elastic FWI.

The goal of this study is to recover the changes in water saturation ($S_{w}$)
by considering an elastic earth. To minimize the effects of crosstalk that can
occur due to the time-lapse nature of the problem, the TL-FWI is formulated in
terms of porosity ($\phi$), clay content ($C$), and water saturation, as
proposed by Hu et al., (2021) in the context of FWI. This parameterization
(called PCS) is chosen to facilitate the time-lapse inversion as we can assume
that porosity and clay content are constant with time in the reservoir. To
perform this technique, a connection between the elastic and rock-physical
properties is required. In this work, we use Gassmann’s equation for this
purpose. Using the PCS parameterization, we combine the rock-physics inversion
and FWI in a single inversion scheme. In this way, saturation can be directly
inverted from the raw seismic data (shot gather). By performing TL-FWI with
the PCS parameterization, we can also easily obtain the time-lapse images of
elastic properties using the underlying rock-physics model.

This paper is organized by first presenting FWI and the required formulation
to include Gassmann’s model for estimating the porosity, clay content, and
water saturation. This section is followed by a brief discussion on TL-FWI.
Then, we study and discuss the efficiency of this method by using numerical
models.

## Theory

### Full-waveform inversion

Full-waveform inversion is a local optimization process that minimizes
residuals between the observed and estimated wavefields at the receivers
locations for different source positions Virieux and Operto, (2009). The
objective function can be written as

$\chi(\mathbf{m})=\frac{1}{2}\|\mathbf{W}_{d}\left(\mathbf{R}\mathbf{u}(\mathbf{m})-\mathbf{d}\right)\|_{2}^{2}+\chi_{REG},\
$ (1)

where $\mathbf{R}$ maps the wavefield ($\mathbf{u}$) to the receivers
locations. Vectors $\mathbf{m}$ and $\mathbf{d}$ are the model parameter and
observed data, respectively. $\mathbf{W}_{d}$ is a weighting operator on the
data misfit and $\chi_{REG}$ is a regularization term. FWI is highly ill-
posed, where an infinite number of models matches the data. A variety of
regularization methods such as the Tikhonov, total variation, prior
information, and parameter-relation techniques have been proposed to condition
the FWI problem and reduce the ill-posedness Virieux and Operto, (2009);
Asnaashari et al., (2015). In this study, we use the Tikhonov regularization.
This regularization method is presented in Appendix A in addition to total
variation, prior information, and parameters-relation techniques. In equation
1, the wavefield is the solution of a partial differential equation (PDE)
which is discretized to perform the forward modeling. For the 2D isotropic
elastic case in time domain, this PDE can be written as

$\begin{cases}\rho\dot{v}_{x}=\partial_{x}\tau_{xx}+\partial_{z}\tau_{xz},\\\
\\\ \rho\dot{v}_{z}=\partial_{z}\tau_{zz}+\partial_{x}\tau_{xz},\\\ \\\
\dot{\tau}_{xx}=(\lambda+2G)\partial_{x}v_{x}+\lambda\partial_{z}v_{z}+s,\\\
\\\
\dot{\tau}_{zz}=\lambda\partial_{x}v_{x}+(\lambda+2G)\partial_{z}v_{z}+s,\\\
\\\ \dot{\tau}_{xz}=G(\partial_{z}v_{x}+\partial_{x}v_{z}),\\\ \end{cases}$
(2)

where $\lambda$ and $G$ denote lamé’s parameter and where $\rho$ and $s$ are
density and the source function. In equation 2, the particle velocities in
$x$\- and $z$-directions ($v_{x}$ and $v_{z}$) as well as ($\tau_{xx}$ and
$\tau_{zz}$) and shear stresses ($\tau_{xz}$) constitute vector $\mathbf{u}$.
Spatial partial derivative operators are indicated by $\partial_{x}$ and
$\partial_{z}$ for $x-$ and $z-$directions and the temporal differentiation is
denoted by overhead dot ($\dot{\ooalign{$\sqsubset\mkern 3.0mu$\cr$\mkern
3.0mu\sqsupset$\cr}}$).

In the case of equation 2, the vector of model parameter, $\mathbf{m}$ in
equation 1, is

$\mathbf{m}=[\lambda,G,\rho]^{T},$ (3)

where T is the transpose operator. Although, the forward modeling is performed
using Lamé parameters and density, FWI can be performed using different
parameterizations. For example, the DV parameterization is a common
formulation. Switching between parameters can be carried out using the chain
rule, as shown at the end of this section.

#### Optimization

To minimize the cost function (equation 1), a variety of optimization
algorithms such as steepest-descent, conjugate-gradient, $\ell$-BFGS, and
Newton algorithms can be employed Virieux and Operto, (2009). These algorithms
minimize the cost function in the vicinity of a starting model,
$\mathbf{m}_{0}$, as

$\mathbf{m}=\mathbf{m}_{0}+\alpha\boldsymbol{\Delta}\mathbf{m},$ (4)

where $\alpha$ is step size and the search direction,
$\boldsymbol{\Delta}\mathbf{m}$, for gradient descent is

$\boldsymbol{\Delta}\mathbf{m}=-\left[\begin{array}[]{c}\nabla_{\mathbf{m}^{(1)}}\chi\\\
\nabla_{\mathbf{m}^{(2)}}\chi\\\ \nabla_{\mathbf{m}^{(3)}}\chi\\\
\end{array}\right],$ (5)

where $\nabla_{\mathbf{m}^{(n)}}\chi$ is the gradient of the cost function
with respect to model parameter $\mathbf{m}^{(n)}$.

One of the most problematic challenges of multiparameter FWI is crosstalk
between parameters. Crosstalk is due to the fact that different parameters can
have a coupled effect on the seismic wavefield and different combination of
changes in these parameters can cause similar seismic response Operto et al.,
(2013); Yang et al., (2018). To address this problem, different strategies
have been deployed (Métivier et al.,, 2013; Operto et al.,, 2013; Lavoué et
al.,, 2014; Fabien-Ouellet et al.,, 2017; Keating and Innanen,, 2019). Operto
et al., (2013) listed four main strategies: the choice of appropriate
parameterization, the use of the Hessian as it can measure the level of
coupling between parameters, a data driven strategy i.e., using multi-
components data, and a sequential inversion where the dominant parameters are
estimated before the secondary parameters. To further reduce the crosstalk
problem, it is also proposed to scale the gradient Kamei and Pratt, (2013);
Lavoué et al., (2014). Hence equation 5 can be rewritten as

$\boldsymbol{\Delta}\mathbf{m}=-\left[\begin{array}[]{c}\nabla_{\mathbf{m}^{(1)}}\chi\\\
\xi\nabla_{\mathbf{m}^{(2)}}\chi\\\ \beta\nabla_{\mathbf{m}^{(3)}}\chi\\\
\end{array}\right],$ (6)

where $\xi$ and $\beta$ scale the gradient of different parameter classes.

By using a quasi-Newton method, $\ell$-BFGS, we consider an approximation of
the Hessian into the solution (Nocedal and Wright,, 2006) which helps reducing
crosstalk. However the most important point of this study is the chosen
parameterization. Considering the fact that the time-lapse changes in $P$-wave
velocity due to fluid substitution (more affected than $S$-wave velocity and
density) are approximately $10\%$ of the baseline model Zhou and Lumley,
(2021), the challenges of FWI are more pronounced in time-lapse studies than
in conventional FWI. Ma et al., (2016) studied the possibility of deploying
TL-FWI for multiparameter acoustic studies. Although they estimate the model
parameters very well, crosstalk prevents from obtaining an acceptable time-
lapse image of the subsurface for any of two time-variable parameters.
Thereby, an efficient parameterization is required to minimize the crosstalk
in multiparameter TL-FWI. Inspired by Hu et al., (2021), we study the
feasibility of performing TL-FWI using the PCS parameterization. Recalling
that the porosity and clay content can be reasonably considered constant with
time, water saturation is the only time-variable parameter.

In the next section, we provide the required equations to adapt the FWI
problem to estimate the PCS parameters using Gassmann’s equation. In this
case, the vector of model parameter can be written as

$\mathbf{m}=[\phi,C,S_{w}]^{T}.$ (7)

Using the aforementioned assumption that porosity and clay content remain
constant over time for reservoirs, the changes in water saturation can be the
only parameter to be sought in the time-lapse inversion.

It is worth recalling that the chain rule is used to change the parameters for
inversion. As an example, $\mathbf{q}$ parameterization,
$\mathbf{q}=[q_{1},q_{2},q_{3}]^{T}$, can be changed to a new
parameterization, $\mathbf{p}=[p_{1},p_{2},p_{3}]^{T}$, as

$\frac{\partial\chi}{\partial q_{i}}=\frac{\partial\chi}{\partial
p_{1}}\frac{\partial p_{1}}{\partial q_{i}}+\frac{\partial\chi}{\partial
p_{2}}\frac{\partial p_{2}}{\partial q_{i}}+\frac{\partial\chi}{\partial
p_{3}}\frac{\partial p_{3}}{\partial q_{i}},$ (8)

for $i\in(1,2,3)$. FWI can be performed in any parameterization while the
forward modeling is implemented using equation 2.

This study is carried out using PyFWI Mardan et al., 2023a  which is an open-
source package that implements FWI in terms of DV parameterization. In the
next section, we present the methodology to transfer the earth model from PCS
parameterization to DV parameterization for forward modeling. In Appendix B,
we present the required equations to obtain the gradient of the cost function
in terms of PCS parameterization for performing the inversion.

### Porous media homogenization

To perform TL-FWI with a rock-physics parameterization, a rock-physics model
is needed to map the rock-physics properties to the elastic properties.
Although there are numerous empirical rock-physics models that satisfy this
requirement, we employ Gassmann’s equation (Gassmann,, 1951) which is valid
for most of consolidated rocks (Dupuy et al., 2016b, ). Gassmann’s relation is
used to compute the variation of bulk modulus during the fluid substitution.
In practice, Gassmann’s model requires some information on the fluid and
mineral properties. In the context of TL-FWI, this information should be
available as time-lapse FWI is usually performed in well studied fields.

The effective properties of the fluid (bulk modulus, $K$, and density) can be
calculated using weighted arithmetic averages (Voigt,, 1889),

$\displaystyle K_{f}$ $\displaystyle=S_{w}K_{w}+(1-S_{w})K_{h},$ (9)
$\displaystyle\rho_{f}$ $\displaystyle=S_{w}\rho_{w}+(1-S_{w})\rho_{h},$

where the subscripts $f$, $w$, and $h$ respectively denote the effective
property of the fluid, water, and hydrocarbon.

The density of the solid frame of rock can be calculated using weighted
average of the density of the grains as

$\displaystyle\rho_{s}=$ $\displaystyle C\rho_{c}+(1-C)\rho_{q},$ (10)

where subscripts $s$, $c$, and $q$ denote solid medium, clay, and quartz,
respectively. The bulk and shear moduli of the solid can be estimated using
Voigt-Reuss-Hill method

$\displaystyle K_{s}=$
$\displaystyle\frac{1}{2}\left[\left(CK_{c}+(1-C)K_{q}\right)\right.$ (11)
$\displaystyle+$
$\displaystyle\left.\left(\frac{1}{C/K_{c}+(1-C)/K_{q}}\right)\right],$
$\displaystyle G_{s}=$
$\displaystyle\frac{1}{2}\left[\left(CG_{c}+(1-C)G_{q}\right)\right.$
$\displaystyle\left.+\left(\frac{1}{C/G_{c}+(1-C)/G_{q}}\right)\right].$

The effective dry moduli of the rock, $K_{D}$ and $G_{D}$ can be computed
using various effective medium theories Berryman, (1995); Mavko et al.,
(2020). Following Pride, (2005), we estimate the dry moduli with

$\displaystyle K_{D}$ $\displaystyle=K_{s}\frac{1-\phi}{1+cs\phi},$ (12)
$\displaystyle G_{D}$ $\displaystyle=G_{s}\frac{1-\phi}{1+\frac{3}{2}cs\phi},$

where $cs$ is the general consolidation parameter which defines the degree of
consolidation among grains. This parameter can be considered as a free
parameter in predicting velocities Lee, (2005). For further studies about
$cs$, readers are referred to Lee, (2005) where the author discusses the
importance of this parameter and how it can be calculated.

Finally, the effective properties of the porous medium is calculated with the
following expressions Gassmann, (1951); Voigt, (1889),

$\displaystyle K=\frac{\phi
K_{D}+\left(1-(1+\phi)K_{D}/K_{s}\right)K_{f}}{\phi(1+\Delta)},$ (13)
$\displaystyle G=G_{D},$ $\displaystyle\rho=(1-\phi)\rho_{s}+\phi\rho_{f},$

where

$\Delta=\frac{1-\phi}{\phi}\frac{K_{f}}{K_{s}}\left(1-\frac{K_{D}}{(1-\phi)K_{s}}\right)=\frac{1-\phi}{\phi}\frac{K_{f}}{K_{s}}\left(1-\frac{1}{1+cs\phi}\right).$
(14)

The elastic properties can be calculated using equation 13 as

$\displaystyle V_{P}=\sqrt{\frac{K+\frac{4}{3}G}{\rho}},$ (15) $\displaystyle
V_{S}=\sqrt{\frac{G}{\rho}},$ $\displaystyle\rho=\rho.$

For a detailed discussion on Gassmann’s model and the provenance of equations
9-15 see Dupuy et al., 2016b .

### Time-lapse full-waveform inversion

Time-lapse FWI can be performed by considering different combinations of FWI
runs. These combinations can be performed sequentially such as for the cases
of the cascaded, cross-updating, and weighted-average methods Watanabe et al.,
(2005); Maharramov et al., (2016); Mardan et al., 2022b ; Mardan et al., 2023b
, or in parallel as in the case of the independent and central-difference
methods Plessix et al., (2010); Zhou and Lumley, (2021). In addition to these
methods, Maharramov et al., (2016) proposed a joint inversion of seismic data
from different vintages by introducing the cost function as

$\begin{split}\chi_{tl}(\mathbf{m}_{b},\mathbf{m}_{m})=&\|\mathbf{Ru}(\mathbf{m}_{b})-\mathbf{d}_{b}\|_{2}^{2}+\|\mathbf{Ru}(\mathbf{m}_{m})-\mathbf{d}_{m}\|_{2}^{2}\\\
&+\delta\|\mathbf{m}_{m}-\mathbf{m}_{b}\|_{2}^{2},\end{split}$ (16)

where subscripts $b$ and $m$ denote baseline and monitor models. The
difference between the estimated baseline and monitor models is used as a
regularization term and the weight of this term in the cost function is
controlled by parameter $\delta$ which is constant during the inversion. This
parameter, $\delta$, can vary for different surveys and should be picked by
trial and error. More details on this parameter can be found in Maharramov et
al., (2016). The flowchart of this method (hereafter called simultaneous TL-
FWI) is presented in Figure 1. Simultaneous TL-FWI gets an initial model and
observed seismic data from two vintages. In the first step, the baseline data
are inverted and the estimated baseline is used as the initial model for time-
lapse FWI. TL-FWI minimizes the residuals between the estimated and observed
data from the baseline and monitor models regularized by the difference
between the estimated models.

Figure 1: Flowchart of simultaneous TL-FWI. After performing FWI on baseline
data, a joint time-lapse FWI inverts baseline and monitor data where the
inversion process is regularized by the difference between estimated baseline
and monitor models.

Simultaneous TL-FWI is one of the most accurate strategies for time-lapse
studies (Mardan et al., 2023b, ). Incorporating a regularization term for
minimizing the difference between two estimates (baseline and monitor models)
helps this method to produce less artifacts.

The methodology we propose is as follows. We first estimate the baseline
model, using standard FWI in the PCS parametrization. Then we perform a
single-parameter inversion of $S_{w}$ with jointly inverting the baseline and
monitor data (equation 16). We justify the last step by making the hypothesis
that porosity and clay content remains constant after fluid injection. In the
case of inversion with the DV parameterization, the TL-FWI is a multi-
parameter problem. Furthermore, it requires the additional step of inverting
the elastic properties to obtain $S_{w}$. In this work, we use $\ell$-BFGS to
do so, which minimizes the following cost function

$\chi_{rp}=\frac{1}{2}\|F(\phi,C,S_{w})-\mathbf{m}_{DV}\|^{2},$ (17)

where $F$ is a function that calculates $V_{P}$, $V_{S}$, and $\rho$ based on
equations 9-15 and $\mathbf{m}_{DV}$ is vector of estimated elastic properties
using the DV parameterization.

## Numerical analysis

In this section we assess the feasibility of the proposed strategy with two
synthetic models. For both analysis, FWI is performed to estimate all model
properties (elastic or rock-physics properties) simultaneously. In this
regard, the multi-scale strategy Bunks et al., (1995) is employed to avoid
cycle-skipping. We filter the observed and estimated data using a lowpass
Butterworth filter with cut-off frequencies of 10, 20, 25, 35, 40, and 55 Hz.
The source for forward modeling is a Ricker wavelet with a central frequency
of 30 Hz. Perfectly matched layers (PML) are used at the boundaries of the
models to minimize waves reflected from the boundaries Berenger, (1994). We
use the Tikhonov regularization (equation 18) in the $x$-direction
($\mathbf{B}_{z}=0$) to decrease the ill-posedness of FWI problem and make the
estimated models smoother in the $x$-direction.

For numerical studies, we use two synthetic models. For the first experiment,
a simple layered model is used. Then, TL-FWI is performed using the Marmousi
model which is a more realistic case. During the inversion, the gradient is
scaled using equation 6. The scaling weights $\xi$ and $\beta$ are obtained by
trial and error, by performing FWI runs on the simple layered model. Then we
use the same weights for the Marmousi model. The rock-physical properties used
in this study are shown in Table 1.

| Bulk modulus (GPa) | Shear modulus (GPa) | Density (g/cm3)
---|---|---|---
Quartz | 37.00 | 44.00 | 2.65
Clay | 21.00 | 10.0 | 2.55
Water | 2.25 | 0.00 | 1.00
Gas | 0.04 | 0.00 | 0.10
Table 1: Elastic properties in this study for minerals and fluids where $cs$
is set to 20 Dupuy et al., 2016a ; Hu et al., (2021).

### Layered model

To investigate the performance of the methodology, a simple model with three
flat reflectors is first used (Figure 2). Figure 2a-2c shows the true models
for porosity, clay content and water saturation for the baseline and Figure
2d-2f shows the true monitor model. The initial model of the inversion is
shown in Figure 2g-2i. Finally, the true time-lapse model is shown in Figure
2j-2l. The elastic properties of this model are shown in Figure 3. Comparing
Figure 2j-2l with Figure 3j-3l shows that with a $25\%$ increment of water
saturation, P-wave velocity and density increase by $3\%$ and $0.8\%$,
respectively, while the S-wave velocity decreases by $0.4\%$.

Figure 2: Rock-physics properties of the layered model. (a-c) True baseline,
(d-f) true monitor, and (g-i) initial models. (j-l) Time-lapse changes between
porosity (a and d), clay content (b and e), and water saturation (c and f).
The time-lapse model is presented with two color scales which show the changes
in value and percentage. Dashed line in (l) is used for 1D profile in Figure 5
and 8.  Figure 3: Elastic properties of the layered model. (a-c) True
baseline, (d-f) true monitor, and (g-i) initial models. (j-l) Time lapse
changes between $V_{P}$ (a and d), $V_{S}$ (b and e), and $\rho$ (c and f).
The time-lapse model is presented with two color scales which show the changes
in value and percentage.

To perform this study, seven isotropic sources are used on the surface for
forward modeling. Receivers are located on the surface and in two imaginary
wells at both sides of the model. Noise-free pressure data are inverted to
recover the time-lapse anomaly. The rock-physics properties estimated by using
Gassmann’s equation are presented in Figure 4a-4c. The parameters in this
model are used to obtain an indirect estimate of $P$\- and $S$-wave velocities
as well as density, as shown in Figure 4d-4f. A 1D plot obtained at $X=0.5$ km
(dashed line in Figure 2l) is also shown in Figure 5 for both rock-physics
properties, Figure 5a-5c, and elastic properties, Figure 5d-5f.

Figure 4: The estimated baseline model of (a) porosity, (b) clay content, (c)
water saturation, (d) $P$-wave velocity, (e) $S$-wave velocity, and (f)
density estimated using the PCS parameterization.  Figure 5: 1D assessment of
the results for (a) porosity, (b) clay content, (c) water saturation, (d)
$P$-wave velocity, (e) $S$-wave velocity, and (f) density at $X=0.5$ km
(center of the 2D section). True value is shown with solid line, estimated
model is shown with dashed line and dotted line presents the initial model.

Results of the first step obtained with this simple synthetic model shows that
FWI and the proposed formulation have the potential to recover the rock-
physics properties. However, reconstructing the clay content seems problematic
and the interfaces in this parameter are not recovered sharply (Figure 4b and
5b). Although Hu and Innanen, (2021) showed that the $C$ model can be
estimated better by considering a parameter-relation regularization, this type
of regularization is not taken into account in this study. The parameter-
relation regularization allows us to regularize the inversion based on the
relation between two parameters in the field. So in general, we can expect
that the parameter-relation and prior-model regularization Asnaashari et al.,
(2013) should improve the accuracy of the estimated clay content. These two
regularization methods are provided in Appendix A for interested readers.

After estimating the baseline model (Figure 4), TL-FWI can be performed using
equation 16. We first compare the efficiency of the DV and PCS
parameterization to monitor the rock-physics properties. Using the DV
parameterization, changes in terms of the rock-physics properties are
estimated in two steps. We first estimate the elastic properties of the
baseline and the monitor models using equation 16. In the second step,
equation 17 is employed to estimate the rock-physics properties of the
baseline and the monitor models and differences are presented in Figure 6a-6c.
This is a naive method to perform multiparameter TL-FWI and leads to
overestimation of porosity and inaccurate time-lapse estimate. We cannot get
an accurate time-lapse image using this method due to the crosstalk between
parameters. It is worth reminding that the crosstalk problem is more severe in
time-lapse inversion, because the magnitude of changes in reservoirs is
comparable to generated artifacts. To improve the results, we follow another
strategy. In this strategy, the rock-physics properties of the baseline are
estimated using equation 17. Then, by assuming that porosity and clay content
are not variable with time, the elastic properties of the monitor model are
inverted to estimate the saturation. For this purpose, the estimated rock-
physics model of the baseline is used as initial model to estimate the rock-
physics properties of the monitor model. Despite the inversion of the baseline
model, we perform single-parameter inversion for estimating the rock-physics
properties of the monitor model during which only $S_{w}$ is updated. This
strategy leads to the time-lapse estimate presented in Figure 6d-6f where a
better time-lapse image of the saturation is obtained. This strategy is used
hereafter in this study to compute the time-lapse rock-physics properties
using the DV parameterization. By using the PCS parameterization and
considering time-independent porosity and clay content, direct TL-FWI is
performed to only monitor the water saturation. The estimated time-lapse
models of porosity, clay content, and water saturation are presented in Figure
6g-6i. Comparing Figure 6 to Figure 2j-2l shows that the PCS parameterization
can improve the accuracy of time-lapse studies for monitoring the saturation
of the fluids in the subsurface. The accuracy of the results is measured using
root-mean-square error (RMSE) and direct monitoring of $S_{w}$ using the PCS
parameterization leads to a $37.3\%$ higher accuracy in comparison to using
the DV parameterization.

Figure 6: Estimated time-lapse model for $\phi$, $C$, and $S_{w}$ with (a-c)
the DV parameterization while all model parameters can vary with time, (d-f)
the DV parameterization while $S_{w}$ is the only time-dependent parameter,
and (g-i) the PCS parameterization.

In addition to the rock-physics properties, the PCS parameterization can be
used to recover time-lapse images of the elastic properties. For this case, we
use equation 16 to obtain the baseline and monitor models in terms of
porosity, clay content and water saturation. The results are then used for
rock-physics modeling (equation 9-15) to obtain the elastic properties of the
baseline and monitor models. While the DV parameterization leads to time-lapse
images presented in Figure 7a-7c, the PCS parameterization improves the time-
lapse estimates (Figure 7d-7f). The DV parameterization does not allow
recovering the changes appropriately. This can be explained by considering
that $25\%$ changes in the saturation (Figure 2l) is now spread out between
three elastic properties and it causes only $3\%$, $0.4\%$ and $0.8\%$
variation respectively in $V_{P}$, $V_{S}$, and $\rho$ (Figure 3j-3l).
Comparing the estimated time-lapse changes in Figure 7 shows that using the
PCS parameterization, we can recover the elastic properties more accurately
than using the DV parameterization. The time-lapse elastic properties
estimated using the PCS parameterization are $53\%$ more accurate than with
the DV parameterization.

Figure 7: Estimated time-lapse model for $V_{P}$, $V_{S}$, and $\rho$ with
(a-c) the DV parameterization and (d-f) the PCS parameterization.

Figure 8a-8c shows 1D profiles of the true and estimated changes in the rock-
physics properties along the dashed line in Figure 2l. The estimated time-
lapse elastic properties from PCS and DV parameterizations are also compared
in Figure 8d-8f. This comparison shows that the results obtained using the DV
parameterization (dotted lines) are far from accurate. For 1D analysis, the
PCS parameterization leads to $43\%$ higher accuracy (Figure 8a-8c). For
monitoring the elastic properties, the PCS parameterization yields $57\%$
higher accuracy (Figure 8d-8f).

Figure 8: 1D assessment of the time-lapse estimates for (a) porosity, (b) clay
content, (c) water saturation, (d) $P$-wave velocity, (e) $S$-wave velocity,
and (f) density at $X=0.5$ km. True value is shown with solid line, estimated
model with PCS parameterization is shown with dashed line and dotted line
presents the estimated model with the DV parameterization.

### Marmousi model

In this section, we aim at assessing the methodology using a more realistic
model. For this purpose, a portion of the Marmousi model is chosen (Figure 9).
The true baseline and monitor models are shown in Figure 9a-9f. This model
contains two sandstone reservoirs. The monitor model (Figure 9d-9f) is created
by reducing the water saturation in the reservoir indicated by line A1 and
increasing the water saturation in the reservoir indicated by line A2 (Figure
9l).

The elastic properties of this model are shown in Figure 10. Basically, the
density of the saturated rock increases by replacing the hydrocarbon by water.
Although the $P$-wave velocity has a reverse relationship with density, the
hydrocarbon-water replacement increases the $P$-wave velocity due to the
increase of the saturated bulk modulus. However, the fluid shear modulus is
nil and $S$-wave velocity is decreased due to the increase in density. Figure
9j-9l and Figure 10j-10l show the time-lapse changes in rock-physics and
elastic properties, respectively. The relative changes are shown in Figure
11a-11c and Figure 11d-11f for the rock-physics and elastic properties. In
this case, the maximum changes in water saturation occurs in the reservoir A2,
which shows a $40\%$ change with respect to the baseline. This change leads to
maximum change of $6.2\%$, $2.31\%$, and $4.46\%$ in elastic properties,
$V_{P}$, $V_{S}$, and $\rho$ respectively. Hence, we can expect that changes
can be better tracked using the PCS parameterization.

Figure 9: Rock-physics properties of the Marmousi model. (a-c) True baseline,
(d-f) true monitor, and (g-i) initial models. (j-l) Time-lapse changes between
porosity (a and d), clay content (b and e), and water saturation (c and f).
Lines A1 and A2 are used for 1D studies in Figure 14.  Figure 10: Elastic
properties of the Marmousi model. (a-c) True baseline, (d-f) true monitor, and
(g-i) initial models. (j-l) Time-lapse changes between $V_{P}$ (a and d),
$V_{S}$ (b and e), and $\rho$ (c and f).  Figure 11: Absolute value of percent
changes relative to the baseline in terms of studied (a-c) rock-physics and
(d-f) elastic properties.

FWI of the baseline data is carried out using the PCS parameterization with
values of respectively 5 and 8 for $\xi$ and $\beta$, obtained from the
layered model, studied in the previous section. The initial model for this
inversion is presented in Figure 9g-9i. The results of the inversion of the
baseline are presented in Figure 12. As the receivers are located only on the
surface for this model, it can be seen that the quality of the model update
decreases toward the edges. Finally, the results of TL-FWI are shown in Figure
13. As shown in Figure 13a-13c, the DV parameterization allows recovering
changes in only reservoir A1. On the other hand, the time-lapse water
saturation anomalies are detected in both reservoirs using the PCS
parameterization, as shown in Figure 13d-13f.

Figure 13g-13l shows the changes in elastic properties. In this case, the PCS
parameterization, shown in Figure 13j-13l, yields significant improvement over
the DV parameterization, shown in Figure 13g-13i. This improvement leads to
$11.7\%$ higher accuracy by using the PCS parameterization. As can be seen in
Figure 13g, the changes in $P$-wave velocity are not recovered in reservoir
$A_{2}$. Although the time-lapse changes in $S$-wave velocity can be detected,
the value in reservoir $A_{1}$ is overestimated (Figure 13h). Last but not
least, the time-lapse anomalies of density are not detectable from Figure 13i.
Despite inaccuracies in the estimated time-lapse elastic properties, Figure
13j-13l show that the PCS parameterization can lead to a better estimation of
elastic properties.

Figure 12: The estimated models of (a) porosity, (b) clay content, (c) water
saturation, (d) $P$-wave velocity, (e) $S$-wave velocity, and (f) density
estimated using the PCS parameterization.  Figure 13: Estimated time-lapse
model for $\phi$, $C$, and $S_{w}$ with (a-c) the DV and (d-f) PCS
parameterizations. Estimated time-lapse model for $V_{P}$, $V_{S}$, and $\rho$
with (g-i) the DV and (j-l) the PCS parameterizations.

Finally, a 1D assessment of the results is carried out along lines A1 and A2,
shown in Figure 9l. The 1D presentation of time-lapse changes (Figure 14)
shows that TL-FWI with the DV parameterization is not able to recover the
changes in the reservoirs appropriately, while using the PCS parameterization
can provide accurate time-lapse images of the model for both the elastic and
rock-physics properties. In this 1D analysis, the PCS parameterization leads
to $24.51\%$ higher accuracy for estimating the elastic properties and
$19.57\%$ for estimating the rock-physics properties.

Figure 14: 1D assessment of the time-lapse estimates for saturation along the
line (a) A1, (b) A2. Time-lapse estimates for elastic properties along the
line (c-e) A1 and (f-h) A2. True value is shown with solid line, estimated
model with the PCS parameterization is shown with dashed line and the dotted
line presents the estimated model with the DV parameterization.

## Discussion

The potential of the PCS parameterization for estimating time-lapse changes in
elastic and rock-physics properties is shown in this study. The direct
estimation of changes in rock-physics properties provides more accurate
results rather than the indirect estimation. However, the estimates presented
in Figure 6 show clearly that an effective inversion strategy and optimization
method can improve the accuracy of estimated changes using indirect inversion.
As the computation time of the rock-physics modeling is negligible, global
inversion methods such as Monte-Carlo methods might be a better choice for the
second step of the inversion (Dupuy et al., 2016b, ). However, the direct
estimation of rock-physics properties aid to avoid the challenges of picking
an appropriate optimization method and parameters for the rock-physics
inversion.

The estimated time-lapse changes for the two presented synthetic models show
that using the PCS parameterization can also enhance the accuracy for
monitoring the elastic properties (Figure 7 and Figure 13). This is due to the
fact that using the PCS parameterization, the crosstalk between parameters of
different vintages can be avoided. In this study, we assumed that the only
variable parameter with time is fluid saturation. However, pore pressure is
another important parameter that should be considered. Although the changes in
the pore pressure can be negligible in some fields due to the properties such
as high permeability (Dupuy et al.,, 2021), further studies should be
conducted to analyze the efficiency of simultaneous monitoring of saturation
and pore pressure.

In this work, we have assumed a simple rock-physics model in which only
saturation effects on elastic properties are modeled. In reality, this model
will contain errors in describing the actual subsurface properties either from
the calibration and subsurface assumptions or even the model itself. However,
time-lapse studies are carried out in well studied fields. The calibration of
the rock-physics model should not pose too much of a problem in this scenario.
Besides, the geological information and well-log data of these fields can be
employed to regularize the inversion problem and increase the accuracy of FWI
results.

## Conclusions

Time-lapse seismic data have been proven as an important source of information
for reservoir monitoring. Full-waveform inversion is a powerful tool that
relies on seismic data to recover the subsurface properties. However this
method generally suffers from crosstalk between parameters in multiparameter
studies. This problem becomes more dramatic in time-lapse studies, because the
time-lapse changes can be caused by any combination of different parameters in
different vintages.

In this study, we formulated time-lapse full waveform inversion in terms of
porosity, clay content, and water saturation using Gassmann’s equation, to
assess its performance for reducing crosstalk and improve estimates of
saturation changes. Using this formulation, we can rely on the assumption that
water saturation is the only variable parameter that changes in time and
thereby, it is possible to avoid crosstalk caused by time-lapse study. We
showed that the PCS parameterization improves the time-lapse estimate. By
using Gassmann’s equation and the new parameterization, it is also shown that
the time-lapse estimate in elastic properties can be obtained indirectly with
higher accuracy ($11.7\%$ higher accuracy in case of the Marmousi model).

## References

  * Anagaw and Sacchi, (2012) Anagaw, A. Y., and M. D. Sacchi, 2012, Edge-preserving seismic imaging using the total variation method: Journal of Geophysics and Engineering, 9, 138–146.
  * Asnaashari et al., (2013) Asnaashari, A., R. Brossier, S. Garambois, F. Audebert, P. Thore, and J. Virieux, 2013, Regularized seismic full waveform inversion with prior model information: Geophysics, 78, R25–R36; doi: 10.1190/geo2012-0104.1.
  * Asnaashari et al., (2015) ——–, 2015, Time-lapse seismic imaging using regularized full-waveform inversion with a prior model: Which strategy?: Geophysical Prospecting, 63, 78–98; doi: 10.1111/1365-2478.12176.
  * Berenger, (1994) Berenger, J.-P., 1994, A perfectly matched layer for the absorption of electromagnetic waves: Journal of Computational Physics, 114, 185–200; doi: 10.1006/jcph.1994.1159.
  * Berryman, (1995) Berryman, J. G., 1995, Mixture theories for rock properties: Rock physics and phase relations: A handbook of physical constants, 3, 205–228.
  * Buland and Ouair, (2006) Buland, A., and Y. E. Ouair, 2006, Bayesian time-lapse inversion: Geophysics, 71, R43–R48; doi: 10.1190/1.2196874.
  * Bunks et al., (1995) Bunks, C., F. M. Saleck, S. Zaleski, and G. Chavent, 1995, Multiscale seismic waveform inversion: Geophysics, 60, 1457–1473; doi: 10.1190/1.1443880.
  * (8) Dupuy, B., S. Garambois, A. Asnaashari, H. M. Balhareth, M. Landrø, A. Stovas, and J. Virieux, 2016a, Estimation of rock physics properties from seismic attributes — part 2: Applications: Geophysics, 81, M55–M69; doi: 10.1190/geo2015-0492.1.
  * (9) Dupuy, B., S. Garambois, and J. Virieux, 2016b, Estimation of rock physics properties from seismic attributes — part 1: Strategy and sensitivity analysis: Geophysics, 81, M35–M53; doi: 10.1190/geo2015-0239.1.
  * Dupuy et al., (2021) Dupuy, B., A. Romdhane, P. Eliasson, and H. Yan, 2021, Combined geophysical and rock physics workflow for quantitative CO2 monitoring: International Journal of Greenhouse Gas Control, 106, 103217; doi: 10.1016/j.ijggc.2020.103217.
  * Esser et al., (2018) Esser, E., L. Guasch, T. van Leeuwen, A. Y. Aravkin, and F. J. Herrmann, 2018, Total variation regularization strategies in full-waveform inversion: SIAM Journal on Imaging Sciences, 11, 376–406; doi: 10.1137/17M111328X.
  * Fabien-Ouellet et al., (2017) Fabien-Ouellet, G., E. Gloaguen, and B. Giroux, 2017, Time domain viscoelastic full waveform inversion: Geophysical Journal International, 209, 1718–1734; doi: 10.1093/gji/ggx110.
  * Gassmann, (1951) Gassmann, F., 1951, Über die elastizität poröser medien: Vierteljahrsschrift der Naturforschenden Gesellschaft in Zurich, 96, 1–23.
  * Hicks et al., (2016) Hicks, E., H. Hoeber, M. Houbiers, S. P. Lescoffit, A. Ratcliffe, and V. Vinje, 2016, Time-lapse full-waveform inversion as a reservoir-monitoring tool — A North Sea case study: The Leading Edge, 35, 850–858; doi: 10.1190/tle35100850.1.
  * Hu and Innanen, (2021) Hu, Q., and K. Innanen, 2021, Elastic full-waveform inversion with rock-physics constraints: First International Meeting for Applied Geoscience & Energy, Society of Exploration Geophysicists, 662–666.
  * Hu et al., (2021) Hu, Q., S. Keating, K. A. Innanen, and H. Chen, 2021, Direct updating of rock-physics properties using elastic full-waveform inversion: Geophysics, 86, MR117–MR132; doi: 10.1190/geo2020-0199.1.
  * Kamei and Pratt, (2013) Kamei, R., and R. G. Pratt, 2013, Inversion strategies for visco-acoustic waveform inversion: Geophysical Journal International, 194, 859–884; doi: 10.1093/gji/ggt109.
  * Keating and Innanen, (2019) Keating, S., and K. A. Innanen, 2019, Parameter crosstalk and modeling errors in viscoacoustic seismic full-waveform inversion: Geophysics, 84, R641–R653; doi: 10.1190/geo2018-0410.1.
  * Landrø et al., (2003) Landrø, M., H. H. Veire, K. Duffaut, and N. Najjar, 2003, Discrimination between pressure and fluid saturation changes from marine multicomponent time-lapse seismic data: Geophysics, 68, 1592–1599; doi: 10.1190/1.1620633.
  * Lang and Grana, (2019) Lang, X., and D. Grana, 2019, Rock physics modelling and inversion for saturation-pressure changes in time-lapse seismic studies: Geophysical Prospecting, 67, 1912–1928.
  * Lavoué et al., (2014) Lavoué, F., R. Brossier, L. Métivier, S. Garambois, and J. Virieux, 2014, Two-dimensional permittivity and conductivity imaging by full waveform inversion of multioffset GPR data: A frequency-domain quasi-Newton approach: Geophysical Journal International, 197, 248–268; doi: 10.1093/gji/ggt528.
  * Lee, (2005) Lee, M. W., 2005, Proposed moduli of dry rock and their application to predicting elastic velocities of sandstones: US Department of the Interior, US Geological Survey.
  * Lescoffit et al., (2016) Lescoffit, S. P., M. Houbiers, C. Henstock, E. Hicks, K.-M. Nilsen, H. Hoeber, A. Ratcliffe, and V. Vinje, 2016, Time-lapse full-waveform inversion applied to permanent reservoir monitoring data from Grane, a Norwegian North Sea field: SEG Technical Program Expanded Abstracts 2016, Society of Exploration Geophysicists, 1304–1308.
  * Lumley, (2001) Lumley, D. E., 2001, Time-lapse seismic reservoir monitoring: Geophysics, 66, 50–53; doi: 10.1190/1.1444921.
  * Ma et al., (2016) Ma, Y., M. Maharramov, R. Clapp, and B. Biondi, 2016, Multiparameter full-waveform inversion in the isotropic acoustic media with application to time-lapse seismic inverse problem: SEG Technical Program Expanded Abstracts 2016, OnePetro, 5557–5561.
  * Maharramov et al., (2016) Maharramov, M., B. L. Biondi, and M. A. Meadows, 2016, Time-lapse inverse theory with applications: Geophysics, 81, R485–R501; doi: 10.1190/geo2016-0131.1.
  * Mallick, (2007) Mallick, S., 2007, Amplitude-variation-with-offset, elastic-impedence, and wave-equation synthetics—a modeling study: Geophysics, 72, C1–C7.
  * (28) Mardan, A., B. Giroux, and G. Fabien-Ouellet, 2022a, Time-lapse full-waveform inversion for monitoring the fluid saturation: Second EAGE Conference on Seismic Inversion, European Association of Geoscientists & Engineers, 1–5.
  * (29) ——–, 2022b, Time-lapse seismic full waveform inversion using improved cascaded method: Second EAGE Conference on Seismic Inversion, European Association of Geoscientists & Engineers, 1–5.
  * (30) ——–, 2023a, PyFWI: A Python package for full-waveform inversion and reservoir monitoring; doi: 10.2139/ssrn.4330227.
  * (31) ——–, 2023b, Weighted-average time-lapse seismic full-waveform inversion: Geophysics, 88, R25–R38; doi: 10.1190/geo2022-0090.1.
  * (32) Mardan, A., B. Giroux, G. Fabien-Ouellet, and M. Saberi, 2022c, Direct monitoring of fluid saturation using time-lapse full-waveform inversion: Presented at the 2nd International Meeting for Applied Geoscience & Energy Expanded Abstracts, Society of Exploration Geophysicists.
  * Mavko et al., (2020) Mavko, G., T. Mukerji, and J. Dvorkin, 2020, The rock physics handbook: Cambridge university press.
  * Métivier et al., (2013) Métivier, L., R. Brossier, J. Virieux, and S. Operto, 2013, Full waveform inversion and the truncated Newton method: SIAM Journal on Scientific Computing, 35, B401–B437; doi: 10.1137/120877854.
  * Naeini et al., (2017) Naeini, E. Z., N. Kamath, I. Tsvankin, and T. Alkhalifah, 2017, Full waveform inversion for reservoir characterization-a synthetic study: 79th EAGE Conference and Exhibition 2017, European Association of Geoscientists & Engineers, 1–5.
  * Nocedal and Wright, (2006) Nocedal, J., and S. J. Wright, 2006, Numerical optimization: Springer, New York, NY.
  * Operto et al., (2013) Operto, S., Y. Gholami, V. Prieux, A. Ribodetti, R. Brossier, L. Metivier, and J. Virieux, 2013, A guided tour of multiparameter full-waveform inversion with multicomponent data: From theory to practice: The Leading Edge, 32, 1040–1054; doi: 10.1190/tle32091040.1.
  * Pan et al., (2019) Pan, W., K. A. Innanen, Y. Geng, and J. Li, 2019, Interparameter trade-off quantification for isotropic-elastic full-waveform inversion with various model parameterizations: Geophysics, 84, R185–R206; doi: 10.1190/GEO2017-0832.1.
  * Plessix et al., (2010) Plessix, R., S. Michelet, H. Rynja, H. Kuehl, C. Perkins, J. de Maag, and P. Hatchell, 2010, Some 3D applications of full waveform inversion: Presented at the 72nd EAGE Conference and Exhibition-Workshops and Fieldtrips, European Association of Geoscientists & Engineers.
  * Pride, (2005) Pride, S., 2005, in Hydrogeophysics: Water science and technology library: Springer, 253–284.
  * Queißer and Singh, (2013) Queißer, M., and S. C. Singh, 2013, Full waveform inversion in the time lapse mode applied to CO2 storage at Sleipner: Geophysical Prospecting, 61, 537–555; doi: 10.1111/j.1365-2478.2012.01072.x.
  * Raknes et al., (2013) Raknes, E. B., W. Weibull, and B. Arntsen, 2013, Time-lapse full waveform inversion: Synthetic and real data examples: SEG Technical Program Expanded Abstracts 2013, Society of Exploration Geophysicists, 944–948.
  * Romdhane et al., (2022) Romdhane, A., B. Dupuy, E. Querendez, and P. Eliasson, 2022, Toward quantitative CO2 monitoring at Sleipner, Norway: Geophysical Monitoring for Geologic Carbon Storage, 383–402.
  * Tarantola, (1986) Tarantola, A., 1986, A strategy for nonlinear elastic inversion of seismic reflection data: Geophysics, 51, 1893–1903; doi: 10.1190/1.1442046.
  * Vedanti and Sen, (2009) Vedanti, N., and M. K. Sen, 2009, Seismic inversion tracks in situ combustion: A case study from Balol oil field, India: Geophysics, 74, B103–B112; doi: 10.1190/1.3129262.
  * Virieux et al., (2017) Virieux, J., A. Asnaashari, R. Brossier, L. Métivier, A. Ribodetti, and W. Zhou, 2017, in 6\. An introduction to full waveform inversion: Society of Exploration Geophysicists, R1–1–R1–40.
  * Virieux and Operto, (2009) Virieux, J., and S. Operto, 2009, An overview of full-waveform inversion in exploration geophysics: Geophysics, 74, WCC1–WCC26; doi: 10.1190/1.3238367.
  * Voigt, (1889) Voigt, W., 1889, Ueber die beziehung zwischen den beiden elasticitätsconstanten isotroper körper: Annalen der physik, 274, 573–587.
  * Watanabe et al., (2005) Watanabe, T., S. Shimizu, E. Asakawa, and T. Matsuoka, 2005, in Differential waveform tomography for time-lapse crosswell seismic data with application to gas hydrate production monitoring: Society of Exploration Geophysicists, 2323–2326.
  * Yan et al., (2019) Yan, H., B. Dupuy, A. Romdhane, and B. Arntsen, 2019, CO2 saturation estimates at Sleipner (North Sea) from seismic tomography and rock physics inversion: Geophysical Prospecting, 67, 1055–1071; doi: 10.1111/1365-2478.12693.
  * Yang et al., (2018) Yang, P., R. Brossier, L. Métivier, J. Virieux, and W. Zhou, 2018, A time-domain preconditioned truncated Newton approach to visco-acoustic multiparameter full waveform inversion: SIAM Journal on Scientific Computing, 40, B1101–B1130; doi: 10.1137/17M1126126.
  * Zhou and Lumley, (2021) Zhou, W., and D. Lumley, 2021, Central-difference time-lapse 4D seismic full-waveform inversion: Geophysics, 86, R161–R172; doi: 10.1190/geo2019-0834.1.

##  Appendix A Regularization in seismic FWI

By taking advantage of Tikhonov regularization, smoothness of the estimated
model can be achieved Asnaashari et al., (2013). This regularization can be
presented as

$\chi_{1}=\lambda_{1}\left(\mathbf{m}^{T}\mathbf{B}_{x}^{T}\mathbf{B}_{x}\mathbf{m}+\mathbf{m}^{T}\mathbf{B}_{z}^{T}\mathbf{B}_{z}\mathbf{m}\right),$
(18)

where $\mathbf{B}$ and $\lambda$ denote respectively the first-order spatial
derivative operator and regularization hyperparameter. The subscripts $x$ and
$z$ specify the direction.

Total variation regularization (TV) has also shown promising results in FWI to
decrease artifacts and preserve the edges and the discontinuities Anagaw and
Sacchi, (2012); Esser et al., (2018). This regularization can be achieved by
using

$\chi_{2}=\lambda_{2}\|\sqrt{(\mathbf{B}_{x}\mathbf{m})^{2}+(\mathbf{B}_{z}\mathbf{m})^{2}}\|_{1}.$
(19)

Also, Asnaashari et al., (2013) state that where nonseismic information
exists, it should be used as prior information for FWI inversion. Hence, they
proposed regularizing the cost function using the prior information
($\mathbf{m}_{p}$) as

$\chi_{3}=\lambda_{3}\|\mathbf{W}_{m}(\mathbf{m}-\mathbf{m}_{p})\|_{2}^{2},$
(20)

where $\mathbf{W}_{m}$ is a weighting operator in the model space.

Finally, in the case of multiparameter FWI, if information about the
relationship between different parameters exists, then regularization can be
considered to force the optimization algorithm to follow Hu and Innanen,
(2021),

$\chi_{4}=\lambda_{4}\|\mathbf{m}_{1}-f(\mathbf{m}_{2})\|_{2}^{2},$ (21)

where $f$ is a function to estimates $\mathbf{m}_{1}$ from $\mathbf{m}_{2}$.

##  Appendix B The gradient of cost function in PCS parameterization

In equation 13, we showed how to calculate the elastic properties of the
porous medium using $K$, $G$, and $\rho$. To perform FWI with PCS
parameterization, it is required to estimate the gradient of these parameters
in terms of PCS parameterization. Gradient of $K$, $G$, and $\rho$ with
respect to $\phi$ can be obtained as

$\displaystyle\begin{split}\frac{\partial
K}{\partial\phi}&=\frac{K_{D}+\phi\frac{\partial
K_{D}}{\partial\phi}-\frac{K_{f}}{K_{s}}\left(K_{D}+(1+\phi)\frac{\partial
K_{D}}{\partial\phi}\right)}{\phi(1+\Delta)}-\\\ &\frac{\left[\phi
K_{D}+\left(1-\frac{(1+\phi)K_{D}}{K_{s}}\right)K_{f}\right]\left(1+\Delta+\phi\frac{\partial\Delta}{\partial\phi}\right)}{\phi^{2}(1+\Delta)^{2}},\\\
\frac{\partial G}{\partial\phi}&=\frac{\partial G_{D}}{\partial\phi},\\\
\frac{\partial\rho}{\partial\phi}&=\rho_{f}-\rho_{s},\end{split}$ (22)

where

$\displaystyle\begin{split}\frac{\partial
K_{D}}{\partial\phi}&=-K_{s}\frac{1+cs}{(1+cs\phi)^{2}},\\\
\frac{\partial\Delta}{\partial\phi}&=-\frac{K_{f}}{K_{s}\phi^{2}}\left(1-\frac{1}{1+cs\phi}\right)+\frac{K_{f}(1-\phi)}{K_{s}\phi}\frac{cs}{(1+cs\phi)^{2}},\\\
\frac{\partial
G_{D}}{\partial\phi}&=-G_{s}\frac{1+\frac{3}{2}cs}{(1+\frac{3}{2}cs\phi)^{2}}.\end{split}$
(23)

With respect to $C$, we have

$\displaystyle\begin{split}\frac{\partial K}{\partial
C}&=\frac{\phi\frac{\partial K_{D}}{\partial
C}+\frac{K_{f}(1+\phi)}{K_{s}^{2}}\left(\frac{\partial K_{s}}{\partial
C}K_{D}-\frac{\partial K_{D}}{\partial C}K_{s}\right)}{\phi(1+\Delta)}-\\\
&\frac{\left[\phi
K_{D}+\left(1-\frac{(1+\phi)K_{D}}{K_{s}}\right)K_{f}\right]\phi\frac{\partial\Delta}{\partial
C}}{\phi^{2}(1+\Delta)^{2}},\\\ \frac{\partial G}{\partial C}&=\frac{\partial
G_{D}}{\partial C},\\\ \frac{\partial\rho}{\partial
C}&=(1-\phi)\frac{\partial\rho_{s}}{\partial C},\end{split}$ (24)

where

$\displaystyle\begin{split}\frac{\partial\Delta}{\partial
C}&=\frac{(\phi-1)K_{f}}{\phi
K_{s}^{2}}\left(1-\frac{1}{1+cs\phi}\right)\frac{\partial K_{s}}{\partial
C},\\\ \frac{\partial K_{s}}{\partial C}&=K_{c}-K_{q},\qquad&&\frac{\partial
G_{s}}{\partial C}=G_{c}-G_{q},\\\ \frac{\partial K_{D}}{\partial
C}&=\frac{1-\phi}{1+cs\phi}\frac{\partial K_{s}}{\partial
C},\qquad&&\frac{\partial G_{D}}{\partial
C}=\frac{1-\phi}{1+\frac{3}{2}cs\phi}\frac{\partial G_{s}}{\partial C},\\\
\frac{\partial\rho_{s}}{\partial C}&=(\rho_{c}-\rho_{q}).\end{split}$ (25)

Finally, with respect to $S_{w}$, we have

$\displaystyle\begin{split}\frac{\partial K}{\partial
S_{w}}&=\frac{\left(1-\frac{(1+\phi)K_{D}}{K_{s}}\right)\frac{\partial
K_{f}}{\partial S_{w}}}{\phi(1+\Delta)}-\\\ &\frac{\left[\phi
K_{D}+\left(1-\frac{(1+\phi)K_{D}}{K_{s}}\right)K_{f}\right]\phi\frac{\partial\Delta}{\partial
S_{w}}}{\phi^{2}(1+\Delta)^{2}},\\\ \frac{\partial G}{\partial S_{w}}&=0,\\\
\frac{\partial\rho}{\partial S_{w}}&=\phi(\rho_{w}-\rho_{h}),\end{split}$ (26)

where

$\displaystyle\begin{split}\frac{\partial K_{f}}{\partial
S_{w}}&=K_{w}-K_{h},\\\ \frac{\partial\Delta}{\partial
S_{w}}&=\frac{\Delta}{K_{f}}\frac{\partial K_{f}}{\partial S_{w}}.\end{split}$
(27)

Gathering the previous terms, the gradient of the cost function with respect
to model parameters in PCS parameterization can be calculated as

$\left[\begin{array}[]{c}\frac{\partial\chi}{\partial\phi}\\\
\frac{\partial\chi}{\partial C}\\\ \frac{\partial\chi}{\partial
S_{w}}\end{array}\right]=\left[\begin{array}[]{c c c}\frac{\partial
V_{P}}{\partial\phi}&\frac{\partial
V_{S}}{\partial\phi}&\frac{\partial\rho}{\partial\phi}\\\ \frac{\partial
V_{P}}{\partial C}&\frac{\partial V_{S}}{\partial
C}&\frac{\partial\rho}{\partial C}\\\ \frac{\partial V_{P}}{\partial
S_{w}}&\frac{\partial V_{S}}{\partial S_{w}}&\frac{\partial\rho}{\partial
S_{w}}\end{array}\right]\left[\begin{array}[]{c}\frac{\partial\chi}{\partial
V_{P}}\\\ \frac{\partial\chi}{\partial V_{S}}\\\
\frac{\partial\chi}{\partial\rho}\end{array}\right].$ (28)