# Bayesian Time-Lapse Full Waveform Inversion
using Hamiltonian Monte Carlo

P. D. S. de Lima  0000-0002-7353-536X Departamento de Física Teórica e
Experimental,
Universidade Federal do Rio Grande do Norte, 59078-970 Natal-RN, Brazil
School of Physics, Trinity College Dublin, Dublin 2, Ireland  M. S. Ferreira
0000-0002-0856-9811 School of Physics, Trinity College Dublin, Dublin 2,
Ireland  Centre for Research on Adaptive Nanostructures and Nanodevices
(CRANN) & Advanced Materials and Bioengineering Research (AMBER) Centre,
Trinity College Dublin, Dublin 2, Ireland  G. Corso  0000-0003-1748-4040
Departamento de Biofísica e Farmacologia,
Universidade Federal do Rio Grande do Norte, 59078-970 Natal-RN, Brazil  J. M.
de Araújo  0000-0001-8462-4280 Departamento de Física Teórica e Experimental,
Universidade Federal do Rio Grande do Norte, 59078-970 Natal-RN, Brazil

###### Abstract

Time-lapse images carry out important information about dynamic changes in
Earth’s interior which can be inferred using different Full Waveform Inversion
(FWI) schemes. The estimation process is performed by manipulating more than
one seismic dataset, associated with the baseline and monitors surveys. The
time-lapse variations can be so minute and localised that quantifying the
uncertainties becomes fundamental to assessing the reliability of the results.
The Bayesian formulation of the FWI problem naturally provides confidence
levels in the solution, but evaluating the uncertainty of time-lapse seismic
inversion remains a challenge due to the ill-posedness and high dimensionality
of the problem. The Hamiltonian Monte Carlo (HMC) can be used to effectively
sample over high dimensional distributions with affordable computational
efforts. In this context, we propose a probabilistic Bayesian sequential
approach for time-lapse FWI using the HMC method. Our approach relies on the
integration of the baseline survey information as prior knowledge in the
monitor estimation. We compare the proposed methodology with a parallel scheme
in perfect and perturbed acquisition geometry scenarios. We also investigate
the correlation effect between baseline and monitor samples in the propagated
uncertainties. The results show that our strategy provides accurate times-
lapse estimates with errors of similar magnitude to the parallel methodology.

##  I Introduction

Extracting information from noisy and limited data is paramount in subsurface
characterization procedures. In particular, the availability of time-lapse
data, that is, multiple seismic data acquired at different calendar times,
offers not only the opportunity to determine static structures but also
dynamic changes in geological models [1]. In this context, the first seismic
survey is called baseline while the remaining ones are referred to as
monitors. Constructing time-lapse images is crucial to control the oil
production behaviour [2] in reservoirs or to monitor the CO2 storage [3, 4] in
the injection process.

Regardless of its application, time-lapse models can be constructed using the
Full Waveform Inversion (FWI) [5, 6]. This high-resolution technique takes
into account both the phase and amplitude of the wavefield in the model
parameter estimation. However, additional challenges are present in the time-
lapse inversion due to the small and spatially localized nature of variations
combined with the non-repeatability of the acquisition geometry and climate
factors [7].

Several deterministic time-lapse FWI methodologies have been proposed in the
literature. The parallel strategy [8, 9] independently inverts both baseline
and monitor data from the same initial model. The double-difference [10, 11,
12, 13] inverts directly from data difference while the jointly-inversion [14,
15, 16, 17] obtains both models simultaneously. Another approach is the
sequential difference [18, 19, 20, 21], which uses the inverted baseline model
as the initial model in the monitor inversion. Generalizations of the
sequential scheme have been recently designed to overcome the non-
repeatability problem [22, 23]. However, the uncertainty regarding the time-
lapse estimates cannot be properly inferred using these deterministic
approaches [24]. In fact, the ill-posed nature of the FWI problem produces
solutions with nontrivial propagated errors, such as noise in data, model
parametrization and simplification of the physical theory [25]. Therefore, the
uncertainty quantification process is fundamental for evaluating result
reliability and facilitating the decision-making process.

Bayesian approaches address the challenge of uncertainty quantification by
transforming the seismic inversion into a statistical inference problem, where
the results become probability density functions which are sampled from the
model space. Markov Chain Monte Carlo (MCMC) has been used over decades in
geophysical inversion problems [24, 26, 27] for this task. However, MCMC
methods are inefficient in estimating probability distributions in high-
dimensional spaces, due to the well-known curse of dimensionality effect. This
dimensionality issue is substantially suppressed when the Hamiltonian Monte
Carlo (HMC) method is used [28]. The HMC is a sampling technique that explores
the geometry of target distribution to avoid the random walk behaviour of
standard MCMC methods [29]. Originally formulated for quantum chromodynamics
[30], the HMC has been popularized in the seismic inversion community in the
last decade [31, 32, 33, 34, 35, 36, 37]. Alternatively to MCMC-based methods,
the ensemble-based [38, 39, 40, 41] and variational [42, 43, 44, 45, 46, 47]
approaches have also been employed to address similar problems.

Advances in Bayesian time-lapse estimation have been performed in [48, 49],
where the uncertainties are estimated directly from data changes using data
compression techniques [50]. However, this approach is valid only when the
acquisition geometry is repeatable, a strong assumption for real-data
applications. In Ref. [47], Zhang and Curtis performed a comparison between
the separate and jointly Bayesian time-lapse inversion using a stochastic
Stein Variational Gradient Descent (sSVGD) method [46] by imposing a prior
distribution directly in time-lapse estimates. In this paper, we explore a
different avenue by considering the posterior information regarding the
baseline survey to design a prior distribution to be used in the monitor
estimation. Results for this sequential strategy in a simplified scenario
where the baseline model is known was reported in [51]. Here, we follow a
general approach and demonstrate the benefits of employing informative prior
information in time-lapse studies by comparing the sequential and parallel
approaches in terms of estimates, propagated uncertainties and sample
correlation. Moreover, we also consider in our investigations the
nonrepeatability effects associated with different source locations in the
surveys.

In the following section, we describe the Bayesian formulation of the FWI
problem and extensions to the time-lapse case considering the parallel and
sequential strategies. We also present the HMC as a sampling method and
details involving the tuning strategy to the mass matrix in this section. In
section 4, we compare the statistical estimates obtained using both parallel
and sequential strategies. We also compare the Bayesian inversion results with
deterministic ones. We conclude by discussing the limitations of sequential
strategy and possible improvements.

##  II Methods

###  II.1 Bayesian Full Waveform Inversion

The solution of a FWI problem from a Bayesian perspective is not a single
model but a probability density function of models [25], called the posterior
distribution $\rho(\mathbf{m}|\mathbf{d})$. The posterior distribution
combines the prior knowledge of the model $\mathbf{m}$ with the information
acquired through data measurements $\mathbf{d}^{\text{obs}}$ via Bayes’s
theorem [52]:

$\rho(\mathbf{m}|\mathbf{d})=\frac{\rho(\mathbf{m})\rho(\mathbf{d}|\mathbf{m})}{\rho(\mathbf{d})}\,,$
(1)

where $\rho(\mathbf{m})$ is the prior probability distribution which
encapsulates the previous model knowledge and the likelihood function
$\rho(\mathbf{d}|\mathbf{m})$ that encompasses the conditional probability of
an observation being explained by a certain model. The denominator in equation
(1) is called the Bayesian evidence (or partition function):

$\rho(\mathbf{d})=\int_{\mathbb{M}}\rho(\mathbf{d}|\mathbf{m})\rho(\mathbf{m})\,d\mathbf{m}\,,$
(2)

which normalizes the posterior distribution.

In this regard, the set of plausible model parameters $\mathbf{m}$ spans a
model space $\mathbb{M}$ and the correspondent mapping to the data space
$\mathbb{D}$ is done by a nonlinear forward operator
$\mathbf{d}^{\text{mod}}(\mathbf{m})=\mathbf{G}(\mathbf{m})$. The error
between observations and data predictions from a proposal model is called
residuals $\Delta\mathbf{d}=\mathbf{d}^{\text{mod}}-\mathbf{d}^{\text{obs}}$.
In light of the central limit theorem, we assume that the residual
distribution is Gaussian. Therefore, the likelihood function can be written
as:

$\rho(\mathbf{d}|\mathbf{m})=\exp{\left[-\frac{1}{2}\Delta\mathbf{d}^{\text{T}}(\Sigma^{\text{noise}})^{-1}\Delta\mathbf{d}\right]}\,,$
(3)

where $\Sigma^{\text{noise}}$ is the noise covariance matrix. In our synthetic
studies, we assume that the data residual error variance is known, even though
the residual error distribution in more realistic cases can be estimated [53]
or may take different forms [54, 55, 56, 57]. We notice that in the
uncorrelated data case, the matrix $\Sigma^{\text{noise}}$ has only diagonal
elements. The statistical analysis of the posterior distribution is typically
evaluated by computing the posterior mean

$\mathbf{m}^{\text{mean}}=\int_{\mathbb{M}}\mathbf{m}\>\rho(\mathbf{m}|\mathbf{d})\,d\mathbf{m}\,,$
(4)

and the posterior covariance

$\Sigma_{\textbf{m}}=\int_{\mathbb{M}}(\mathbf{m}-\mathbf{m}^{\text{mean}})^{\text{T}}(\mathbf{m}-\mathbf{m}^{\text{mean}})\rho(\mathbf{m}|\mathbf{d})\,d\mathbf{m}$
(5)

of the model parameter $\mathbf{m}$. However, in practical FWI problems we
have only $n$ representative model samples
$(\mathbf{m}^{(1)},\mathbf{m}^{(2)}\ldots,\mathbf{m}^{(n)})$ as estimation of
the posterior distribution $\rho(\mathbf{m}|\mathbf{d})$.

###  II.2 Bayesian Time-Lapse Full Waveform Inversion

As far as time-lapse problems are concerned, we need to deal with two
different data spaces, each one associated with different seismic
acquisitions, namely the baseline $\mathbf{d}^{\text{obs}}_{\text{B}}$ and the
monitor $\mathbf{d}^{\text{obs}}_{\text{M}}$ data. Deterministically, this
information can be used to construct the baseline $\mathbf{m}_{\text{B}}$ and
monitor $\mathbf{m}_{\text{M}}$ models, although we are primarily interested
in the time-lapse
$\mathbf{m}_{\text{TL}}=\mathbf{m}_{\text{M}}-\mathbf{m}_{\text{B}}$ model
change. Accordingly, in the probabilistic realm we need to handle two
probability distributions to infer time-lapse changes, the baseline
$\rho(\mathbf{m}_{\text{B}}|\mathbf{d}_{\text{B}})\equiv\rho_{\text{B}}(\mathbf{m}|\mathbf{d})$
and the monitor
$\rho(\mathbf{m}_{\text{M}}|\mathbf{d}_{\text{M}})\equiv\rho_{\text{M}}(\mathbf{m}|\mathbf{d})$
posterior distributions.

The straightforward Bayesian inference approach is to estimate the uncertainty
on the time-lapse changes $\mathbf{m}_{\text{TL}}$ directly from data
difference
$\mathbf{d}_{\text{TL}}=\mathbf{d}_{\text{M}}-\mathbf{d}_{\text{B}}$ [58].
Although this strategy avoids separately sampling each distribution, thus
reducing the computational cost, it cannot be applied in realistic scenarios
with nonrepeatable acquisitions.

Looking at the two distributions separately, one possible strategy is to
estimate the baseline $\rho_{\text{B}}(\mathbf{m}|\mathbf{d})$ and monitor
$\rho_{\text{M}}(\mathbf{m}|\mathbf{d})$ using the same prior information,
which means
$\rho_{\text{B}}(\mathbf{m})=\rho_{\text{M}}(\mathbf{m})=\rho(\mathbf{m})$.
Alternatively, one can follow a joint Bayesian inversion [47], where baseline
and time-lapse probabilities are obtained simultaneously. However, as the
time-lapse change is small and spatially localized and the time interval
between the acquisition is long enough to construct well-resolved velocity
models, it is natural to use all information learned from the baseline
inversion as priors for the monitor estimation. Therefore, we follow a
probabilistic version of the sequential strategy, where the prior probability
used in the monitor inversion is derived from the baseline posterior
distribution.

Following this sequential scheme, we first draw $n_{\text{B}}$ samples
$(\mathbf{m}^{(1)}_{\text{B}},\mathbf{m}^{(2)}_{\text{B}},\ldots,\mathbf{m}^{(n_{\text{B}})}_{\text{B}})$
of the baseline posterior $\rho_{\text{B}}(\mathbf{m}|\mathbf{d})$ which are
then used to design a prior distribution to the monitor, that is, we set
$\rho_{\text{M}}(\mathbf{m})=\rho_{\text{B}}(\mathbf{m}|\mathbf{d})$. Using
this prior, we collect $n_{\text{M}}$ samples
$(\mathbf{m}^{(1)}_{\text{M}},\mathbf{m}^{(2)}_{\text{M}},\ldots,\mathbf{m}^{(n_{\text{M}})}_{\text{M}})$
of the monitor posterior distribution
$\rho_{\text{M}}(\mathbf{m}|\mathbf{d})$. These samples are then used to
access the time-lapse statistical information. In particular, the correlation
between baseline and monitor samples can be assessed by computing the
covariance matrix:

$\Sigma_{\text{MB}}=\int_{\mathbb{M}}(\mathbf{m}_{\text{M}}-\mathbf{m}_{\text{M}}^{\text{mean}})^{\text{T}}(\mathbf{m}_{\text{B}}-\mathbf{m}_{\text{B}}^{\text{mean}})\rho_{\mathbf{MB}}(\mathbf{m}|\mathbf{d})\,d\mathbf{m}\,,$
(6)

where
$\rho_{\mathbf{MB}}(\mathbf{m}|\mathbf{d})\equiv\rho(\mathbf{m}_{\text{M}},\mathbf{m}_{\text{B}}|\mathbf{d}_{\text{M}},\mathbf{d}_{\text{B}})$
is the associated joint probability distribution. In the remainder of the
paper, our focus will be on the comparison of the parallel and sequential
Bayesian inversions using the HMC as the sampling method.

###  II.3 Hamiltonian Monte Carlo

The Hamiltonian Monte Carlo method incorporates the geometry of the target
distribution to efficiently sample over high-dimensional model spaces [28]. In
the HMC context, the model parameters $\mathbf{m}$ are interpreted as a set of
particles following trajectories of a classical mechanical system. To complete
the set of canonical variables, one defines the momenta
$\mathbf{p}\in\mathbb{P}$ which is randomly sampled according to a Gaussian
distribution with zero mean and a covariance given by the matrix mass
$\mathbf{M}$. In this way, the Hamiltonian system is subjected to a fake
potential energy that mimics the negative value of the log-posterior
distribution (Eq. 1). Accordingly, the model space $\mathbb{M}$ is extended to
a (fake) phase space $\mathbb{Z}=\mathbb{M}\times\mathbb{P}$, such that the
posterior distribution is sampled over the canonical distribution

$\rho(\mathbf{m},\mathbf{p})=\rho(\mathbf{p})\rho(\mathbf{m}|\mathbf{d})\propto\exp{(-H(\mathbf{m},\mathbf{p}))}\,,$
(7)

with a Hamiltonian $H(\mathbf{m},\mathbf{p})$ given by

$H(\mathbf{m},\mathbf{p})=\frac{1}{2}\mathbf{p}^{\text{T}}\mathbf{M}^{-1}\mathbf{p}-\log{\rho(\mathbf{m}|\mathbf{d})}\,.$
(8)

To sample the distribution (Eq. 7), one starts with a model $\mathbf{m}_{0}$
drawn from the prior distribution $\rho(\mathbf{m})$, which is then combined
with initial random momenta $\mathbf{p}_{0}$ to generate the initial state of
the system. Subsequently, the system evolves following the dynamics of
Hamilton’s equations [59]:

$\frac{d\mathbf{m}}{d\tau}=\frac{\mathbf{p}}{\mu}\,,\quad\frac{d\mathbf{p}}{d\tau}=\nabla\log{\rho(\mathbf{m}|\mathbf{d})}\,,$
(9)

during a (artificial) time period $\tau$. The evolved state
$(\mathbf{m}_{\tau},\mathbf{p}_{\tau})$ is accepted with a probability given
by the Metropolis-Hasting criteria [60]

$\mathrm{min}\left[1,\exp{(H(\mathbf{m}_{0},\mathbf{p}_{0})-H\left(\mathbf{m}_{\tau},\mathbf{p}_{\tau})\right)}\right]\,.$
(10)

Since we are interested in the estimated posterior (Eq. 1), the canonical
distribution is marginalized by discarding the momenta samples. Each sample
$\mathbf{m}_{\tau}$ is stored and becomes the new initial position of the
system, and it is again evolved according to the dynamical equations (Eq. 9)
using new random momentum. This procedure is a single HMC step, which is then
repeated $n_{\text{HMC}}$ times, producing a Markov chain with $n$ samples of
the canonical distribution (Eq. 7). We notice that only a fraction of the
samples $(n\leq n_{\text{HMC}})$ is accepted, such that the acceptance rate
$n/n_{\text{HMC}}$ is typically used as a control parameter to tune the HMC
efficiency [61]. Therefore, it is recommended that the HMC parameters be
adjusted to maintain high values ($n/n_{\text{HMC}}>0.6$) of the acceptance
rate.

####  II.3.1 Mass Matrix Tuning

The mass matrix $\mathbb{M}$ has a strong impact on speeding up the phase
space exploration and needs to be carefully adjusted [62]. Improper choices
can lead to high acceptance rates without an effective exploration of the
phase space. On the other hand, assigning distinct mass values to different
model parameters can be used to prioritize certain regions of the space [33,
34], working similarly as regularizers in deterministic inversion [63, 20]. In
particular, the optimal mass matrix choice for linear problems corresponds to
the inverse of posterior covariance. However, the latter matrix is too large
to be fully computed in practical problems and is unknown before the
inversion. Thus, convenient choices are necessary to perform Bayesian
inversion using HMC in high-dimensional problems.

We follow the strategy proposed in [37], where each particle initially has a
pseudo-mass $\mu$. After $k$ HMC steps, the mass matrix is updated to
$\mathbf{M}^{\prime}=\boldsymbol{\Gamma}^{-1}\mathbf{M}$, where
$\boldsymbol{\Gamma}=\text{diag}(\gamma_{1},\gamma_{2},\ldots,\gamma_{\text{dim}(\mathbb{M})})$
and

$\gamma_{i}=\gamma_{\text{min}}+(\gamma_{\text{max}}-\gamma_{\text{min}})\left(\frac{z-z_{\text{water}}}{z_{\text{max}}-z_{\text{water}}}\right)\,,$
(11)

where $z_{\text{water}}$ and $z_{\text{max}}$ are the water layer depth and
maximum model depth, respectively. The expression (Eq. 11) adjusts the mass of
the particles according to its depth $z$ in the physical domain. This approach
allows us to improve the HMC convergence using fewer forward calculations in
comparison with the standard identity matrix choice. In addition, it also
avoids additional Hessian approximation calculations [37]. In our study, we
fix the initial mass to $\mu=10$ and the minimum and maximum gamma values to
$\gamma_{\text{min}}=1.0$ and $\gamma_{\text{max}}=2.0$.

##  III Results

###  III.1 Experimental Setup

Figure 1: (a) The true baseline model with the red (magenta) star represents
one of the ten sources used in the baseline (monitor perturbed) acquisition
and the orange square represents some of the 200 receivers employed in our
experiments. (b) The true time-lapse model. (c) The mean prior model used in
the baseline estimation for both inversion strategies.

To compare both probabilistic time-lapse inversion approaches, our analysis is
focused on the Marmousi model, which has been extensively used as a laboratory
for time-lapse inversion studies [20, 64, 23, 46]. The domain is a target zone
(Figure 1a) of the original model with 3.0 km of maximum offset and 1.5 km of
depth, with a spatial resolution of 20 m in both directions. Therefore, the
sampling model space has over $\sim 10^{4}$ dimensions for each survey.

We proceed as [20] and consider a time-lapse velocity model (Figure 1b) which
has $\sim 0.1$ km/s of velocity variation inside two gas reservoirs. These
$\sim 5\%$ realistic changes represent the partial gas replacement by water.
Initially, we simplify our analysis, considering a repeatable acquisition
geometry (red star in Figure 1a) between baseline and monitor surveys. After
that, we consider a nonrepeatable scenario where the sources in the monitor
survey (magenta star in Figure 1a) are horizontally shifted by 100 m.

The synthetic baseline and monitor data were generated using 10 isotropic
explosive sources (Ricker wavelet with 10 Hz of central frequency) at 40 m
depth and recorded during 3.5 s by 200 receivers spaced every 20 m, deployed
at 120 m depth. In accordance with Ref. [45], we also contaminate our seismic
data with 0.1 standard deviation uncorrelated Gaussian noise. The
corresponding forward modelling was performed by simulating a constant density
acoustic wave propagation following a finite difference scheme with
derivatives of eight-order in space and second-order in time. In addition, the
gradient of the log-posterior (Eq. 1) used in Hamilton’s equations (Eq. 9) was
computed using the adjoint state method [65, 6].

###  III.2 Baseline Estimation and Monitor Prior Design

Figure 2: The mean (a) and the standard deviation (b) estimates of the
baseline model using the HMC method.

The baseline posterior distribution $\rho_{\text{B}}(\mathbf{m}|\mathbf{d})$
was estimated from a prior probability $\rho_{\text{B}}(\mathbf{m})$ with mean
model illustrated in Figure 1c. This uninformative prior consists of
independent and uniform distributions for each model parameter that assumes
only the tomographic knowledge of the geological structure and extreme
velocity values. About $800$ baseline samples were generated using 1200 HMC
steps without discard, that is, no burn-in period was considered.

The Bayesian inversion results for the baseline survey are illustrated in
Figure 2, where we focus on analysing the first statistical moment and the
dispersion associated with respective posterior distribution. The main
structures of the geological model are presented in the mean model (Figure
2a), particularly in the shallow part ($z<1$ km), while only large-scale
features were estimated at the deeper region ($z>1$ km). As demonstrated by
[34] and [37], the noise data has a strong impact in assessing correct
velocities mainly as the depth increases. Moreover, due to the limitations of
the acquisition setup, it is difficult to detect dynamic structures around the
bottom edge of the model. In fact, these features of the estimated baseline
model are mapped as uncertainties in Figure 2b, where the standard deviation
increases with depth and is higher in regions between discontinuities of the
model, as in the region of reservoirs that we have particular interest.
Similar uncertainty results using different sampling methods can be found in
[66, 44, 45, 37].

Figure 3: Comparison between deterministic and probabilistic inversions. The
maximum a posteriori (MAP) baseline (a) and the model obtained using the
deterministic inversion (b). Inverted time-lapse models using the
deterministic version of the parallel (left) and sequential strategies (right)
with perfect (c, d) and perturbed (e, f) acquisition geometry.

As previously mentioned, deterministic methods fail to offer a precise
framework for quantifying uncertainties of estimated velocity models.
Nevertheless, we utilize them to compare the inversion results achieved
through deterministic time-lapse workflows (parallel and sequential) with
those obtained by employing their Bayesian probabilistic versions. For this
purpose, we use the L-BFGS method [67] to minimize the misfit function
(argument of the likelihood function in Eq. 3) in both cases, where the same
experimental setup described at the beginning of this section was used. First,
we compare the sample model that maximizes the posterior distribution
$\rho_{\text{B}}(\mathbf{m}|\mathbf{d})$, that is, the maximum a posteriori
(MAP) baseline model (Figure 3a), with those obtained using the L-BFGS method
(Figure 3b) for the deterministic inversion. The two models (Figure 3a and b)
are quite similar mainly in the region of shorter offsets when $z<1.2$ km,
even though the L-BFGS uses more information (diagonal Hessian approximation)
than the HMC to deliver smoother models.

For the parallel Bayesian inversion, the monitor posterior distribution
$\rho_{\text{M}}(\mathbf{m}|\mathbf{d})$ was obtained from the same prior
employed in the baseline estimation. However, the gain of information
regarding the geological model due to baseline inversion can be useful to
geometrically constrain the monitor estimation once the time-lapse changes are
small and spatially localised. Therefore we use the baseline estimates to
construct a prior distribution for the monitor survey in the sequential
approach. Although the baseline posterior distribution
$\rho_{\text{B}}(\mathbf{m}|\mathbf{d})$ possesses information concerning
higher statistical moments, we define the monitor prior distribution
$\rho_{\text{M}}(\mathbf{m})$ as a Gaussian approximation of the latter, that
is,

$\displaystyle\rho_{\text{M}}(\mathbf{m})$
$\displaystyle\equiv\rho_{\text{B}}(\mathbf{m}|\mathbf{d})\,,$
$\displaystyle\approx\exp{\left[-\frac{1}{2}(\mathbf{m}-\mathbf{m}_{\text{B}}^{\text{mean}})^{\text{T}}\Sigma_{\text{B}}^{-1}(\mathbf{m}-\mathbf{m}_{\text{B}}^{\text{mean}})\right]}\,,$
(12)

where $\mathbf{m}_{\text{B}}^{\text{mean}}$ is the sample mean and
$\Sigma_{\text{B}}$ is the sample baseline covariance matrix. In our
calculations, only the diagonal elements of $\Sigma_{\text{B}}$ were
considered. Additionally, we also use the prior baseline information to design
a momentum distribution $\rho(\mathbf{p})$ where initial masses for model
parameters located in the region around the reservoirs (black squares in
Figure 1b) were differently selected. We emphasize that our sequential
approach differs from [47] since the authors use the last few baseline samples
as starting points to the monitor inversion instead of properly designing a
new prior distribution.

###  III.3 Repeatable Acquisition Geometry

Figure 4: The mean (left) and the standard deviation (right) estimates of the
time-lapse model with perfect acquisition geometry using the parallel (a, b)
and sequential (c, d) Bayesian inversions.

We compare the time-lapse inversion results using the parallel and sequential
approaches in the perfect acquisition geometry scenario for both deterministic
and probabilistic cases. The time-lapse posterior distribution was constructed
using random pairs of baseline and monitor samples in both cases. Once again,
we limit ourselves to investigating only the mean and standard deviation
behaviour of the distribution, which are shown in Figure 4. In the two
strategies, the reservoirs are present in the time-lapse mean model, but the
right reservoir is well-resolved in comparison with the left one. The small-
scale structures (decoherent pixels with non-zero velocity values) that appear
in the parallel scheme (Figure 4a) can be associated with differences between
baseline and monitor samples that do not cancel out on average. This kind of
inversion artefact was recently reported by [47] and also appears in the
deterministic inversion (Figure 3c) of previous studies [20, 7]. On the other
hand, it is possible to identify clearly the velocity changes in the reservoir
region when we employ the sequential strategy (Figure 4c), even though it is
difficult to recover the left-side reservoir in that case. These aspects of
the sequential procedure also emerge in deterministic cases (Figure 3d) and
can be driven by affordable regularisation terms in the misfit function [20].

Figure 5: Estimated covariance between baseline and monitor samples in the
parallel (left) and sequential (right) strategies with perfect (a, b) and
perturbed (c, d) acquisition geometry.

The mean estimates (Figure 4a, c) are statistically meaningless without proper
uncertainty quantification of results. The standard deviation
$\sigma_{\text{TL}}$ of the time-lapse changes (Figures 4b, d) was computed
using the well-known first-order error propagation formula

$\sigma_{\text{TL}}=\sqrt{\sigma_{\text{M}}^{2}+\sigma_{\text{B}}^{2}-2\sigma_{\text{MB}}}\,,$
(13)

where $\sigma_{\text{M}}$ and $\sigma_{\text{B}}$ are the standard deviation
associated with the monitor and baseline inversion, respectively, and
$\sigma_{\text{MB}}$ stands for the diagonal element of the covariance matrix
between samples. In the light of Eq. 13, errors in the baseline estimation
will always be reflected as uncertainties in time-lapse images. Moreover, the
propagated uncertainty also depends on the correlation degree between baseline
and monitor samples. In this regard, uncertainties are lower in the parallel
strategy, especially in the deeper part of the model due to the strong
correlation with the baseline samples (Figure 5a). This occurs because the
sampling is performed from the same prior distribution. In the sequential
strategy, monitor samples have lower uncertainty but weaker correlation
(Figure 5b). Therefore, the time-lapse images have similar uncertainties for
both strategies.

###  III.4 Nonrepeatable Acquisition Geometry

Figure 6: The mean (left) and the standard deviation (right) estimates of the
time-lapse model with perturbed acquisition geometry using the parallel (a, b)
and sequential (c, d) Bayesian inversions.

We now consider the performance of the Bayesian inversion approaches when the
acquisition geometry of the monitor survey differs from the baseline one. We
first note that the correlation maps between monitor and baseline samples
(Figure 5c and d) are not sensitive to the nonrepeatable geometry regardless
of the strategy used. This is evidence that the correlation between different
surveys depends more on prior choice in model space $\mathbb{M}$ than slight
changes in data space $\mathbb{D}$. Such similarities also appear in the
standard deviation maps (Figure 6c and d) for the parallel and sequential
schemes, indicating that the estimated uncertainties to the time-lapse model
are less affected by changes in source locations. However, the estimated mean
time-lapse model in the parallel strategy (Figure 6a) presents more artefacts
when compared to the perfect geometry case (Figure 4a), where the former has
the presence of negative anomalies along the model that are paramount in the
deeper part ($z<1.0$ km). As it is well known in the deterministic framework
[20, 7], the sequential scheme and their extensions are more robust to
nonrepeatability issues and as shown in Figure 6c, the reservoirs can be
better identified using this approach.

##  IV Discussion

We have demonstrated in this work that exploring the posterior information of
the baseline survey can benefit the estimation of well-resolved time-lapse
changes. The occurrence of small-scale structures when the inversion is
performed in parallel can be misunderstood as being geomechanical changes that
blur the true velocity variation in the over and underburden regions. This
behaviour for the time-lapse estimates (Figure 4a, b) also were identified by
[47] using the sSVGD method. On the other hand, we argue that the knowledge
regarding the baseline model can be used to design a suitable prior
distribution in the monitor estimation, enabling better reservoir
identifications.

The strategies presented here depend on the quality and convergence of the
baseline estimation. The correlation between samples (Figure 5) shows that
this dependence occurs differently in both methodologies, and it is crucial to
decrease the uncertainties in the parallel when compared with the sequential
inversion. Therefore, we observe a trade-off between the two approaches in
terms of mean estimates and uncertainties. In fact, the results cannot be
useful in practice when the uncertainties have the same magnitude order of the
velocity changes. However, we also argue that the error in the baseline
estimates can be alleviated by considering (1) a burn-in period that
eliminates the initial chain model dependence and produces sample models with
lower dispersion; (2) more samples that improve the convergences at the price
of increasing the computational cost; (3) well-log information to better
constrain vertical resolution the parameters in the reservoir’s region or (4)
data with reduced noise. Moreover, the sequential strategy is more robust to
nonrepeatable errors, even though both inversion methods produce similar
uncertainty estimations in perfect and perturbed acquisition geometry.
Therefore, quantifying the uncertainty can help us discriminate inversion
artefacts from true changes and make better decisions concerning reservoir
management.

##  V Conclusion

In this study, we propose a Bayesian sequential approach to time-lapse
estimation using the Full Waveform Inversion (FWI). The idea involves
integrating the baseline posterior information as prior knowledge in the
monitor evaluation. Moreover, we use the Hamiltonian Monte Carlo (HMC) as a
sampling method combined with an appropriate mass matrix design. We compare
the proposed strategy with the parallel Bayesian inversion, where baseline and
monitor estimations are conducted separately. These comparisons are performed
with perfect and perturbed acquisition geometry, using a time-lapse Marmousi
model as a study case. We concluded that the sequential strategy produces
decorrelated samples between surveys and can provide more accurate time-lapse
estimates with uncertainties of the same order as the parallel methodology.

## Authors Contributions

P. D. S. de Lima initiated the study, performed the simulations, analyzed the
results, co-wrote, and reviewed the manuscript. M. S. Ferreira, G. Corso and
J. M. de Araújo analyzed the results, co-wrote, and reviewed the manuscript.

###### Acknowledgements.

The authors gratefully acknowledge support from Shell Brasil through the “New
Computationally Scalable Methodologies for Target-Oriented 4D Seismic in Pre-
Salt Reservoirs” project at Universidade Federal do Rio Grande do Norte (UFRN)
and the strategic importance of the support given by ANP through the R&D levy
regulation. The authors thank the High-Performance Computing Center (NPAD) at
UFRN for providing computational resources that made the present work
possible. P. D. S. de Lima acknowledges the CAPES/PRINT program (process no.
88887.838338/2023-00). G. Corso acknowledges the CNPq Brazilian research
agency for funding (grant 307907/2019-8). J. M. de Araújo acknowledges the
CNPq Brazilian research agency for funding (grant 311589/2021-9).

## Data Availability Statement

Data associated with this research are available and can be obtained by
contacting the corresponding author.

## References

  * Lumley [2001] D. E. Lumley, Geophysics 66, 50 (2001).
  * de Oliveira Werneck  _et al._ [2022] R. de Oliveira Werneck, R. Prates, R. Moura, M. M. Gonçalves, M. Castro, A. Soriano-Vargas, P. Ribeiro Mendes Júnior, M. M. Hossain, M. F. Zampieri, A. Ferreira, A. Davólio, D. Schiozer, and A. Rocha, Journal of Petroleum Science and Engineering 210, 109937 (2022).
  * Li  _et al._ [2021a] D. Li, S. Peng, Y. Guo, Y. Lu, and X. Cui, International Journal of Greenhouse Gas Control 108, 103336 (2021a).
  * Nakata  _et al._ [2022] R. Nakata, U.-G. Jang, D. Lumley, T. Mouri, M. Nakatsukasa, M. Takanashi, and A. Kato, Geophysical Research Letters 49, e2022GL098734 (2022).
  * Tarantola [1984] A. Tarantola, Geophysics 49, 1259 (1984).
  * Virieux and Operto [2009] J. Virieux and S. Operto, Geophysics 74, WCC1 (2009).
  * Zhou and Lumley [2021a] W. Zhou and D. Lumley, Geophysics 86, 1 (2021a).
  * Lumley  _et al._ [2005] D. Lumley, D. C. Adams, M. Meadows, S. Cole, and R. Wright, in  _SEG Technical Program Expanded Abstracts 2003_ (2005) pp. 1394–1397.
  * Plessix  _et al._ [2010] R. Plessix, G. Baeten, J. W. de Maag, M. Klaassen, Z. Rujie, and T. Zhifei, in  _SEG Technical Program Expanded Abstracts 2010_ (2010) pp. 930–934.
  * Watanabe  _et al._ [2005] T. Watanabe, S. Shimizu, E. Asakawa, and T. Matsuoka, in  _SEG Technical Program Expanded Abstracts 2004_ (2005) pp. 2323–2326.
  * Zhang and Huang [2013] Z. Zhang and L. Huang, Geophysics 78, R259 (2013).
  * Yang  _et al._ [2015] D. Yang, M. Meadows, P. Inderwiesen, J. Landa, A. Malcolm, and M. Fehler, Geophysics 80, M129 (2015).
  * Yang  _et al._ [2016] D. Yang, F. Liu, S. Morton, A. Malcolm, and M. Fehler, Geophysics 81, R225 (2016).
  * Zheng  _et al._ [2012] Y. Zheng, P. Barton, and S. Singh, in  _SEG Technical Program Expanded Abstracts 2011_ (2012) pp. 4195–4200.
  * Maharramov and Biondi [2014] M. Maharramov and B. Biondi, in  _SEG Technical Program Expanded Abstracts 2014_ (2014) pp. 954–959.
  * Maharramov and Biondi [2015] M. Maharramov and B. Biondi, in  _77th EAGE Conference and Exhibition 2015_ (2015) pp. 1–5.
  * Maharramov  _et al._ [2016] M. Maharramov, B. L. Biondi, and M. A. Meadows, Geophysics 81, R485 (2016).
  * Routh  _et al._ [2012] P. Routh, G. Palacharla, I. Chikichev, and S. Lazaratos, in  _SEG Technical Program Expanded Abstracts 2012_ (2012) pp. 1–6.
  * Asnaashari  _et al._ [2012] A. Asnaashari, R. Brossier, S. Garambois, F. Audebert, P. Thore, and J. Virieux, in  _SEG Technical Program Expanded Abstracts 2012_ (2012) pp. 1–5.
  * Asnaashari  _et al._ [2015] A. Asnaashari, R. Brossier, S. Garambois, F. Audebert, P. Thore, and J. Virieux, Geophysical Prospecting 63, 78 (2015).
  * Kamei and Lumley [2017] R. Kamei and D. Lumley, Geophysical Journal International 209, 1239 (2017).
  * Zhou and Lumley [2021b] W. Zhou and D. Lumley, Geophysics 86, R161 (2021b).
  * Mardan  _et al._ [2023] A. Mardan, B. Giroux, and G. Fabien-Ouellet, Geophysics 88, R25 (2023).
  * Sambridge and Mosegaard [2002] M. Sambridge and K. Mosegaard, Reviews of Geophysics 40, 3 (2002).
  * Tarantola [2004] A. Tarantola,  _Inverse Problem Theory and Methods for Model Parameter Estimation_ (Society for Industrial and Applied Mathematics, USA, 2004).
  * Sambridge [2013] M. Sambridge, Geophysical Journal International 196, 357 (2013).
  * Sen and Stoffa [2013] M. K. Sen and P. L. Stoffa,  _Global Optimization Methods in Geophysical Inversion_ (Cambridge University Press, 2013).
  * Betancourt [2018] M. Betancourt, A conceptual introduction to hamiltonian Monte Carlo (2018), arXiv:1701.02434 [stat.ME] .
  * Brooks  _et al._ [2011] S. Brooks, A. Gelman, G. Jones, and X.-L. Meng,  _Handbook of Markov Chain Monte Carlo_ (Chapman and Hall / CRC Press, 2011).
  * Duane  _et al._ [1987] S. Duane, A. Kennedy, B. J. Pendleton, and D. Roweth, Physics Letters B 195, 216 (1987).
  * Sen and Biswas [2017] M. K. Sen and R. Biswas, Geophysics 82, R119 (2017).
  * Fichtner and Simutė [2018] A. Fichtner and S. Simutė, Journal of Geophysical Research: Solid Earth 123, 2984 (2018).
  * Fichtner  _et al._ [2018] A. Fichtner, A. Zunino, and L. Gebraad, Geophysical Journal International 216, 1344 (2018).
  * Gebraad  _et al._ [2020] L. Gebraad, C. Boehm, and A. Fichtner, Journal of Geophysical Research: Solid Earth 125, e2019JB018428 (2020).
  * Aleardi and Salusti [2020] M. Aleardi and A. Salusti, Geophysics 85, R177 (2020).
  * Aleardi  _et al._ [2020] M. Aleardi, A. Salusti, and S. Pierini, Near Surface Geophysics 18, 515 (2020).
  * de Lima  _et al._ [2023] P. D. S. de Lima, G. Corso, M. S. Ferreira, and J. M. de Araújo, Physica A: Statistical Mechanics and its Applications 617, 128618 (2023).
  * Eikrem  _et al._ [2019] K. S. Eikrem, G. Nævdal, and M. Jakobsen, Geophysical Prospecting 67, 379 (2019).
  * Thurin  _et al._ [2019] J. Thurin, R. Brossier, and L. Métivier, Geophysical Journal International 219, 1613 (2019).
  * Huang  _et al._ [2020] X. Huang, K. S. Eikrem, M. Jakobsen, and G. Nævdal, Geophysics 85, C125 (2020).
  * Gineste  _et al._ [2019] M. Gineste, J. Eidsvik, and Y. Zheng, Geophysics 85, R29 (2019).
  * Zhang  _et al._ [2020] X. Zhang, C. Roy, A. Curtis, A. Nowacki, and B. Baptie, Geophysical Journal International 222, 1639 (2020).
  * Zhang and Curtis [2020a] X. Zhang and A. Curtis, Journal of Geophysical Research: Solid Earth 125, e2019JB018589 (2020a).
  * Zhang and Curtis [2020b] X. Zhang and A. Curtis, Geophysical Journal International 222, 406 (2020b).
  * Zhang and Curtis [2021] X. Zhang and A. Curtis, Geophysics 86, A45 (2021).
  * Zhang  _et al._ [2023] X. Zhang, A. Lomas, M. Zhou, Y. Zheng, and A. Curtis, Geophysical Journal International 234, 546 (2023).
  * Zhang and Curtis [2023] X. Zhang and A. Curtis, Bayesian variational time-lapse full-waveform inversion (2023), arXiv:2308.08805 [physics.geo-ph] .
  * Kotsi  _et al._ [2020a] M. Kotsi, A. Malcolm, and G. Ely, in  _SEG Technical Program Expanded Abstracts 2020_ (2020) pp. 845–849.
  * Kotsi  _et al._ [2020b] M. Kotsi, A. Malcolm, and G. Ely, Geophysical Journal International 222, 1245 (2020b).
  * Berti  _et al._ [2023] S. Berti, M. Aleardi, and E. Stucchi, Geophysical Prospecting (2023).
  * De Lima  _et al._ [2023] P. D. De Lima, M. S. Ferreira, G. Corso, and J. M. de Araújo, in  _5th EAGE Conference on Petroleum Geostatistics_ (2023) pp. 1–5.
  * Bayes and Price [1763] M. Bayes and M. Price, Philosophical Transactions of the Royal Society of London 53, 370 (1763).
  * Bodin  _et al._ [2012] T. Bodin, M. Sambridge, N. Rawlinson, and P. Arroucau, Geophysical Journal International 189, 1536 (2012).
  * Brossier  _et al._ [2010] R. Brossier, S. Operto, and J. Virieux, Geophysics 75, R37 (2010).
  * Métivier  _et al._ [2016] L. Métivier, R. Brossier, Q. Mérigot, E. Oudet, and J. Virieux, Geophysical Journal International 205, 345 (2016).
  * Liu  _et al._ [2016] Y. Liu, J. Teng, T. Xu, Y. Wang, Q. Liu, and J. Badal, Geophysical Journal International 209, 106 (2016).
  * Carvalho  _et al._ [2021] P. T. C. Carvalho, S. L. E. F. da Silva, E. F. Duarte, R. Brossier, G. Corso, and J. M. de Araújo, Geophysical Journal International 229, 35 (2021).
  * Buland and Ouair [2006] A. Buland and Y. E. Ouair, Geophysics 71, R43 (2006).
  * Lemos [2018] N. A. Lemos,  _Analytical Mechanics_ (Cambridge University Press, 2018).
  * Metropolis  _et al._ [1953] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, The Journal of Chemical Physics 21, 1087 (1953).
  * Leimkuhler and Reich [1994] B. Leimkuhler and S. Reich,  _Simulating Hamiltonian Systems_ (Cambridge University Press, 1994).
  * Fichtner  _et al._ [2021] A. Fichtner, A. Zunino, L. Gebraad, and C. Boehm, Geophysical Journal International 227, 941 (2021).
  * Asnaashari  _et al._ [2013] A. Asnaashari, R. Brossier, S. Garambois, F. Audebert, P. Thore, and J. Virieux, Geophysics 78, R25 (2013).
  * Li  _et al._ [2021b] D. Li, S. Peng, X. Huang, Y. Guo, Y. Lu, and X. Cui, International Journal of Greenhouse Gas Control 110, 103417 (2021b).
  * Plessix [2006] R. Plessix, Geophysical Journal International 167, 495 (2006).
  * Izzatullah  _et al._ [2021] M. Izzatullah, T. van Leeuwen, and D. Peter, Geophysical Journal International 227, 1523 (2021).
  * Liu and Nocedal [1989] D. C. Liu and J. Nocedal, Mathematical programming 45, 503 (1989).