# A practical and efficient approach for Bayesian reservoir inversion:
Insights from the Alvheim field data

Karen S. Auestad  The Tien Mai<EMAIL_ADDRESS>Mina Spremic  Jo Eidsvik
<EMAIL_ADDRESS>
###### Abstract

Stochastic reservoir characterization, a critical aspect of subsurface
exploration for oil and gas reservoirs, relies on stochastic methods to model
and understand subsurface properties using seismic data. This paper addresses
the computational challenges associated with Bayesian reservoir inversion
methods, focusing on two key obstacles: the demanding forward model and the
high dimensionality of Gaussian random fields. Leveraging the generalized
Bayesian approach, we replace the intricate forward function with a
computationally efficient multivariate adaptive regression splines method,
resulting in a 34 acceleration in computational efficiency. For handling high-
dimensional Gaussian random fields, we employ a fast Fourier transform (FFT)
technique. Additionally, we explore the preconditioned Crank–Nicolson method
for sampling, providing a more efficient exploration of high-dimensional
parameter spaces. The practicality and efficacy of our approach are tested
extensively in simulations and its validity is demonstrated in application to
the Alvheim field data.

###### keywords:

Bayesian inversion; MCMC; Generalized Bayesian approach; spatial correlation;
fast Fourier transform (FFT); Stochastic reservoir.

organization= Department of Mathematical Sciences, Norwegian University of
Science and Technology,city=Trondheim, postcode=7034, country=Norway

##  1 Introduction

Stochastic reservoir characterization from seismic data involves using
stochastic methods and techniques to model and characterize subsurface
reservoirs based on seismic data, see e.g. [8] and [24, 22] for an overview of
the problem. Seismic data are measurements of the reflections of sound waves
or vibrations sent into the subsurface layers of the Earth. These reflections
provide information about the geological structure and properties of the
subsurface, including the presence and characteristics of hydrocarbon
reservoirs. The goal of stochastic reservoir characterization is to provide a
range of plausible subsurface models that honor the available data and capture
the uncertainty associated with the reservoir properties. This information is
crucial for making informed decisions in the exploration and development of
oil and gas reservoirs, as it helps in assessing risks, optimizing well
placement, and designing efficient production strategies.

Stochastic reservoir modeling entails the representation of reservoir
properties, such as oil, gas, and clay, as random fields. This involves
assigning probability distributions to these properties to capture the spatial
variability within the reservoir. Bayesian inversion has emerged as a well-
suited method for reservoir characterization challenges, as it enables the
integration of prior knowledge with observational data. In many cases, there
is a need for a priori knowledge of Earth parameters since the information
obtained from measurements alone may be insufficient. This a priori knowledge
is encapsulated by a prior distribution, allowing for the incorporation of
trends and uncertainties in our initial beliefs [36, 46, 17, 21]. In Bayesian
inversion, a likelihood function, grounded in geophysical principles, informs
us about which Earth parameter values align with the observed data. By
amalgamating the prior distribution and the likelihood function, a posterior
distribution of Earth parameters can be formed. Although, examining posterior
densities in Bayesian inversion can be analytically challenging, various
methods are available for approximating the posterior distribution [39, 31,
32, 23].

The practical implementation of Bayesian reservoir inversion methods often
grapples with computational challenges, primarily stemming from two key
sources. Firstly, the forward model or function integral to the inversion
process is frequently demanding in terms of computational resources. Secondly,
in order to faithfully represent the spatial variability within a reservoir, a
random field, typically a Gaussian random field, is utilized. However, the
utilization of such random fields results in a high-dimensional distribution
that necessitates extensive sampling and evaluation efforts.

To illustrate geophysical inversion, we examine seismic amplitude data from
the Alvheim field. In this specific study, a Gaussian random field is employed
to capture reservoir variability extending into a dimension exceeding 44,000.
This considerable dimensionality adds to the computational complexity,
intensifying the challenges associated with sampling and evaluating the
distribution. Consequently, it presents a noteworthy hurdle in the effective
application of Bayesian reservoir inversion methods.

Another significant challenge is the design of an efficient Markov Chain Monte
Carlo (MCMC) method in high-dimensional settings. This challenge is
particularly crucial as Bayesian methods face a daunting aspect: the slow
convergence of the sampling of the posterior distribution, which becomes
increasingly sluggish as the dimensionality escalates. This poses a notable
obstacle in the application of Bayesian reservoir inversion methods,
emphasizing the need for innovative solutions to address the intricacies of
high-dimensional parameter spaces.

In this study, our primary objective is to address and mitigate the
computational challenges associated with Bayesian reservoir inversion methods.
The initial computational complexity arising from the evaluation of the
forward model is effectively tackled by employing a cutting-edge statistical
approach known as the generalized Bayesian approach. This approach, as
detailed in abundant literature [7, 25, 35, 3, 30, 33, 37, 34, 27, 2],
represents a generalization of classical Bayesian methods. Instead of relying
on the true likelihood, this approach allows the use of a fractional
likelihood or a loss-based function. The generalized Bayesian approach has
demonstrated efficacy both theoretically and in practical applications. To
leverage this novel insight, we address the first computational challenge by
substituting the intricate forward function with an estimated function derived
from training data. Specifically, we employ the multivariate adaptive
regression splines (MARS), [18] for this purpose. Notably, our implementation
of this approach yields a noteworthy acceleration in computational efficiency.
In our case study of the Alvheim field data, the reported speedup is
approximately 34 times faster when compared to the use of the true forward
function. This substantial improvement in running time underscores the
practical utility of the generalized Bayesian approach in enhancing the
computational efficiency of Bayesian reservoir inversion methods.

To address the computational challenges posed by the substantial
dimensionality of the Gaussian random field, we adopt a strategy that
leverages the efficiency of the Fast Fourier transform (FFT) for simulating
and evaluating the Gaussian random field, as outlined by the works [10, 14,
45, 1]. This technique involves the circulant embedding of a structured
covariance matrix onto a torus. The implementation of this method not only
significantly reduces the required memory but also yields remarkably fast
computations. Notably, this method proves to be highly efficient, enabling
simulations and evaluations on standard computing hardware. In our case study
focused on the Alvheim data, which involves a Gaussian random field of
dimensionality exceeding 44,000, the computations are easily managed on a
laptop with 16GB RAM. Furthermore, the simulation time for a vector within
this high-dimensional space is impressively less than a second. This
underscores the practical feasibility and efficiency of employing FFT-based
methods to handle and streamline the simulation and evaluation of Gaussian
random fields in large-scale applications.

The impetus behind the present study is rooted in the recognition that a more
tailored sampling methodology, specifically designed for the high
dimensionality, has the potential to yield substantial computational speed
enhancements. Notably, [13] introduced a groundbreaking approach to Markov
Chain Monte Carlo (MCMC) sampling that addresses the notorious ’curse of
dimensionality’ associated with Gaussian priors, given appropriate assumptions
on the likelihood function. This innovative method, termed “preconditioned
Crank–Nicolson” (pCN), represents a modification of the standard random-walk
Metropolis sampling technique. The key feature of pCN involves scaling the
position from the previous iteration before adding a random shift and
generating the proposal. This adjustment effectively mitigates the
computational challenges associated with high-dimensional spaces. By doing so,
pCN provides a more efficient and effective means of exploring the parameter
space, particularly in scenarios where traditional methods may encounter
difficulties due to the curse of dimensionality [26]. The introduction of pCN
thus opens up new avenues for achieving significant computational speedups in
the context of MCMC sampling, offering a promising solution to challenges
posed by increased system dimensionality. Utilizing the pCN method within our
framework, where the proposed approach involves incorporating a physically
informed prior, enables effective exploration of various high-probability
regions within the intricate posterior distribution. This aspect holds
particular significance in the context of Bayesian inverse problems, as
emphasized in previous works such as [39, 31, 32].

In application to the Alvheim field data, our study reveals varying saturation
levels around wells, with more gas near gas-producing wells and more oil near
oil-producing ones. Uncertainty is low near these wells but higher in
transition zones. Compared to [46], our study shows similar high saturation
levels near wells but differs in saturation distributions and uncertainty
levels, with our MCMC method generally yielding lower uncertainty. We also
note higher saturation levels in clay content. Ternary plots highlight
differences, especially in shallow areas, where our results suggest higher oil
saturation and less uncertainty. While [46] use local ensemble-based
techniques, their effectiveness in reducing uncertainty remains unclear.

The subsequent sections of the paper are organized as follows: Section 2
provides an introduction to the Bayesian rock physics problem and offers
insights into the Alvheim field data. Our proposed procedure is outlined in
Section 3. The findings from the Alvheim field real data are presented and
deliberated upon in Section 5. Further, Section 4 encompasses the presentation
of numerical studies. For additional information on the fast Fourier transform
sampling method and more insights into the simulation studies, refer to A, B
and C. We discuss and conclude our work in Section 6.

##  2 Background

###  2.1 The Alvheim field

The Alvheim field, situated in the North Sea, Figure 1, off the coast of
Norway and dating back to the Paleocene era, is a turbiditic oil and gas
reservoir. The field’s characteristics include a complex sand distribution and
hydrocarbon trapping system, featuring multiple submarine fan lobes and
significant variations in depositional facies, net-to-gross ratios, and sand
textures. The sands range from massive, thick-bedded formations to more
diverse, inter-bedded sand-shale intervals. The overall depositional pattern
is influenced by structural topography linked to deeper faults and salt
tectonics.

Figure 1: Location of the Alvheim field in the North Sea. Map from:
https://www.norskpetroleum.no/en/interactive-map-quick-downloads/interactive-
map/

This field poses unique challenges for seismic reservoir characterization due
to its complex geological features. The reservoir sands, attributed to the
Heimdal Formation and capped by Lista Formation shales, exhibit intricate rock
physics properties influenced by a compaction history related to burial depths
of approximately 2 kilometers below the seafloor and temperatures around 70°C.
Diagenetic processes, such as the transformation of smectite-rich shales to
illite and cementation of quartzose sands, contribute to significant elastic
stiffening of the rock frame. The presence of both unconsolidated sands and
cemented sandstones further complicates the fluid sensitivity of seismic
signals, see e.g [4]. Additionally, the distribution of hydrocarbons,
encompassing both gas and oil, is complex and poorly understood, likely
influenced by regional tectonic events. These challenges underscore the
difficulty in accurately characterizing the Alvheim field using seismic
methods, requiring careful consideration of geological intricacies and fluid
dynamics.

(a) The depth at the Alvheim field.

(b) Zero offset data, $R_{0}$.

(c) Gradient data, $G$.

Figure 2: AVO data presented in the inline-crossline coordinate system for the
Alvheim field. Four wells are marked by coloured circles, where the colour
indicates the majority type of hydrocarbons found at that well. The wells are
located at $($inline,crossline$)$-locations $(484,4798)$, $(716,4914)$,
$(1026,4924)$ and $(1300,5150)$.

####  2.1.1 Seismic AVO data

This investigation employs seismic amplitude-versus-offset (AVO) data acquired
from the Alvheim field, covering angle stacks within near ($12^{\circ}$), mid
($22^{\circ}$), and far ($31^{\circ}$) angle ranges. The data have undergone
processing and preconditioning to ensure their appropriateness for
quantitative analysis, following the methodology outlined by [43]. AVO
attributes, such as zero-offset reflectivity ($R_{0}$) and AVO gradient ($G$),
have been extracted and calibrated using scalars derived from upscaled well
log data [5]. Each inline and crossline not only includes two data points
($R_{0}$ and $G$) but also incorporates associated traveltimes, as shown in
Figure 2(a), subsequently converted to depth.

We focuses on the top-reservoir horizon, where, in addition to spatial
variations in the inline-crossline domain, depth information is available for
the 2D slice at each inline-crossline location. Referred to as a 2.5D inverse
problem, the AVO data for the Alvheim field, depicted in Figure 2, encompass
$R_{0}$ and $G$ for each inline and crossline of the top-reservoir. Covering
an area of approximately 12.35km$\times$8.85km, the dataset is subsampled by
selecting measurements from every 4 inlines and crosslines, equivalent to 50m
intervals in both crossline and inline directions. The spatial grid
representation in Figure 2(a) is utilized for modeling, maintaining the
original inline and crossline units without conversion to distance metrics.

####  2.1.2 Well log data

This study utilizes data from four wells, namely 24/6–2, 24/6–4, 25/4–7, and
25/4–8, whose locations and outcomes at the top-reservoir are illustrated in
Figure 2(a). These well details will be referred to in figures throughout the
paper, although the legends will be omitted. All these wells are situated
within the seismic inline-crossline ranges and drilled on structural highs in
the shallower sections of the top-reservoir. The well log data encompass sonic
velocities ($V_{\text{p}}$ and $V_{\text{s}}$), densities, and petrophysical
logs, including fluid saturations (oil, gas, and brine) and clay content, as
outlined in [43] and [46].

In the case of well 24/6–2, commercial gas is found in the Heimdal Formation,
with a gas column of 52m down to the gas-oil contact at 2151m. Additionally,
there is a thin oil column of 17m down to the oil-water contact at 2168m. Well
25/4–7 encountered commercial oil saturation, with an oil column of 48m down
to the oil-water contact at 2133m. Hydrocarbon saturations generally fall
within commercial ranges (0.6–0.9), although lower values are observed in
heterolithic zones where the presence of clay influences reservoir transport
properties. Thin clay laminations may also induce patchy saturations at the
log scale. For Wells 24/6–4 and 25/4–8, constant interpretations are made for
different zones due to missing saturation logs estimated from resistivity
logs.

###  2.2 Bayesian model

####  2.2.1 Probabilistic geophysical inversion

The inverse problem arises when there is the ability to collect data, denoted
as ${\boldsymbol{y}}$, which are associated with an underlying variable of
significance referred to as ${\boldsymbol{x}}$. The relationship between the
variable of interest and the observed data is expressed as:

${\boldsymbol{y}}=h({\boldsymbol{x}})+{\boldsymbol{e}},$ (1)

where $h$ defines a recognized functional relationship and ${\boldsymbol{e}}$
represents stochastic perturbations.

The computational challenge of inversion is multifaceted, primarily stemming
from the intricate interplay of high-dimensional, nonlinear, and uncertain
data. In many real-world scenarios, the variables of interest and observed
data inhabit high-dimensional spaces, necessitating sophisticated algorithms
to explore parameter space efficiently. Nonlinear function $h$ further
compounds the challenge, requiring iterative methods to converge on optimal
solutions. Moreover, uncertainty pervades the inversion process, demanding
robust techniques such as Bayesian inference to quantify and incorporate
uncertainty into parameter estimates. [47].

Bayesian inversion involves combining a likelihood model for observed data
with prior knowledge about the variable of interest to generate a posterior
distribution $f({\boldsymbol{x}}|{\boldsymbol{y}})$. Prior knowledge is
expressed through a prior distribution, denoted $f({\boldsymbol{x}})$. The
likelihood probability density function,
$f({\boldsymbol{y}}|{\boldsymbol{x}})$, indicates the likelihood of obtaining
the observed data given a specific value of ${\boldsymbol{x}}$. According to
Bayes’ rule, the posterior distribution function can be written as:
$f({\boldsymbol{x}}|{\boldsymbol{y}})\propto
f({\boldsymbol{y}}|{\boldsymbol{x}})f({\boldsymbol{x}}),$ where the posterior
$f({\boldsymbol{x}}|{\boldsymbol{y}})$ is referred to as
$\pi({\boldsymbol{x}})$.

####  2.2.2 Observation model

Before delving into the above mentioned models, we provide essential notation,
adopting conventions akin to those in [46]. Reservoir variables are denoted at
locations ${\boldsymbol{u}}\in\mathcal{D}$, with $\mathcal{D}$ representing
the inline and crossline grid used to portray seismic data at the top
reservoir. The grid dimensions are $248\times 178$, yielding $N=44,144$ grid
cells with inline-crossline locations
${\boldsymbol{u}}_{1},\ldots,{\boldsymbol{u}}_{N}$. Seismic AVO data at
locations ${\boldsymbol{u}}\in\mathcal{D}$ are

${\boldsymbol{y}}({\boldsymbol{u}})=[R_{0}({\boldsymbol{u}}),G({\boldsymbol{u}})].$

The main interest lies in the oil and gas saturation, and clay content at each
of the inline-crossline locations. Results will be displayed for the reservoir
variables of interest ${\boldsymbol{r}}({\boldsymbol{u}})$, which we denote by
${\boldsymbol{r}}({\boldsymbol{u}})=[{\boldsymbol{S}}_{\text{g}}({\boldsymbol{u}}),{\boldsymbol{S}}_{\text{o}}({\boldsymbol{u}}),{\boldsymbol{V}}_{\text{clay}}({\boldsymbol{u}})]$.

For the purpose of efficient computations, we work with transformed
saturations and clay content. Variables are transformed back and forth by
applying logistic functions at all locations ${\boldsymbol{u}}\in\mathcal{D}$,

$\displaystyle S_{l}({\boldsymbol{u}})$ $\displaystyle=$
$\displaystyle\frac{\exp\left[x_{l}({\boldsymbol{u}})\right]}{1+\exp\left[x_{\text{g}}({\boldsymbol{u}})\right]+\exp\left[x_{\text{o}}({\boldsymbol{u}})\right]},\,l=\text{o,g}\mbox{(oil,gas)}$
$\displaystyle S_{\text{b}}({\boldsymbol{u}})$ $\displaystyle=$
$\displaystyle\frac{1}{1+\exp\left[x_{\text{g}}({\boldsymbol{u}})\right]+\exp\left[x_{\text{o}}({\boldsymbol{u}})\right]},$

meaning that
$S_{\text{g}}({\boldsymbol{u}})+S_{\text{o}}({\boldsymbol{u}})+S_{\text{b}}({\boldsymbol{u}})=1$
for every location ${\boldsymbol{u}}$, where $x_{\text{o}}({\boldsymbol{u}})$
and $x_{\text{g}}({\boldsymbol{u}})\in\mathbb{R}$. In doing so, we keep brine
as a reference or background saturation. Similar for clay content
$V_{\text{clay}}({\boldsymbol{u}})$,

$V_{\text{clay}}({\boldsymbol{u}})=\frac{\exp\left[x_{\text{clay}}({\boldsymbol{u}})\right]}{1+\exp\left[x_{\text{clay}}({\boldsymbol{u}})\right]},$

where $x_{\text{clay}}({\boldsymbol{u}})\in\mathbb{R}$. Altogether, we then
get the reservoir variables of interest from a logistic function, compactly
denoted by

${\boldsymbol{r}}({\boldsymbol{u}})={\boldsymbol{f}}_{0}\left[{\boldsymbol{x}}({\boldsymbol{u}})\right],\hskip
8.53581pt {\boldsymbol{u}}\in\mathcal{D}.$ (2)

It is note that our Bayesian model adopts fixed hyperparameters in both the
prior and likelihood functions. An alternative approach in a fully Bayesian
formulation involves placing prior distributions on these parameters, as
exemplified in [36]. However, in our specific setting, we opt for a
sensitivity analysis approach, focusing on the examination of the most
relevant input parameters.

At a location ${\boldsymbol{u}}$, we let ${\boldsymbol{y}}({\boldsymbol{u}})$
denote the zero-offset coefficient $R_{0}$ and gradient $G$ at that location.
The model can be written in an additive manner,

${\boldsymbol{y}}({\boldsymbol{u}})={\boldsymbol{h}}_{0}\left[{\boldsymbol{x}}({\boldsymbol{u}})\right]+{\boldsymbol{\epsilon}}({\boldsymbol{u}}),\quad{\boldsymbol{\epsilon}}({\boldsymbol{u}})\sim
N(\mbox{\boldmath{$0$}},{\boldsymbol{\Omega}}_{0}),$ (3)

where
${\boldsymbol{h}}_{0}\left[{\boldsymbol{x}}({\boldsymbol{u}})\right]={\boldsymbol{g}}_{0}\left\\{{\boldsymbol{f}}_{0}\left[{\boldsymbol{x}}({\boldsymbol{u}})\right]\right\\}$,
and ${\boldsymbol{g}}_{0}[{\boldsymbol{r}}({\boldsymbol{u}})]$ represents the
rock physics forward model function applied to the reservoir variables
${\boldsymbol{r}}({\boldsymbol{u}})$ of primary interest. The $2\times 2$
matrix ${\boldsymbol{\Omega}}_{0}$ characterizes the variances in $R_{0}$ and
$G$ and the correlation between the two. Through the utilization of seismic
reflection data in proximity to the four wells and the corresponding well log
data in the reservoir zones, we define the parameters
$\mbox{Var}(R_{0})=0.003$, $\mbox{Var}(G)=0.03$, and
$\mbox{Corr}(R_{0},G)=-0.6$. Notably, these parameter values closely align
with those employed in [16].

Therefore, the likelihood model is characterized by a Gaussian distribution,
with the mean derived from the nonlinear geophysical forward model and a
spatially independent variance. The applicability of geophysical forward model
relationships is assumed for each inline and crossline, considering the known
depth of the top reservoir at each grid location. The likelihood model assumes
conditional independence, implying that data at different locations are
independent given the rock properties. Combining data from all inline-
crossline locations into a length $2N$ vector ${\boldsymbol{y}}$, the
probability density function $p({\boldsymbol{y}}|{\boldsymbol{x}})$ follows a
Gaussian distribution with mean ${\boldsymbol{h}}_{0}({\boldsymbol{x}})$ and
covariance matrix ${\boldsymbol{\Omega}}$. Assuming conditional independence
between AVO data at different inline-crossline grid locations, the $N\times N$
covariance matrix ${\boldsymbol{\Omega}}$ is block diagonal, with the $2\times
2$ matrix ${\boldsymbol{\Omega}}_{0}$ repeating along the block diagonal
elements. The function ${\boldsymbol{h}}({\boldsymbol{x}})$ is constructed
based on the same rock physics relations applied to various spatially varying
input reservoir variables.

The realistic forward model, denoted as ${\boldsymbol{h}}_{0}(\cdot)$, is
constructed by integrating rock physics and geophysical models with geological
properties. This forward model takes as input random fields representing clay
content and fluid saturation. To address the cementation effect presumed to
occur at a specific depth, we incorporate rock physics models based on
mechanical and chemical compaction. It is also noted that the forward model is
considered to be independent of locations and depth. Depending on the depth of
the input reservoir variables, we employ either an unconsolidated sand model
or a contact cement model, as outlined in [6]. Additionally, Gassmann fluid
substitution is implemented, followed by the Shuey approximation, to derive
the amplitude-versus-offset (AVO) attributes. Further details regarding the
forward model can be found in [46].

####  2.2.3 Physical-informed priors

The prior distributions are formulated based on spatial random field models
for the reservoir variables ${\boldsymbol{r}}({\boldsymbol{u}})$ as defined in
equation (2). Independent Gaussian random fields are employed to represent the
spatial structures within the transformed saturation variables
${\boldsymbol{x}}_{\text{g}}=[x_{\text{g}}({\boldsymbol{u}}_{1}),\ldots,x_{\text{g}}({\boldsymbol{u}}_{N})]$,
${\boldsymbol{x}}_{\text{o}}=[x_{\text{o}}({\boldsymbol{u}}_{1}),\ldots,x_{\text{o}}({\boldsymbol{u}}_{N})]$,
and
${\boldsymbol{x}}_{\text{clay}}=[x_{\text{clay}}({\boldsymbol{u}}_{1}),\ldots,x_{\text{clay}}({\boldsymbol{u}}_{N})]$.
To be more specific, given the spatial discretization of locations, we obtain
Gaussian multivariate distributions that,

${\boldsymbol{x}}_{\text{g}}\sim
N({\boldsymbol{\mu}}_{\text{g}},{\boldsymbol{\Sigma}}_{\text{g}}),\hskip
8.53581pt{\boldsymbol{x}}_{\text{o}}\sim
N({\boldsymbol{\mu}}_{\text{o}},{\boldsymbol{\Sigma}}_{\text{o}}),\hskip
8.53581pt{\boldsymbol{x}}_{\text{clay}}\sim
N({\boldsymbol{\mu}}_{\text{clay}},{\boldsymbol{\Sigma}}_{\text{clay}}).$ (4)

In this context, the mean vectors ${\boldsymbol{\mu}}_{\text{g}}$,
${\boldsymbol{\mu}}_{\text{o}}$, and ${\boldsymbol{\mu}}_{\text{clay}}$ have
dimensions of $N$, while the covariance matrices
${\boldsymbol{\Sigma}}_{\text{g}}$, ${\boldsymbol{\Sigma}}_{\text{o}}$, and
${\boldsymbol{\Sigma}}_{\text{clay}}$ are sized $N\times N$. These parameters
are derived from physical data, specifically well log information gathered at
various depths.

Figure 3: Prior mean saturations of gas, oil and brine (from left to right).

The mean values contribute to establishing a smoothly varying background model
with depth. Predominantly high gas saturation is anticipated in the very
shallow inline-crossline locations, elevated oil saturation is expected at
relatively shallow positions, and brine is highly probable to be predominant
at deeper locations. The standard deviation for each variable is assumed to
remain constant with depth.

Regarding the covariance matrix, we incorporate spatial correlation using a
Gaussian correlation function with an effective correlation length of $15$
grid cells. This correlation function is applied to $x_{\text{o}}$,
$x_{\text{g}}$, and $x_{\text{clay}}$ and is specified based on prior
assumptions about the geological depositional environment, particularly
considering geographically separated lobe structures.

The prior mean values of oil, gas and brine are depicted in Figure 3.
Generally, the visualization illustrates that in shallow regions, the prior
mean saturations of gas and oil are elevated. Specifically, the prior mean
saturation of oil is relatively high in the vicinity of the well marked by a
violet circle in Figure 2(a), while the prior mean saturation of gas is
comparatively low in this region. In the deepest areas, such as the top-left
corner, both gas and oil saturations exhibit lower values. In the intermediate
terrain between the middle ridge and the elevated region containing the
predominantly oil-filled well, the prior mean saturation for oil surpasses
that of gas. Moreover, in the shallowest regions, gas exhibits a slightly
higher prior mean saturation than oil.

##  3 Proposed method

###  3.1 Approximate forward function via MARS

Evaluating the forward model in the Alvheim case at a value ${\boldsymbol{x}}$
takes a relatively long time. At each iteration in the MCMC algorithm,
$h({\boldsymbol{x}}^{p})$ needs to be calculated to find the likelihood at the
proposed sample ${\boldsymbol{x}}^{p}$. The goal is to find a good
approximation that considerably reduces computation time.

The idea here is to use some ensemble data to train an approximate model that
is faster to evaluate in compared to the underlying forward model. This idea
maybe similar those approaches based on Kalman filtering but we take advantage
of these ensemble data for training an relatively good approximate function
for the purpose of reducing the cost of calculating the acceptance rate within
the MCMC algorithm. Let ${\boldsymbol{x}}^{i}\sim
N({\boldsymbol{\mu}},{\boldsymbol{\Sigma}}),{\boldsymbol{y}}^{i}={\boldsymbol{h}}({\boldsymbol{x}}^{i}),i=1,...,n$
be some ensemble data. Based on this ensemble data, we train a function
$\hat{{\boldsymbol{h}}}$ that is not only a good approximation of
${\boldsymbol{h}}$ but also faster to evaluate. For this purpose, we make use
of the multivariate adaptive regression splines (MARS), introduced in [18].

MARS constitutes a non-parametric regression method utilized to elucidate
intricate associations among variables. It operates by delineating the
connection between predictor variables and the response variable through a
sequence of basis functions, typically piecewise linear in nature, amalgamated
to compose the ultimate regression model. The algorithm initiates with a
constant prediction and discerns pivotal breakpoints in the predictor
variables, segmenting the data into distinct sections where varied linear
functions are employed. During the forward pass, MARS integrates basis
functions into the model predicated on the identified breakpoints,
encompassing simple linear terms, hinge functions, or interaction terms.
Subsequently, a backward pass ensues wherein MARS prunes redundant basis
functions to streamline the model, thus averting overfitting and enhancing
interpretability. The resultant MARS model embodies a blend of these basis
functions, yielding a piecewise linear regression model adept at capturing the
dataset’s nonlinearities. MARS proves advantageous in datasets exhibiting
nonlinear or interactive relationships, offering flexibility in modeling
without imposing specific functional forms. Its utility extends to scenarios
where conventional linear models fall short in complexity and more intricate
machine learning models risk overfitting. Particularly well-suited for
regression tasks, MARS stands as a potent tool for data exploration and
predictive analytics [28].

The MARS approximation takes the following form:

$\widehat{h}({\boldsymbol{x}})=\beta_{0}+\sum_{d=1}^{D}\beta_{d}g_{d}({\boldsymbol{x}}),$

where $\beta_{1},\beta_{2},...,\beta_{D}$ represent coefficients jointly
determined by minimizing the residual sum of squares. The linear basis
functions $g_{d}({\boldsymbol{x}})$ are in the set
$\mathcal{C}=\\{(x_{j}-c_{j})_{+},(c_{j}-x_{j})_{+}\\}_{c_{j}\in\\{x_{1j},x_{2j},...x_{nj}\\}}$
for each covariate $x_{j}$ and observed value $x_{ij}$ in the dataset, where
$(a)_{+}=\max(a,0),\forall a\in\mathbb{R}$. Despite $g_{d}$ relying solely on
one or multiple covariates, it functions across all covariates collectively.
The model encompasses $D$ terms in total.

Building a MARS model consists of two main steps, a forward model building
procedure and a backward pruning procedure. In the forward step, functions
$g_{d}({\boldsymbol{x}})$ are chosen successively until a large, often
overfitted model is obtained, which then is pruned to a smaller model. To go
through these steps in detail, let $\mathcal{D}$ be the set of chosen
combinations of basis functions for the model. Initially,
$\mathcal{D}=\\{g_{0}\\}$, with $g_{0}=1$. Next, all functions in
$\mathcal{D}$ are successfully multiplied by the pairs $\mathcal{C}$ and the
term

$\widehat{\beta}_{D+1}g_{r}({\boldsymbol{x}})(x_{j}-c_{j})_{+}+\widehat{\beta}_{D+2}g_{r}({\boldsymbol{x}})(c_{j}-x_{j})_{+},\quad
g_{r}(x)\in\mathcal{D}$

that produces the largest decrease in the training error is added to the
model. This greedy step is repeated until a maximum number of terms is reached
or some other stopping criterion is fulfilled.

After the stopping criterion is met, the model is pruned by removing the term
that leads to the smallest increase in residual square error. This produces an
estimated best model $\widetilde{f}_{\lambda}$, for each number of terms
$\lambda$. For each $\widetilde{f}_{\lambda}$ the generalized cross-validation

$GCV_{f}(\lambda)=\frac{\sum_{i=1}^{n}(f_{i}-\widetilde{f}_{\lambda}({\boldsymbol{x}}_{i}))^{2}}{\left(1-\frac{\mathcal{D}(\lambda)}{n}\right)^{2}},$

is calculated, where $f_{i}$ is the response at data point $i$.
$\mathcal{D}(\lambda)$ is the effective number of parameters as described in
[28]. Finally, the model that minimises $GCV_{h}(\lambda)$ is chosen to be
$\widehat{f}$. In this work, we make use of the MARS implementation from the
R-package “earth” [38].

###  3.2 MCMC Algorithms

####  3.2.1 Physical-informed proposal via preconditioned Crank–Nicolson
algorithm

The primary challenge encountered in Bayesian methods lies in the slow
convergence of posterior distribution sampling, a problem exacerbated as the
dimension increases. To address this, a groundbreaking solution was presented
in [13], introducing a transformative approach to Markov Chain Monte Carlo
(MCMC) method that effectively eliminates the “curse of dimensionality” for
Gaussian priors, given appropriate assumptions on the likelihood. This
innovative technique, known as the “preconditioned Crank–Nicolson” (pCN)
algorithm, exhibits dimension-independent sampling efficiency within its
resulting Markov chain. The pCN algorithm, a modification of standard random-
walk Metropolis sampling, scales the previous iteration’s position before
adding a random shift to generate a new proposal. Notably, the pCN proposal
kernel arises from discretizing an Ornstein–Uhlenbeck process with an
invariant measure, ensuring reversibility with respect to the prior. This
unique feature simplifies the acceptance ratio to a straightforward ratio of
the likelihood.

The physical-informed prior distribution in our framework is a Gaussian
distribution with non-zero mean ${\boldsymbol{\mu}}$ and covariance matrix
${\boldsymbol{\Sigma}}$. Then the pCN proposal distribution, outlined in
Algorithm 1, at step $k$, is as

$q({\boldsymbol{x}}^{\prime}|{\boldsymbol{x}}^{(k)})=\phi\left({\boldsymbol{\mu}}+\sqrt{1-s^{2}}({\boldsymbol{x}}^{(k)}-{\boldsymbol{\mu}}),s^{2}{\boldsymbol{\Sigma}}\right),$
(5)

with $s>0$ being a tuning parameter. This proposal is reversible with respect
to the prior distribution [40]. It means that the acceptance probability
reduces to the likelihood ratio

$A({\boldsymbol{x}}^{\prime},{\boldsymbol{x}}^{(k)})=\min\Bigg{\\{}1,\frac{f({\boldsymbol{y}}|{\boldsymbol{x}}^{\prime})}{f({\boldsymbol{y}}|{\boldsymbol{x}}^{(k)})}\Bigg{\\}}.$

This results in notable acceleration for problems discretized on an
$N$-dimensional grid. Moreover, the pCN method exhibits sampling efficiency
that is independent of dimension, meaning that unlike the standard random
walk, the number of steps required with pCN does not grow as the dimension $N$
increases. Consequently, pCN efficiently explores the posterior support. The
tuning of the parameter $s$ in pCN, as discussed in [13], is aimed at
achieving an acceptance rate of approximately 23.4%, [19].

Algorithm 1 Proposed sampling procedure

  * 1.

Choose the stepsize $s\in(0,1)$. Set $k=0$ and draw ${\boldsymbol{x}}^{(0)}$
from the prior, see (4).

  * 2.

Draw ${\boldsymbol{\xi}}\sim N(\mathbf{0},{\boldsymbol{\Sigma}})$ (using FFT
method).
Then propose the new point as

${\boldsymbol{x}}^{\prime}={\boldsymbol{\mu}}+\sqrt{1-s^{2}}({\boldsymbol{x}}^{(k)}-{\boldsymbol{\mu}})+s{\boldsymbol{\xi}}.$

  * 3.

Set ${\boldsymbol{x}}^{(k+1)}={\boldsymbol{x}}^{\prime}$ with probability
$A({\boldsymbol{x}}^{\prime},{\boldsymbol{x}}^{(k)})$.
Otherwise set ${\boldsymbol{x}}^{(k+1)}={\boldsymbol{x}}^{(k)}$.
Increment $k$ by one and return to step 2.

Utilizing the pCN method within our framework, where the proposed approach
involves incorporating a physically informed prior, enables effective
exploration of various high-probability regions within the intricate posterior
distribution. This aspect holds particular significance in the context of
Bayesian inverse problems, as emphasized in previous works such as [39, 31,
32].

####  3.2.2 Other MCMC samplers

Choosing an effective proposal distribution poses a significant and frequently
challenging undertaking in the application of the Metropolis-Hastings (MH)
algorithm. An optimal proposal distribution should adequately cover the tails
of the posterior distribution while also bearing similarity to the posterior
itself, as outlined by [20]. Additionally, the proposal distribution should
facilitate straightforward sampling. The selection of the proposal
distribution plays a crucial role in influencing the efficiency of the MH
algorithm. In addressing our specific problem, we will explore various
suggestions for proposal distributions, including the random walk, pCN, and
Metropolis adjusted Langevin algorithm (MALA).

Given the physics-informed prior as the Gaussian distribution,
$\mathcal{N}({\boldsymbol{\mu}},{\boldsymbol{\Sigma}})$ , in Section 2.2.3,
with $s>0$ being a tuning parameter, the random walk proposal is as

$q_{1}({\boldsymbol{x}}^{p}|{\boldsymbol{x}}_{t})=\phi({\boldsymbol{x}}_{t},s^{2}{\boldsymbol{I}}),$

random walk with covariance from the prior

$q_{2}({\boldsymbol{x}}^{p}|{\boldsymbol{x}}_{t})=\phi({\boldsymbol{x}}_{t},s^{2}{\boldsymbol{\Sigma}}),$

the preconditioned Crank-Nicolson proposal in (5) denoted by $q_{3}$, and the
MALA

$q_{4}({\boldsymbol{x}}^{p}|{\boldsymbol{x}}_{t})=\phi\left({\boldsymbol{x}}_{t}+\frac{s^{2}}{2}\nabla\log\left(\pi({\boldsymbol{x}}_{t})\right),s^{2}I\right).$

####  3.2.3 Sampling and evaluating large Gaussian random fields via FFT

Let us emphasize that the random variable dimensionality is $N=44,144$. The
proposed methodology involves obtaining high-dimensional Gaussian random
fields from a distribution of this size using the fast Fourier transform (FFT)
algorithm, as detailed in [14]. Similar techniques were previously employed in
expedited 3D elastic inversion, as seen in [9].

The FFT algorithm, embedding the covariance matrix onto a torus, avoids the
need to construct the full prior covariance matrix for the entire domain,
enabling efficient realization generation in the Fourier domain. The general
idea is to use structure of circulant matrices and perform matrix-vector
multiplication in a lower-dimensional space in the Fourier domain. Additional
details can be found in A. In the context of sampling high-dimensional
Gaussian random fields, the approach overcomes memory limitations by
leveraging circulant matrix structures and Fourier domain transformations, as
documented in [10, 45, 14, 1]. This method facilitates rapid sample generation
and proves computationally advantageous for estimating log-density in
scenarios requiring the evaluation of probability density functions, such as
calculating acceptance rates in MCMC algorithms.

##  4 Numerical study: small-scale data

In this section, in Subsection 4.1, we initially conduct numerical
investigations to assess the efficacy of MARS in approximating the true
forward function. Subsequently, in Subsection 4.2, we scrutinize comparisons
between MCMC samples derived from the true model and those from the
approximated model within a limited zone of the Alvheim field. Finally, in
Subsection 4.3, we assess the efficiency of various proposal distributions.

###  4.1 Evaluating approximation of the forward model

To approximate the forward model in the Alvheim case, a training data set is
created. The data set has $20,000$ observations of covariates $x_{g}$,
$x_{o}$, $x_{c}$ and $x_{d}$. The value of data point $i$ at covariate $j$ is
denoted by $x_{ij}$, where $j=g,o,c,d$ corresponds to gas, oil, clay and depth
respectively. Data point $i$ is sampled as follows: First, a value for
$x_{id}$ is sampled from a uniform distribution on the interval between the
minimum and maximum depth at the top reservoir of the Alvheim field shown in
Figure 2(a). Next $x_{id}$ is used to find the expected value of the prior and
draw $x_{ig}$, $x_{io}$ and $x_{ic}$ from the prior distribution. The test
data are created in the same manner as the training data. It has $44144$ test
data points, which is the size of one sample ${\boldsymbol{x}}$ at the Alvheim
field.

Recall that there are two types of seismic AVO data, the zero offset ($R_{0}$)
and the seismic gradient ($G$). That is,
${\boldsymbol{h}}:\mathbb{R}^{4}\to\mathbb{R}^{2}$. To approximate
${\boldsymbol{h}}$, two models are created. A model for the zero offset,
$\widehat{h_{R_{0}}}$, and one for the seismic gradient $\widehat{h_{G}}$ such
that $\widehat{{\boldsymbol{h}}}$ collectively denotes the two models
$\widehat{h_{R_{0}}}$ and $\widehat{h_{G}}$ with
$\widehat{h_{R_{0}}}:\mathbb{R}^{4}\to\mathbb{R}$ and
$\widehat{h_{G}}:\mathbb{R}^{4}\to\mathbb{R}$.

As candidates for approximating ${\boldsymbol{h}}$, the MARS models
$\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ and
$\widehat{{\boldsymbol{h}}}_{\text{MARS}_{\text{OF}}}$ are trained on the
$20,000$ data points. The large MARS models for the two AVO properties
$h_{R_{0}}$ and $h_{G}$ have $39$ and $40$ terms of which $26$ and $24$ are
interaction terms respectively. The pruned MARS model for $h_{R_{0}}$ has $25$
terms of which $17$ are interaction terms, while the pruned MARS model for
$h_{G}$ has $25$ of which $13$ are interaction terms. In addition the kernel
regression models $\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{1000}}$ and
$\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{4000}}$ are trained, as described in
B, with $1000$ and $4000$ of the training data respectively. They are trained
on subsets of the training data set because of computation time. As the NPKR
models are weighted averages of the training data, the number of training data
affects the prediction time. The subsets of the training data are random
subsets, however one could try to select these points in a more optimal way to
improve the NPKR models as discussed in B.

The accuracy of the predictions are examined by evaluating the sample
correlation

$\text{Corr}\left(f,\widehat{f}\right)=\frac{\sum_{i=1}^{n}\left(f-\bar{f}\right)\left(\widehat{f}_{i}-\bar{\widehat{f}}\right)}{\sqrt{\sum_{i=1}^{n}\left(f_{i}-\bar{f}\right)^{2}}\sqrt{\sum_{i=1}^{n}\left(\widehat{f}_{i}-\bar{\widehat{f}}\right)^{2}}}.$
(6)

and MSE

$\text{MSE}\left(f,\widehat{f}\right)=\frac{1}{n}\sum_{i=1}^{n}\left(f_{i}-\widehat{f}_{i}\right)^{2}.$
(7)

Precision metrics, such as correlation and MSE are the mean of these metrics
for the two models. That is, the correlation between a model
$\widehat{{\boldsymbol{h}}}$ and the response ${\boldsymbol{h}}$ is the mean
of $\text{Corr}\left(h_{R_{0}},\widehat{h_{R_{0}}}\right)$ and
$\text{Corr}\left(h_{G},\widehat{h_{G}}\right)$ calculated using equation (6).
Similarly, the MSE is the mean of
$\text{MSE}\left(h_{R_{0}},\widehat{h_{R_{0}}}\right)$ and
$\text{MSE}\left(h_{G},\widehat{h_{G}}\right)$ from equation (7). Moving on
$\widehat{{\boldsymbol{h}}}_{\text{MARS}}$,
$\widehat{{\boldsymbol{h}}}_{\text{MARS}_{\text{OF}}}$,
$\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{1000}}$ and
$\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{4000}}$ are referred to as one
model, although they are two models, as described in the second paragraph of
this section.

The computation time for predicting the test data, correlation and MSE are
reported in Table 1. The computation time is the average of predicting the
test data 50 times. For comparison, the computation time of the true forward
model used on average $2.059$ seconds to compute
${\boldsymbol{h}}({\boldsymbol{x}})$. The fastest model was
$\widehat{{\boldsymbol{h}}}_{\text{MARS}}$, which was on average $34$ times
proximately faster than the true forward model. The larger MARS model was also
very fast, although slower than the smaller MARS model. The larger MARS model
is marginally better than $\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ when
considering the correlation and MSE, however, it is much more complex. By
setting the maximum number of terms in the MARS model after pruning to a lower
number, the computation time is reduced. This is at a cost of increasing MSE
and decreasing correlation.

The correlation achieved by $\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{4000}}$
approaches that of the MARS models, the MARS models demonstrate an MSE that is
half that of $\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{4000}}$, along with
computation times that are, on average, over 200 times faster. Notably,
$\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{1000}}$ displays the poorest
performance in terms of both MSE and correlation. The results in Table 1
reveal significantly slower computational speeds for the NPKR (kernel
regression) models compared to the MARS models, including the forward function
${\boldsymbol{h}}$. Moreover, the computational time increases with the
augmentation of data points in the NPKR model, evident in the slower
computation of $\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{4000}}$ compared to
$\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{1000}}$. Beyond the reported
computation times in Table 1, an additional test involved predicting test data
using an NPKR model trained with 100 data points, which averaged 0.674 seconds
for prediction. Although this time is considerably shorter than that of
$\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{4000}}$ and
$\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{1000}}$, it still remains
significantly longer than the computation times observed for the MARS models.

$\widehat{{\boldsymbol{h}}}$ | correlation | MSE$\times 10^{-5}$ | computation time [sec]
---|---|---|---
$\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ | 0.995 | 3.8 | 0.06
$\widehat{{\boldsymbol{h}}}_{\text{MARS}_{\text{OF}}}$ | 0.996 | 3.4 | 0.11
$\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{1000}}$ | 0.983 | 14.0 | 5.69
$\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{4000}}$ | 0.991 | 7.2 | 22.7
Table 1: Correlation, MSE and computation time for
$\widehat{{\boldsymbol{h}}}_{\text{MARS}}$,
$\widehat{{\boldsymbol{h}}}_{\text{MARS}_{\text{OF}}}$,
$\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{1000}}$ and
$\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{4000}}$.

Figure 4 shows a scatter plot of the ${\boldsymbol{h}}$ versus
$\widehat{{\boldsymbol{h}}}_{\text{MARS}}$,
$\widehat{{\boldsymbol{h}}}_{\text{MARS}_{\text{OF}}}$ on the test data.
Overall the figures indicate that the errors are slightly smaller for the
predictions of the seismic gradient than for the zero offset. Neither
$\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ nor
$\widehat{{\boldsymbol{h}}}_{\text{MARS}_{\text{OF}}}$ seems to make grave
systematic errors, except perhaps for the lowest values of $h_{R_{0}}$. The
predictions in Figures 4(a) and 4(c) look almost identical. This is confirmed
by Figure 5, where the predictions of
$\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ are plotted against the predictions
of $\widehat{{\boldsymbol{h}}}_{\text{MARS}_{\text{OF}}}$. The $14$ more terms
in model $\widehat{{\boldsymbol{h}}}_{\text{MARS}_{\text{OF}}}$ seems
excessive when predictions are close to identical.

(a) True forward model values for $h_{R_{0}}$ plotted against
$\widehat{{\boldsymbol{h}}}_{\text{MARS}}$.

(b) True forward model values for $h_{G}$ plotted against
$\widehat{{\boldsymbol{h}}}_{\text{MARS}}$.

(c) True forward model values for $h_{R_{0}}$ plotted against
$\widehat{{\boldsymbol{h}}}_{\text{MARS}_{\text{OF}}}$.

(d) True forward model values for $h_{G}$ plotted against
$\widehat{{\boldsymbol{h}}}_{\text{MARS}_{\text{OF}}}$.

Figure 4: True forward model plotted as a function of predicted forward models
on test data. Figure 5: Predictions of $h_{R_{0}}$ from
$\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ plotted as a function of
predictions of $h_{R_{0}}$ from
$\widehat{{\boldsymbol{h}}}_{\text{MARS}_{\text{OF}}}$. The black line is
$\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ =
$\widehat{{\boldsymbol{h}}}_{\text{MARS}_{\text{OF}}}$.

Section 5 gives results from MCMC simulations at the Alvheim field using
$\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ as a replacement of
${\boldsymbol{h}}$. This model was chosen because it was the fastest, and had
marginally higher MSE and lower correlation than
$\widehat{{\boldsymbol{h}}}_{\text{MARS}_{\text{OF}}}$. Considering
correlation and MSE, both MARS models would serve as good substitutes for the
forward model, while reducing the computation time considerably.

In pursuit of computational efficiency, the NPKR models investigated in this
study do not serve as optimal replacements for the forward model. An inherent
drawback of the NPKR model is its diminished performance with an increased
volume of training data, thereby decelerating model execution. The reduction
in training data quantity directly correlates with decreased computation time
for predictions. For instance, conducting predictions on 50 test data points
utilizing an NPKR model trained on 100 data points resulted in an average
duration of 0.674 seconds, which is swifter compared to the forward model
(${\boldsymbol{h}}$). Nonetheless, this pace remains significantly slower than
that achieved by MARS models, and it comes with the trade-off of elevated MSE
and diminished correlation. Implementing a more strategic data sampling
methodology may enhance correlation and mitigate MSE; however, constructing a
model based on four thousand data points exhibited a correlation proximate to
MARS models while incurring substantially higher MSE.

###  4.2 Comparison of the surrogate and the exact forward models

In Section 4.1, the analysis demonstrates that utilizing the MARS surrogate
resulted in a speed enhancement of approximately $32$ times when employing the
true forward function. Furthermore, the MARS surrogate exhibited a low MSE and
a strong correlation. However, given its nature as an approximation, the MCMC
samples are drawn from a posterior distribution which serves as an
approximation to the true posterior distribution, denoted as $\pi$. For
clarity, the approximated posterior, where ${\boldsymbol{h}}$ is substituted
by $\widehat{{\boldsymbol{h}}}_{\text{MARS}}$, is labeled as $\widehat{\pi}$.
To assess the resemblance between $\widehat{\pi}$ and $\pi$, MCMC samples are
generated using both ${\boldsymbol{h}}$ and
$\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ as forward models. Additionally,
MCMC samples employing $\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ as the
forward model and $\widetilde{{\boldsymbol{\Omega}}}$, as defined in equation
(11), as the covariance matrix of the likelihood model, are examined to
determine whether samples from this posterior are more akin to those from
$\pi$ than those from $\widehat{\pi}$. The posterior with the forward function
$\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ and the covariance matrix
$\widetilde{{\boldsymbol{\Omega}}}$ in the likelihood model is denoted as
$\widetilde{\pi}$. The details regarding the approximation
$\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ and the covariance matrix
$\widetilde{{\boldsymbol{\Omega}}}$ are provided in B.2.

Figure 6: The black square is a smaller area where MCMC with the true forward
function is performed.

The MH algorithm with proposal distribution $q_{2}$ is used to sample the
posteriors $\pi$, $\widehat{\pi}$ and $\widetilde{\pi}$. The tuning parameter
$s$ is tuned such that the acceptance rate is approximately $23.4\%$, [19].
The value $s=0.027$ was used, which gave the acceptance rate $23.5\%$ when
using the true forward model, $23.4\%$ when using
$\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ and $23.5\%$ when using
$\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ and
$\widetilde{{\boldsymbol{\Omega}}}$ as covariance matrix of the likelihood.
The initial sample was the mean of the prior distribution. The number of
iterations was $10^{6}$ where every $10$th sample is saved. The initial
portion of the MCMC samples is discarded as burn-in, ensuring that the
subsequent $50,000$ samples are presumed to be drawn from the stationary
distributions of the Markov chains.

Due to computation time, this was only done for a smaller area located at
inlines$\times$crosslines $[928,1160]\times[4788,5020]$, where every other
inline and crossline is omitted, such that the location of the measurements
are approximately $100$ meters apart. This means that the smaller area is
discretized by $30\times 30=900$ grid points. Hence
${\boldsymbol{x}}_{g},{\boldsymbol{x}}_{o},{\boldsymbol{x}}_{c}\in\mathbb{R}^{900}$.
Figure 6 shows the smaller area at the Alvheim field where the MCMC was
performed. In the vicinity of the well, the samples are subject to the
influence of the likelihood model for the well-log data. Therefore, it is of
particular interest to scrutinize the similarities and dissimilarities, with a
focus on locations farther away from the well.

(a) Mean gas saturation from MCMC samples from $\pi$.

(b) Mean gas saturation from MCMC samples from $\widehat{\pi}$.

(c) Mean gas saturation from MCMC samples from $\widetilde{\pi}$.

(d) Uncertainty of gas saturation from MCMC samples from $\pi$.

(e) Uncertainty of gas saturation from MCMC samples from $\widehat{\pi}$.

(f) Uncertainty of gas saturation from MCMC samples from $\widetilde{\pi}$.

Figure 7: Mean gas saturations and uncertainties from MCMC samples from three
posteriors. The posterior where $\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ is
surrogate is $\widehat{\pi}$. In the last posterior, $\widetilde{\pi}$, both
the surrogate $\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ and the error
correcting covariance matrix $\widetilde{{\boldsymbol{\Omega}}}$ are used. The
MH algorithm with proposal distribution $q_{2}$ with $s=0.027$ was used in all
three cases. Uncertainty is expressed as the difference between the 90th and
10th quantiles.

(g) Mean oil saturation from MCMC samples from $\pi$.

(h) Mean oil saturation from MCMC samples from $\widehat{\pi}$.

(i) Mean oil saturation from MCMC samples from $\widetilde{\pi}$.

(j) Uncertainty of oil saturation from MCMC samples from $\pi$.

(k) Uncertainty of oil saturation from MCMC samples from $\widehat{\pi}$.

(l) Uncertainty of oil saturation from MCMC samples from $\widetilde{\pi}$.

Figure 8: Mean oil saturations and uncertainties from MCMC samples from three
posteriors. The posterior where $\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ is
surrogate is $\widehat{\pi}$. In the last posterior, $\widetilde{\pi}$, both
the surrogate $\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ and the error
correcting covariance matrix $\widetilde{{\boldsymbol{\Omega}}}$ are used. The
MH algorithm with proposal distribution $q_{2}$ with $s=0.027$ was used in all
three cases. Uncertainty is expressed as the difference between the 90th and
10th quantiles.

(a) Mean brine saturation from MCMC samples from $\pi$.

(b) Mean brine saturation from MCMC samples from $\widehat{\pi}$.

(c) Mean brine saturation from MCMC samples from $\widetilde{\pi}$.

(d) Uncertainty of brine saturation from MCMC samples from $\pi$.

(e) Uncertainty of brine saturation from MCMC samples from $\widehat{\pi}$.

(f) Uncertainty of brine saturation from MCMC samples from $\widetilde{\pi}$.

Figure 9: Mean brine saturations and uncertainties from MCMC samples from
three posteriors. The posterior where
$\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ is surrogate is $\widehat{\pi}$. In
the last posterior, $\widetilde{\pi}$, both the surrogate
$\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ and the error correcting covariance
matrix $\widetilde{{\boldsymbol{\Omega}}}$ are used. The MH algorithm with
proposal distribution $q_{2}$ with $s=0.027$ was used in all three cases.
Uncertainty is expressed as the difference between the 90th and 10th
quantiles.

(a) Mean clay content from MCMC samples from $\pi$.

(b) Mean clay content from MCMC samples from $\widehat{\pi}$.

(c) Mean clay content from MCMC samples from $\widetilde{\pi}$.

(d) Uncertainty of clay content from MCMC samples from $\pi$.

(e) Uncertainty of clay content from MCMC samples from $\widehat{\pi}$.

(f) Uncertainty of clay content from MCMC samples from $\widetilde{\pi}$.

Figure 10: Mean clay content and uncertainties from MCMC samples from three
posteriors. The posterior where $\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ is
surrogate is $\widehat{\pi}$. In the last posterior, $\widetilde{\pi}$, both
the surrogate $\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ and the error
correcting covariance matrix $\widetilde{{\boldsymbol{\Omega}}}$ are used. The
MH algorithm with proposal distribution $q_{2}$ with $s=0.027$ was used in all
three cases. Uncertainty is expressed as the difference between the 90th and
10th quantiles.

Figure 8 illustrates mean gas saturations and corresponding uncertainties from
three sets of posterior samples, which exhibit remarkable similarity overall.
A notable distinction is the wider light green area at approximately inline
$1050$ and crossline $4955$ in Figure 7(a) compared to Figures 7(b) and 7(c),
with a corresponding increase in uncertainty. Additionally, at inline $1000$
and crossline $4950$, uncertainties in Figures 7(d) and 7(e) align more
closely, while Figure 7(f) displays lower uncertainty in this region. In the
yellow area around inline $1100$ and crossline $4920$, Figures 7(b) and 7(c)
are similar, with Figure 7(a) showing slightly less intensity. Similarly,
differences in uncertainty patterns are observed, especially in the bottom
right corner, where samples using the surrogate model are more alike and
differ slightly from those using the exact forward model.

In Figure 8, mean oil saturations and uncertainties demonstrate consistent
patterns across all cases. However, discrepancies arise, particularly in the
lower left corner, where samples from $\pi$ indicate lower oil saturation
compared to $\widehat{\pi}$ and $\widetilde{\pi}$. The use of the error-
correcting covariance matrix improves the approximation to $\pi$, as evident
in Figure 7(i) compared to Figure 7(h). Conversely, Figure 7(i) exhibits a
small green area in the lower right corner not present in Figures 7(h) and
7(g). Differences in uncertainty are also notable, with higher uncertainty in
Figure 7(l) compared to Figures 7(j) and 7(k). Areas where uncertainties in
samples from $\widehat{\pi}$ and $\widetilde{\pi}$ resemble each other but
differ from $\pi$ are observed around inline $1060$ and crossline $4900$, with
lower uncertainty in samples from $\pi$. Similarly, at approximately inline
$1060$ and crossline $5020$, oil saturation samples from $\widehat{\pi}$ and
$\widetilde{\pi}$ exhibit greater similarity compared to $\pi$.

The mean brine saturation and corresponding uncertainties illustrated in
Figure 9 exhibit consistency across the three cases. However, a disparity
emerges around inline $1060$ and crossline $5020$, where Figures 9(b) and 9(c)
indicate higher brine saturations compared to Figure 9(a). This corresponds to
the area where differences in oil saturation were observed. The most
significant differences in brine saturation uncertainties occur in the bottom
left and right corners, where Figures 9(c) and 9(d) exhibit stronger
similarities compared to Figures 9(e) and 9(d).

Regarding the mean clay content depicted in Figures 10, they are largely
similar, with slight differences observed towards crossline $4955$. Notably,
the mean clay content in Figure 10(a) is notably lower in this area compared
to Figures 10(b) and 10(c). Samples from $\widehat{\pi}$ in this region
slightly resemble those from $\pi$ compared to their overall similarity.
Uncertainty patterns in Figures 10 are generally similar. However,
uncertainties may be slightly larger in Figure 10(d), particularly around
inline $900\sim 1000$ at higher crosslines, and around inline $1130$ and
crossline $4930$.

Figures 8, 8, 9 and 10 show overall great similarity in the results from the
MCMC samples from $\pi$, $\widehat{\pi}$ and $\widetilde{\pi}$ on the smaller
area in the Alvheim field. This indicates that substituting ${\boldsymbol{h}}$
with $\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ gives a good approximation to
the posterior $\pi$, both when using the covariance matrix
${\boldsymbol{\Omega}}$ and when using $\widetilde{{\boldsymbol{\Omega}}}$ as
covariance matrix in the likelihood model. MARS models were also considered as
good surrogates for the time-consuming forward models in [12] and [11]. A MARS
model does not assume a functional form of the response, which could be the
reason it applied to all three problems.

The use of the adjusted correlation function
$\widetilde{{\boldsymbol{\Omega}}}$ showed minor improvements in some areas,
while in other areas, using ${\boldsymbol{\Omega}}$ gave the means and
uncertainties closest to the mean and uncertainties of the MCMC samples of
$\pi$. The empirical variance and covariance added to the covariance matrix to
adjust for the error between $\widehat{{\boldsymbol{h}}}$ and
${\boldsymbol{h}}$ were very small, such that
$\widetilde{{\boldsymbol{\Omega}}}$ and ${\boldsymbol{\Omega}}$ are quite
similar. That the covariance matrices are very similar could be the reason for
the similar result.

Most of the differences between the exact and approximate samples are found at
the edges of the area. As mentioned, there is a well in that area, such that
samples close to the well are affected by the likelihood model for the well-
log data in addition to the likelihood model for the seismic AVO data.
However, the mean and uncertainty of the saturations are in general very
similar, indicating that $\widehat{\pi}$ and $\widetilde{\pi}$ are similar to
$\pi$. In MCMC, there is a randomness when proposing and accepting samples,
such that the exact same samples might not be considered for the three
different Markov chains. This could lead to some differences in the means and
uncertainties, however, this difference vanishes with an increasing number of
MCMC iterations if $\pi$, $\widehat{\pi}$ and $\widetilde{\pi}$ were the same
distributions.

###  4.3 Comparing the proposals distributions

Next, the performance of the MH algorithm with the four different proposal
distributions, $q_{1}$, $q_{2}$, $q_{3}$ and $q_{4}$ given in Section 3.2.2,
are compared, where the tuning parameter is denoted by $s$. The approximation
$\text{MARS}_{\nabla}$ is used instead of the gradient of the log posterior in
the MALA to reduce computation time.

Here, we focus the approximate posterior $\widehat{\pi}$ on the small area.
That is, $\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ is used as the forward
model and the covariance matrix of the likelihood model is
${\boldsymbol{\Omega}}$. The objective is to compare how efficiently the four
MH algorithms explore the domain of $\widehat{\pi}$. The ESS per computation
time is measured to assess which proposal gives the most information about
$\widehat{\pi}$ per time. The reported ESS is the mean ESS of all $3N$
stochastic processes
$x_{g}^{1},x_{o}^{1},x_{c}^{1},x_{g}^{2},x_{o}^{2},x_{c}^{2},...,$
$x_{g}^{N},x_{o}^{N},x_{c}^{N}$. The tuning parameter $s$ and acceptance
rates, which are connected to the efficiencies of the algorithms are reported.
The four MH algorithms were tuned such that the acceptance rates were
approximately $23.4\%$ when using the proposal distributions $q_{1}$, $q_{2}$
and $q_{3}$ and approximately $57.4\%$ when using the proposal distribution
$q_{4}$. As the acceptance rates and ESS should be computed on MCMC samples
after burn-in, the Markov chain starts at the mean of the posterior samples
whose logistic transformation is shown in Figure 7(a), Figure 7(g) and Figure
10(a). The MH algorithms are run for $10^{6}$ iterations where every $100$th
sample is saved. This gives $10,000$ MCMC samples from $\widehat{\pi}$.

proposal | $s$ | acceptance rate [%] | computation time [sec] | ESS [#] | ESS/time
---|---|---|---|---|---
$q_{1}$ | $0.006$ | $25.5$ | $13546.24$ | $11.01$ | $0.0008$
$q_{2}$ | $0.027$ | $23.3$ | $16547.33$ | $93.88$ | $0.0057$
$q_{3}$ | $0.045$ | $22.8$ | $11408.03$ | $212.64$ | $0.0186$
$q_{4}$ | $0.038$ | $58.9$ | $44418.04$ | $39.72$ | $0.0009$
Table 2: Comparison of the efficiency of four MH algorithms. The acceptance
rate, computation time and ESS are calculated on $10,000$ MCMC samples.

These results, reported in Table 2, show that the MH algorithm with the
proposal distribution $q_{3}$ was the most efficient. The algorithm produced
approximately $0.0186$ independent samples per second. This is much higher
than the ESS per computation time for the MH algorithm with the other
proposals distributions. The second most efficient algorithm was the MH
algorithm with the proposal distribution $q_{2}$. The ESS of the MH algorithm
with proposal distribution $q_{3}$ is twice better than the ESS for $q_{2}$ as
proposal distribution. The MH algorithm with proposal distribution $q_{3}$ was
also about 1 hour and 25 minutes faster than the second most efficient
algorithm. The slowest algorithm was the MALA, which used more than 12 hours.
For comparison, the fastest algorithm used 3 hours and 10 minutes. The MALA
used more than three times as much time as $q_{1}$, but because the samples
were less correlated, $q_{4}$ was slightly more efficient than $q_{1}$ in this
case.

Using the MH algorithm with $q_{3}$ results in faster MCMC sampling compared
to the other three distributions. It is sensible that employing the $q_{4}$
distribution leads to the slowest sampling, as it lacks symmetry like $q_{1}$
and $q_{2}$, and it is not reversible concerning the prior, unlike $q_{3}$,
necessitating gradient evaluations. Similarly, $q_{2}$ is reasonably the
second slowest, given its proposal correlation akin to $q_{3}$ and the
requirement to assess both the prior and likelihood similar to $q_{1}$. While
using $q_{1}$ does not require proposing correlated samples as $q_{3}$ does,
it does entail evaluating the prior, unlike $q_{3}$, where only the likelihood
needs evaluation. The computation time in Table 2 illustrates that proposing
correlated samples is faster than evaluating the proposal distribution. In
this context, the MALA method notably outperforms the MH algorithm with
$q_{1}$ as the proposal distribution, suggesting the significance of utilizing
posterior information as outlined in previous studies. Although MALA’s
computation time is over three times that of the MH algorithm with $q_{1}$,
its reduced sample correlation enables it to convey slightly more information
about $\widehat{\pi}$ per second.

In seismic AVO inversion within the Alvheim field, testing four MH algorithms
revealed that MH with proposal distribution $q_{3}$ was the most efficient in
terms of ESS per computation time. It exhibited both the highest ESS and the
smallest computation time among the four algorithms, making it the clear
choice for efficiency. However, this does not guarantee its efficiency in
other contexts, as finding suitable MCMC methods for various inverse problems
can be challenging and problem-specific. Despite this, $q_{3}$ has proven
effective beyond the Alvheim case, outperforming other methods, like using
$q_{1}$ [44], particularly with increased dimensionality. This efficiency,
irrespective of dimensionality, makes $q_{3}$ preferable in high-dimensional
scenarios like the Alvheim case.

(a) Mean saturation of gas.

(b) Mean saturation of oil.

(c) Mean saturation of brine.

(d) Mean clay content.

(e) Uncertainty of gas saturation.

(f) Uncertainty of oil saturation.

(g) Uncertainty of gas saturation.

(h) Uncertainty of clay content.

Figure 11: Mean and uncertainty of gas, oil, brine and clay saturations from
MCMC samples from $\widehat{\pi}$. The uncertainty is the difference in the
$90$th and $10$th quantile.

##  5 Results with the Alvheim field data

As highlighted from simulation studies in Section 4, preconditioned Crank-
Nicolson algorithm show the best performance in our comparison thus we focus
on this algorithm for our case study on the Alvheim field data. Additional
results for other proposal are given in C.

In this section, MCMC samples from $\widehat{\pi}$ over the entire Alvheim
field are compared to the results of [46]. The two results are compared by
mean saturation of gas, oil and brine, mean clay content, uncertainties and
ternary plots from one deep and one shallow location in the Alvheim field. The
MH algorithm with $q_{3}$ as proposal distribution is used to draw $500,000$
MCMC samples after burn-in where every $50$th sample is saved. That gives
$10,000$ MCMC samples from the posterior $\widehat{\pi}$. The tuning parameter
was $s=0.0088$ and the acceptance rate was $25.4\%$. The Markov chain is
started in two locations to check that the Markov chain converges to the same
values.

The MCMC results are shown in Figure 11. The wells are marked by coloured
circles: red circles indicate wells where mainly gas was found and the violet
circle indicates the well where the primarily oil was found. Figure 11 shows
that near the wells where mainly gas was found, the mean gas saturation is
high and the mean oil saturation is low. Likewise, the mean gas saturation is
low and the mean oil saturation is high near the well where most oil was
retrieved. The uncertainty is generally low near the four wells. The mean gas
and oil saturations are generally low at the deeper areas, such as in the
bottom left corner, the top left corner, in the top middle and between the gas
and oil wells in the area that stretches from approximately (inline,crossline)
= (600,5200) to (1300,4500). The uncertainty of the saturations is generally
low, however higher for the oil saturations than the gas saturations. The
uncertainties are low in the middle of areas which show high oil or gas
saturations, and higher where the saturations transition from high to low
values, such that the uncertainties look like circles enclosing the areas with
high gas or oil saturations.

### Some similarities to Spremić et al., [46]

There are several similarities between our results and those from [46]. The
results in [46, Figure 12] in show high mean saturation of gas and oil near
the wells where mostly gas and oil were found respectively. There is also low
uncertainty close to the four wells. In general, the areas of high mean gas
saturation look similar to each other, but the mean saturation of gas from the
MCMC samples has more defined shapes compared to the mean gas saturation in
[46, Figure 12], which is more square-shaped and blurred around the edges of
high mean saturation areas. Another similarity is the low mean saturation of
gas and oil between the gas and oil wells in the bottom left area, the top
left area and the area that stretches from approximately (inline,crossline) =
(600,5200) to (1300,4500). The uncertainties of the clay content in [46,
Figure 13] and Figure 11 also look very similar.

The mean oil saturations in the two results show some dissimilarities. For
example, in [46, Figure 12] the well marked by the green triangle, there are
modest tendencies toward high oil saturations, whereas in the same area in
Figure 11(b), mean oil saturation values are close to one. There is also an
area to the right of the well marked by the violet circle in Figure 11(b)
which shows high mean oil saturation. This area shows very low mean oil
saturation in [46, Figure 12]. Another example of areas where the MCMC results
show higher oil saturation compared to the results of [46] are two areas below
the well located the furthest to the left at approximately inline $500$ and
crossline $4800$. Conversely above this well at approximately crossline
$4600\sim 4700$, there is an area where the mean oil saturation from the MCMC
results is very low and the mean oil saturation in [46, Figure 12] is high.

### Some differencies to Spremić et al., [46]

Another significant difference is the difference in the uncertainties for the
mean gas and oil saturations. In [46, Figure 12] the uncertainty is high for
the mean oil saturations in the areas described in the previous paragraph.
That is the area below the well marked by the green triangle, the area below
the leftmost well and the area over the leftmost well. In addition, the
uncertainty of the mean oil saturation is high above the well between the
leftmost well and the well marked by the green triangle. The uncertainty of
the gas saturations in [46, Figure 12] are high in the same areas as the mean
gas saturation is high, except near the wells. The uncertainty of the mean gas
saturation from the MCMC results in Figure 11(h) are generally considerably
lower.

There is also a difference in the mean clay content, which shows higher values
for the MCMC results compared to the results of [46]. The mean clay content is
high in the same areas, however, much higher in Figure 11(h),11(d) compared to
the mean clay content in [46, Figure 13].

Ternary plots at two locations in the Alvheim field for the two sets of
samples from the two approximation methods are shown in Figure 12 with
$10,000$ values of MCMC samples from $\widehat{\pi}$ at locations $(800,4800)$
and $(436,5432)$ respectively. Ternary plots for the $100$ approximate
posterior samples from the ensemble based approach at location $(800,4800)$
and $(436,5432)$ are shown in Figure 12(b) and LABEL:fig:ternary_deep_Mina,
respectively. The depth at the location $(800,4800)$ is $2112$ meters and the
depth at the location $(436,5432)$ is $2190$ meters.

The ternary plots in Figure 12 show that the obtained samples in the two
methods give the same composition at the deep location. In the shallow
location, however, there are some differences. The MCMC results show less
uncertainty in the oil and gas saturations. The oil saturations were also
higher in Figure 12(a) compared to Figure 12(b). This agrees with the higher
mean oil saturations in the MCMC results shown in Figure 11(b) compared to the
mean oil saturations in [46, Figure 12]. The uncertainty of the brine
saturations in Figure 12(a) and Figure 12(b) is about the same. The ternary
plots are only from two locations, such that general conclusions cannot be
drawn. However, they indicate that the MCMC samples had more oil and lower
uncertainty in the shallow area compared to the results of [46], which agrees
with results in Figure 11 and [46, Figure 12].

Even though there are dissimilarities between the results from the two
different methods for approximating the posterior $\pi$, the results show the
same tendencies. For most of the areas where the MCMC result shows a high mean
oil saturation and the ensemble-based approach shows a higher mean oil
saturation, the ensemble-based approach also has a higher uncertainty of the
mean oil saturation. It means that the two results are compatible with each
other. There is an exception around inline $1250$ and crossline $4900$ where
the MCMC results show a high mean oil saturation and the ensemble-based
approximation to the posterior shows a low mean oil saturation and low
uncertainty. A possible reason could be that this area is located at the edge
of the area such that there is no AVO data to the right of this area in the
conditioning in the ensemble-based approximation.

(a) Ternary plots for the MCMC samples from $\widehat{\pi}$.

(b) Ternary plots for the approximate posterior samples from [46].

Figure 12: Ternary plots for samples from the two approximations to the
posterior $\pi$ at location $(800,4800)$ in the Alvheim field, where the depth
is $2112$ meters.

In the approach described in [46], the prior samples from the total area are
divided into smaller patches when conditioning on the AVO data, to prevent
false correlations between the AVO data which can occur when using such
methods [46]. Two patch sizes were compared to each other, and when using
larger patches, there was more oil in the results, which resembles the MCMC
results more.

[46] added an iterative loop of the approximate posterior samples to improve
the posterior approximation. The results indicated a possible increase in the
integrated oil saturation, which resembles the tendencies in the results from
MCMC.

Another explanation could be that $\widehat{\pi}$ has less uncertainty than
$\pi$. Figure 10 showed that the approximate posterior samples had slightly
lower uncertainty of the clay content. However, the uncertainty of the clay
content is one of the similarities between the results from the MCMC and the
method described in [46]. The uncertainties of the gas and oil saturations in
Figure 8 and in Figure 8 did not indicate that using the surrogate leads to a
lower uncertainty. If it is the case that $\widehat{\pi}$ has less uncertainty
than $\pi$ it is not clear from the results in Section 4.2.

##  6 Discussion and conclusion

In this paper, we tackle some challenges in Bayesian reservoir inversion
methods including computational complexities due to forward model evaluations
and high-dimensional Gaussian random fields. Innovative statistical approaches
based on the generalized Bayesian method and efficient techniques such as the
fast Fourier transform (FFT) for simulating Gaussian random fields are
employed, making large-scale computations manageable even on standard
hardware. Novel sampling methodologies using the preconditioned Crank-Nicolson
method offer efficient exploration of high-dimensional parameter spaces,
promising substantial computational speed enhancements in MCMC sampling and
streamlining Bayesian reservoir inversion methods for practical applications.

The forward function of the likelihood model in the Alvheim case was
approximated by a MARS model to speed up the MCMC algorithm. The approximation
was about $32$ times faster. It was also confirmed by a comparison between
MCMC samples obtained using the exact forward model. Using the MARS model, the
MH algorithms with different proposal distributions (random walk,
preconditional Crank-Nicolson, MALA) were compared on a smaller part of the
Alvheim field. A comparison of ESS per computation time showed that the
preconditional Crank-Nicolson algorithm was the most efficient algorithm in
the Alvheim case. This algorithm was therefore used with the surrogate to
perform MCMC on the entire Alvheim field. The MCMC results had similarities to
the results in [46], however, the MCMC results showed more oil and lower
uncertainty for the gas and oil saturations.

Another method of employing the MARS approximation involves utilizing samples
derived from $\widehat{\pi}$ to construct an informed independent proposal.
Research by [31] discovered that devising an independent proposal distribution
based on posterior information enhanced efficiency compared to employing a
dependent standard random walk proposal distribution. Concurrently, an
independent proposal could be executed in parallel to expedite computations.
We have performed a study to concurrently sample from $\pi$ and
$\widehat{\pi}$. That is the Metropolis-Hastings (MH) algorithm with proposal
distribution $q_{3}$ was utilized to generate 500,000 MCMC samples from
$\widehat{\pi}$. Moreover, every 1000th iteration involved subjecting the
proposed sample to an acceptance/rejection step utilizing the acceptance
probability computed using the exact forward model rather than the forward
model approximation. Unfortunately, the acceptance rate was zero for exact
MCMC sampling. A challenge arising from combining dependent and independent
proposals lies in their differing optimal tuning. An independent proposal
demonstrates greater efficiency when the acceptance rate nears unity. However,
for the MH algorithm with proposal $q_{3}$ tuned to attain an acceptance rate
of around 23.4%, the proposed samples may be excessively extreme or distant
from the posterior, potentially leading to rejection by the independent
proposal. Consequently, this approach might be more suitable for the MALA,
which boasts a higher optimal acceptance rate.

Since the MARS model is non-parametric, it can in principle model functions of
arbitrary shape. This suggests that MARS can approximate other forward models
in other problems as well. The forward model in the Alvheim example takes four
continuous input parameters. If there are more covariates, one probably need
more terms in the model, which would increase the computation time a bit.
However, the input space of the MARS model which served as a surrogate in [12]
had dimension $15$, and the MARS model still reduced the computation in that
problem as well. The conclusion that MARS is a good surrogate model for any
forward model in any problem can not be made. However, as the model is non-
parametric, the approximation could work for similar problems with complicated
and computationally inefficient forward models in a Bayesian setting with
similar input spaces.

### Acknowledgments

We acknowledge support from the Centre for Geophysical Forecasting - CGF
(Norwegian Research Council grant no. 309960) and the GAMES consortium (grant
no. 294404) at the Norwegian University of Science and Technology (NTNU).

### Conflicts of interest/Competing interests

The authors declare no potential conflict of interests.

##  Appendix A Efficient sampling and evaluation of Gaussian distribution
using FFT

This section outlines the theory and application of circulant matrices and
fast Fourier transform (FFT) in generating Gaussian random field samples and
efficiently evaluating probability density functions.

#### Circulant matrices and FFT

The general idea is to use structure of circulant matrices and perform matrix-
vector multiplication in a lower-dimensional space in the Fourier domain. A
circulant matrix has the following structure

$\boldsymbol{C}=\begin{pmatrix}c_{0}&c_{n-1}&\dots&c_{1}\\\
c_{1}&c_{0}&\dots&c_{2}\\\ c_{2}&c_{1}&\dots&c_{3}\\\
\vdots&\vdots&\ddots&\vdots\\\ c_{n-1}&c_{n-2}&\dots&c_{0}\\\ \end{pmatrix}$

Eigenvalues of a circulant matrix can be obtained by the discrete Fourier
transform (DFT) of the first column of ${\boldsymbol{C}}$, and matrix
${\boldsymbol{C}}$ is fully specified by its one row or column, meaning all
the information that we need is stored in its basis. Matrix of eigenvalues of
a matrix such as ${\boldsymbol{C}}$, is denoted as ${\boldsymbol{\Lambda}}$.

We now assume our covariance matrix ${\boldsymbol{\Sigma}}$ is a circulant
matrix. [nolist, inline]Check the Block circulant assumption

We know there exists a decomposition
${\boldsymbol{\Sigma}}={\boldsymbol{C}}^{1/2}{\boldsymbol{C}}^{1/2},$ since
${\boldsymbol{\Sigma}}$ is assumed to be symmetric, positive definite matrix.
Same is true for the inverse,
${\boldsymbol{\Sigma}}^{-1}={\boldsymbol{C}}^{-1/2}{\boldsymbol{C}}^{-1/2}.$
Let $z\sim N(0,1)$. Using the unitary or Fourier matrix, $F$, one can then
write,

${\boldsymbol{C}}^{1/2}={\boldsymbol{F}}{\boldsymbol{\Lambda}}^{1/2}{\boldsymbol{F}}^{H}.$

To compute the product ${\boldsymbol{C}}^{-1/2}{\boldsymbol{z}}$,

${\boldsymbol{C}}^{1/2}{\boldsymbol{z}}={\boldsymbol{F}}{\boldsymbol{\Lambda}}^{-1/2}{\boldsymbol{F}}^{H}=\frac{1}{\sqrt{n}}(\text{idft2}({\boldsymbol{z}})\odot\sqrt{\text{dft2}({\boldsymbol{c}})})lin$

where ${\boldsymbol{c}}$ is the base of the circulant matrix. Here
dft2$(\cdot)$ and idft2$(\cdot)$ are respectively two-dimensional discrete
Fourier transform and inverse discrete Fourier transform. For more details, we
encourage the readers to have a look at [45].

####  A.0.1 Constructing the distance matrix

Here we explain how the distance matrix can be constructed such that it
represents a base of a circulant matrix. We make use of topological knowledge,
namely a torus. For a more detailed explanation we encourage a curious reader
to go deeper into reading the references, for example [14] provides some
illustrations and detailed explanations.

Assuming our covariance matrix has a circulant structure it is fully defined
by its first row or column. This is called the basis of the matrix and all of
the information that we need is stored there. Covariance matrix is made up of
variances and correlation structure. For the correlation structure, it is
necessary to define how we will measure distances between the points. And this
is a crucial step. We imagine we are on a torus and construct the distances
accordingly.

We measure distance from our corner point, corresponding to the first element
of the matrix, $C_{1,1}$ to all the other points (this will correspond to a
row of a distance matrix, later used in our covariance matrix). And since we
are on a torus, we are looking to measure the minimum distances. We are on a
circulant shaped geometric body, hence, once we are halfway over across the
ellipse/circle (cross-section of an empty torus), we are closer to the points
of choice from the “other side” of the circular body.

One can either choose to implement this as a double for-loop, where one goes
over the grid points on a torus and calculates the distance (and also
correlation function values), or one can represent this in a vector form and
perform calculations through vector and matrix calculations. Details of the
latter can be found in [14]. Additionally, depending on the programming
language one chooses to use, documentation should be thoroughly read. This is
because implementation of the discrete FFT, both one and multi-dimensional,
differs between programming languages and the way normalization is performed
or included can be different. Packages and libraries for performing two-
dimensional FFT are available in both R [41] and Python.

#### Sampling

Depending on the programming language, we can make use of the different
packages available in for example R or python, which enable us to calculate
the one or two-dimensional dft, necessary for computing the vector-matrix
product and obtaining our samples.

Assuming the size of the grid, $\mathcal{D}$, to be $n_{x}\times n_{y}$,
covariance
${\boldsymbol{\Sigma}}={\boldsymbol{\sigma}}^{2}\rho({\boldsymbol{D}})$, where
${\boldsymbol{D}}$ is the distance matrix computed as described in A.0.1. The
algorithm for simulating a sample of a Gaussian random field is shown in
Algorithm 2.

Algorithm 2 Sampling from $N({\boldsymbol{\mu}},{\boldsymbol{\Sigma}})$ using
circulant matrices and FFT.

1:Construct distance matrix ${\boldsymbol{D}}$ as in A.0.1

2:$\tilde{c}=\sigma^{2}\rho({\boldsymbol{D}})$ $\triangleright$ Constructing
the base of a circulant matrix

3:${\boldsymbol{z}}\sim N(0,1)$, $\triangleright$ Generating an independent
sample

4:${\boldsymbol{\Lambda}}^{sr}=\sqrt{\text{Re}(\text{dft2}(\tilde{c}))}$
$\triangleright$ Finding the eigenvalues through DFT

5:${\boldsymbol{x}}^{*}={\boldsymbol{\mu}}+\text{Re}(\text{dft2}({\boldsymbol{\Lambda}}^{sr}\odot\text{idft2}({\boldsymbol{z}})))$
$\triangleright$ Constructing a sample ${\boldsymbol{x}}^{*}\sim
N({\boldsymbol{\mu}},{\boldsymbol{\Sigma}})$

#### Evaluating

We wish to evaluate Gaussian pdf, which is given as,

$p(\mathbf{x};\mathbf{\mu},\mathbf{\Sigma})=\frac{1}{\sqrt{(2\pi)^{k}|\mathbf{\Sigma}|}}\exp\left(-\frac{1}{2}(\mathbf{x}-\mathbf{\mu})^{T}\mathbf{\Sigma}^{-1}(\mathbf{x}-\mathbf{\mu})\right).$

This gives the log-pdf being:

$\log(p(\mathbf{x}))\propto-\frac{1}{2}(\mathbf{x}-\mathbf{\mu})^{T}\mathbf{\Sigma}^{-1}(\mathbf{x}-\mathbf{\mu})$

If we for notational simplicity assume $\mu=0$, we have that
$\log(p(\mathbf{x}))\propto-\frac{1}{2}\mathbf{x}^{T}\mathbf{\Sigma}^{-1}\mathbf{x}$
which can then be written as

$\log(p(\mathbf{x}))\propto-\frac{1}{2}\mathbf{x}^{T}\mathbf{{\boldsymbol{C}}}^{-1/2}\mathbf{{\boldsymbol{C}}}^{-1/2}\mathbf{x}.$
(8)

Thus, We can compute equation (8) in an efficient way, making use of the
theory presented above for sampling. Details can be found in Algorithm 3.
Here, we are evaluating a $N(\mbox{\boldmath{$0$}},{\boldsymbol{\Sigma}})$,
${\boldsymbol{\mu}}=\mbox{\boldmath{$0$}}$, otherwise when we have
$N({\boldsymbol{\mu}},{\boldsymbol{\Sigma}})$, we let
${\boldsymbol{v}}=({\boldsymbol{x}}^{c}-{\boldsymbol{\mu}})$.

Algorithm 3 Evaluating a Gaussian pdf using circulant matrices and FFT.

1:Construct distance matrix ${\boldsymbol{D}}$ as in A.0.1

2:$\tilde{c}=\sigma^{2}\rho({\boldsymbol{D}})$ $\triangleright$ Previously
constructed base of a circulant matrix

3:${\boldsymbol{v}}=({\boldsymbol{x}}^{c}-{\boldsymbol{\mu}})$

4:${\boldsymbol{\Lambda}}^{sr}=\sqrt{\text{Re}(\text{dft2}(\tilde{c}))}$
$\triangleright$ Finding the eigenvalues through DFT

5:${\boldsymbol{u}}=\text{Re}(\text{dft2}({\boldsymbol{\Lambda}}^{-sr}\odot\text{idft2}({\boldsymbol{v}})))$
$\triangleright$ Computing ${\boldsymbol{C}}^{-1/2}{\boldsymbol{v}}$

6:$\log(p(\mathbf{x}))_{FFT}=-\frac{1}{2}{\boldsymbol{u}}^{T}{\boldsymbol{u}}$
$\triangleright$ Equivalent to equation (8)

##  Appendix B Non-parametric kernel regression

NPKR (Non-parametric kernel regression) approximates the response $f$ at a new
data point by a weighted average of the data points in the training set. There
are numerous ways to determine these weights. In this work, the Nadaraya-
Watson estimator

$\widehat{f}({\boldsymbol{x}})=\sum_{i=1}^{n}w_{i}({\boldsymbol{x}})f_{i}$ (9)

will be used with weights

$w_{i}({\boldsymbol{x}})=\frac{K_{{\boldsymbol{b}}}({\boldsymbol{x}}-{\boldsymbol{x}}_{i})}{\sum_{l=1}^{n}K_{{\boldsymbol{b}}}({\boldsymbol{x}}-{\boldsymbol{x}}_{l})}\quad
i=1,2,...,n.$ (10)

The function $K_{{\boldsymbol{b}}}({\boldsymbol{x}}-{\boldsymbol{x}}_{i})$ is
a Gaussian kernel with fixed bandwidths for each predictor
${\boldsymbol{b}}=[b_{1},b_{2},...,b_{p}]$. When using the Gaussian kernel,
the bandwidths are the standard deviation of the density. The weight in
equation (10) decrease with the distance between the new data point and the
data points in the training set. For $p$ covariates, $K_{{\boldsymbol{b}}}$ is
chosen to be the product of $p$ one dimensional Gaussian kernels:

$\displaystyle
K_{{\boldsymbol{b}}}({\boldsymbol{x}}-{\boldsymbol{x}}_{i})=\prod_{j=1}^{p}k_{b_{j}}(x_{j}-x_{ji})=\prod_{j=1}^{p}\frac{1}{b_{j}}k\left(\frac{x_{j}-x_{ji}}{b_{j}}\right)=\prod_{j=1}^{p}\frac{1}{\sqrt{2\pi}b_{j}}e^{-\frac{(x_{j}-x_{ji})^{2}}{2b_{j}^{2}}}.$

The Gaussian kernel is a popular non-compact kernel [28]. The function npreg
from the R package “np” [29] is used to fit the NPKR model. The reader can
refer to [42] for an in-depth exploration and description of the model

The bias-variance trade-off asserts itself when selecting the bandwidths.
Selecting small bandwidths tends to lower the bias and increase the variance
as the nearest points are valued a lot. On the other hand, selecting too large
bandwidths leads to a high bias and low variance, since this tends to include
values far away from the new data points in the estimate. For the NPKR model
used in this work, the bandwidths are found using least squares cross-
validation. This is the only predetermined property of an NPKR model because
the weights need to be determined every time a prediction of a new data point
is made. Much of the work is therefore done at evaluation time [28].

Using an NPKR model has both advantages and disadvantages. An advantage lies
in its flexibility as it does not assume a specific functional form for the
response, which is particularly beneficial in high-dimensional scenarios.
Additionally, the intuitive nature of averaging nearby data points contributes
to its appeal. The model, represented by equations (9) and (10), constitutes a
linear combination of smooth functions, resulting in an overall smooth model.

However, a notable disadvantage of NPKR is its susceptibility to bias at
boundaries, as noted in prior literature [15, 28]. When making predictions
near boundaries, the weighted average tends to favor data points away from the
boundary, potentially leading to inaccurate predictions, especially if the
actual response is decreasing. This bias becomes more pronounced in higher
dimensions due to the increased proportion of data points near the boundary.
Although one approach to mitigate this bias involves using locally weighted
linear regression instead of a weighted average, this method is not explored
in this work. The boundary bias issue also extends to unevenly distributed
training data, where predictions may be skewed towards the direction with more
data points. The reliance of NPKR on the specific data points in the training
set further underscores the importance of evenly spaced training data to
potentially enhance the model’s performance.

###  B.1 Approximating the gradient of the forward model

To use the MALA introduced in Section 3.2.2, the gradient
$\nabla{\boldsymbol{h}}$ is needed. However, the forward model is treated as a
black box. Consequently, the partial derivatives are unavailable. A numerical
approximation to the gradient can be found by adding small perturbations,
$\varepsilon$, to only one of the covariates while the rest are fixed. The
difference in the forward model with and without the perturbation is divided
by $\varepsilon$. For the numerical partial derivative of $h_{R_{0}}$ with
respect to the covariate $x_{g}$ at data point
${\boldsymbol{x}}_{i}=[x_{ig},x_{io},x_{ic},x_{id}]$ this is

$\frac{\partial h_{R_{0}}}{\partial
x_{g}}({\boldsymbol{x}}_{i})=\frac{h_{R_{0}}({\boldsymbol{x}}_{i}+[\varepsilon,0,0,0])-h_{R_{0}}({\boldsymbol{x}}_{i})}{\varepsilon},$

with $\varepsilon=0.0001$. However, the forward function needs to be evaluated
to calculate the numerical partial derivatives, which require a lot of
computation time. In this thesis, the gradients of the four models discussed
so far are tested for approximating the gradient of the forward model,
$\nabla{\boldsymbol{h}}$, in order to reduce the computation time. The
gradient of the MARS models and NPKR models are described in B. Because the
analytical gradient is unavailable, the approximations are compared to the
numerical approximations of the gradient instead. For simplicity, the
numerical approximation to the gradient is denoted by $\nabla{\boldsymbol{h}}$
and referred to as the gradient of the forward model throughout this section.

In addition, MARS models are created with partial derivatives as responses.
That is, six MARS models are fitted with each of the partial derivatives
$\frac{\partial h_{R_{0}}}{\partial x_{g}}$, $\frac{\partial
h_{R_{0}}}{\partial x_{o}}$, $\frac{\partial h_{R_{0}}}{\partial x_{c}}$,
$\frac{\partial h_{G}}{\partial x_{g}}$, $\frac{\partial h_{G}}{\partial
x_{o}}$, and $\frac{\partial h_{G}}{\partial x_{c}}$ as responses. As MCMC is
only performed with
${\boldsymbol{x}}=[{\boldsymbol{x}}_{g},{\boldsymbol{x}}_{o},{\boldsymbol{x}}_{c}]$
the partial derivative of ${\boldsymbol{h}}$ with respect to the depth is not
needed. The six MARS models fitted to the partial derivatives are collectively
denoted as $\text{MARS}_{\nabla}$.

The K-means approximation to the gradient of an NPKR model is again tested.
The procedure is described in B. A simple grid search of the number of
clusters $\in[3,5,10]$ and the number of closest neighbours to average
$\in[3,5,10]$ was performed. The computation time decreased as the number of
clusters increased. Therefore, the number of clusters was set to $10$. When
the number of clusters was $10$, the lowest MSE and the highest correlation
was achieved when using $5$ neighbours for both
$\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{1000}}$ and
$\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{4000}}$.

The accuracy of the approximations to the gradient is the mean correlation

$\text{Corr}\left(\nabla{\boldsymbol{h}},\widehat{\nabla{\boldsymbol{h}}}\right):=\frac{1}{6}\sum_{j\in\\{g,o,c\\}}\left(\text{Corr}\left(\frac{\partial
h_{R_{0}}}{\partial x_{j}},\widehat{\frac{\partial h_{R_{0}}}{\partial
x_{j}}}\right)+\text{Corr}\left(\frac{\partial h_{G}}{\partial
x_{j}},\widehat{\frac{\partial h_{G}}{\partial x_{j}}}\right)\right)$

and mean MSE

$\text{MSE}\left(\nabla{\boldsymbol{h}},\widehat{\nabla{\boldsymbol{h}}}\right):=\frac{1}{6}\sum_{j\in\\{g,o,c\\}}\left(\text{MSE}\left(\frac{\partial
h_{R_{0}}}{\partial x_{j}},\widehat{\frac{\partial h_{R_{0}}}{\partial
x_{j}}}\right)+\text{MSE}\left(\frac{\partial h_{G}}{\partial
x_{j}},\widehat{\frac{\partial h_{G}}{\partial x_{j}}}\right)\right)$

for all the covariates $j$. The correlation, MSE and computation time are
reported in Table 3. Again, the computation time is the average of $50$
computation times. One should keep in mind when reading Table 3, that the
correlation and MSE are calculated using the numerical approximations as
$\nabla{\boldsymbol{h}}$.

As per Table 3, $\text{MARS}_{\nabla}$ emerges as the most effective in
approximating the forward model gradient, displaying superior correlation,
MSE, and computation time. In Chapter 5, MALA employs this approximation. The
gradients of the MARS models, $\nabla\widehat{{\boldsymbol{h}}}_{\text{MARS}}$
and $\nabla\widehat{{\boldsymbol{h}}}_{\text{MARS}_{\text{OF}}}$, rank second
and third highest in correlation, and second and third lowest in MSE,
following $\text{MARS}_{\nabla}$. Notably, in contrast to the 2D case,
$\nabla\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ outperforms
$\nabla\widehat{{\boldsymbol{h}}}_{\text{MARS}_{\text{OF}}}$ marginally in
correlation and MSE. However, computation time for these gradients is
considerably high. Integration of binary search in the partial derivative
extraction algorithm reduces computation time, yet its efficiency varies with
interaction terms in the MARS model. For a MARS model with approximately 40
non-interaction terms, computation time is around 10% that of
$\nabla\widehat{{\boldsymbol{h}}}_{\text{MARS}}$. However, reducing the MARS
model’s degree exacerbates MSE and correlation. Increasing data points in the
NPKR model enhances correlation and reduces MSE for its gradient but extends
computation time. While using the K-means approach expedites approximating
$\nabla\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{4000}}$, it slows down
approximating $\nabla\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{1000}}$.

The computation time of
$\nabla\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{1000}}$ was close to the
computation time of $\nabla\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ and lower
than the computation time for the gradient of the larger MARS model. The
K-means approach with $10$ clusters and $5$ neighbours is not particularly
fast. It was however, a bit better at predicting the gradient of the forward
model compared to $\nabla\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{1000}}$ and
$\nabla\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{4000}}$.

$\widehat{\nabla{\boldsymbol{h}}}$ | correlation | MSE $[10^{-5}]$ | computation time [sec]
---|---|---|---
$\text{MARS}_{\nabla}$ | 0.984 | 0.2 | 0.32
$\nabla\widehat{{\boldsymbol{h}}}_{\text{MARS}}$ | 0.879 | 1.7 | 20.59
$\nabla\widehat{{\boldsymbol{h}}}_{\text{MARS}_{\text{OF}}}$ | 0.872 | 1.9 | 33.73
$\nabla\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{1000}}$ | 0.621 | 6.2 | 22.41
$\widehat{\nabla{\boldsymbol{h}}_{\text{NPKR}_{1000}}}$ (K-means) | 0.679 | 6.1 | 43.58
$\nabla\widehat{{\boldsymbol{h}}}_{\text{NPKR}_{4000}}$ | 0.702 | 4.7 | 62.19
$\widehat{\nabla{\boldsymbol{h}}_{\text{NPKR}_{4000}}}$ (K-means) | 0.844 | 3.7 | 54.28
Table 3: Correlation, MSE and computation time for approximations
$\nabla\widehat{{\boldsymbol{h}}}$.

(a) ACF plot from MCMC computations with $q_{1}$ as proposal distribution.

(b) ACF plot from MCMC computations with $q_{2}$ as proposal distribution.

(c) ACF plot from MCMC computations with $q_{3}$ as proposal distribution.

(d) ACF plot from MCMC computations with $q_{4}$ as proposal distribution.

(e) ACF plot from MCMC computations with $q_{1}$ as proposal distribution.

(f) ACF plot from MCMC computations with $q_{2}$ as proposal distribution.

(g) ACF plot from MCMC computations with $q_{3}$ as proposal distribution.

(h) ACF plot from MCMC computations with $q_{4}$ as proposal distribution.

(i) ACF plot from MCMC computations with $q_{1}$ as proposal distribution.

(j) ACF plot from MCMC computations with $q_{2}$ as proposal distribution.

(k) ACF plot from MCMC computations with $q_{3}$ as proposal distribution.

(l) ACF plot from MCMC computations with $q_{4}$ as proposal distribution.

Figure 13: ACF plot up to lag $k=1000$ of every $100$th value of
$\\{x_{g,t}^{l}\\}_{t=0}^{1000000}$ (first row),
$\\{x_{o,t}^{l}\\}_{t=0}^{1000000}$ (second row) and
$\\{x_{c,t}^{l}\\}_{t=0}^{1000000}$ (third row) at location $l=435$ using four
MH algorithms. The MH algorithm with proposal distribution $q_{1}$, $q_{2}$,
$q_{3}$ and $q_{4}$ was used to get the results in the first, second, third
and fourth columns respectively.

###  B.2 Adjusting for the uncertainty added to the forward model

To adjust for the additional error caused by replacing
${\boldsymbol{h}}({\boldsymbol{x}})$ by
$\widehat{{\boldsymbol{h}}}_{\text{MARS}}({\boldsymbol{x}})$, the relationship
between the variable of interest and the data can be modified to

${\boldsymbol{y}}_{s}=\widehat{{\boldsymbol{h}}}_{\text{MARS}}({\boldsymbol{x}})+\widetilde{{\boldsymbol{\omega}}},\quad\widetilde{{\boldsymbol{\omega}}}\sim\mathcal{N}(0,\widetilde{{\boldsymbol{\Omega}}}_{0}),$

where $\widetilde{{\boldsymbol{\Omega}}}_{0}$ is the covariance matrix
${\boldsymbol{\Omega}}_{0}$ with empirical variance and covariance are added
to it. A homoscedastic independence of the covariates is assumed. Plots of
$|{\boldsymbol{h}}({\boldsymbol{x}})-\widehat{{\boldsymbol{h}}}_{\text{MARS}}({\boldsymbol{x}})|$
against the covariates used to investigate dependency between the absolute
error and the covariates, indicated that this is not an unreasonable
assumption. The idea is that an average sample variance in the approximation
partially compensates for using $\widehat{{\boldsymbol{h}}}_{\text{MARS}}$.
The adjusted covariance matrix is, where the subscript MARS is omitted for
simplicity,

$\widetilde{{\boldsymbol{\Omega}}}_{0}={\boldsymbol{\Omega}}_{0}+\begin{bmatrix}\widehat{\text{Var}}\left[\widehat{h_{R_{0}}}\right]&\widehat{\text{Cov}}\left[\widehat{h_{R_{0}}},\widehat{h_{G}}\right]\\\
\widehat{\text{Cov}}\left[\widehat{h_{R_{0}}},\widehat{h_{G}}\right]&\widehat{\text{Var}}\left[\widehat{h_{G}}\right],\end{bmatrix}$

with

$\widehat{\text{Var}}\left[\widehat{h_{R_{0}}}\right]=\frac{1}{n}\sum_{i=1}^{n}\left(\widehat{h_{{R_{0}}_{i}}}-h_{{R_{0}}_{i}}\right)^{2},\quad\widehat{\text{Var}}\left[\widehat{h_{G}}\right]=\frac{1}{n}\sum_{i=1}^{n}\left(\widehat{h_{G_{i}}}-h_{G_{i}}\right)^{2},$

and

$\widehat{\text{Cov}}\left[\widehat{h_{R_{0}}},\widehat{h_{G}}\right]=\frac{1}{n}\sum_{i=1}^{n}\left(\widehat{h_{{R_{0}}_{i}}}-h_{{R_{0}}_{i}}\right)\left(\widehat{h_{G_{i}}}-h_{G_{i}}\right).$

Here $\widehat{h_{{R_{0}}_{i}}}$ is the prediction for test data point $i$ and
$h_{h_{{R_{0}}_{i}}}$ denotes the response at test data point $i$. The error-
adjusted covariance matrix for the seismic AVO log-likelihood is thus

$\widetilde{{\boldsymbol{\Omega}}}=\begin{bmatrix}\left[\widetilde{{\boldsymbol{\Omega}}}_{0}\right]&0&\ldots&0&0\\\
0&\left[\widetilde{{\boldsymbol{\Omega}}}_{0}\right]&\ldots&0&0\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&\ldots&\left[\widetilde{{\boldsymbol{\Omega}}}_{0}\right]&0\\\
0&0&\ldots&0&\left[\widetilde{{\boldsymbol{\Omega}}}_{0}\right]\\\
\end{bmatrix}.$ (11)

(a) Trace plot from MCMC computations with $q_{1}$ as proposal distribution.

(b) Trace plot from MCMC computations with $q_{2}$ as proposal distribution.

(c) Trace plot from MCMC computations with $q_{3}$ as proposal distribution.

(d) Trace plot from MCMC computations with $q_{4}$ as proposal distribution.

(e) Trace plot from MCMC computations with $q_{1}$ as proposal distribution.

(f) Trace plot from MCMC computations with $q_{2}$ as proposal distribution.

(g) Trace plot from MCMC computations with $q_{3}$ as proposal distribution.

(h) Trace plot from MCMC computations with $q_{4}$ as proposal distribution.

(i) Trace plot from MCMC computations with $q_{1}$ as proposal distribution.

(j) Trace plot from MCMC computations with $q_{2}$ as proposal distribution.

(k) Trace plot from MCMC computations with $q_{3}$ as proposal distribution.

(l) Trace plot from MCMC computations with $q_{4}$ as proposal distribution.

Figure 14: Trace plots of every $100$th value of
$\\{x_{g,t}^{l}\\}_{t=0}^{1000000}$ (first row),
$\\{x_{o,t}^{l}\\}_{t=0}^{1000000}$ (second row) and
$\\{x_{c,t}^{l}\\}_{t=0}^{1000000}$ (third row) at location $l=435$ using four
MH algorithms. The MH algorithm with proposal distribution $q_{1}$, $q_{2}$,
$q_{3}$ and $q_{4}$ was used to get the results in the first, second, third
and fourth columns respectively.

##  Appendix C Additional results in simulation studies

Figure 13 displays ACF plots for every $100$th value of sequences
$\\{x_{g,t}^{l}\\}_{t=0}^{1000000}$, $\\{x_{o,t}^{l}\\}_{t=0}^{1000000}$, and
$\\{x_{c,t}^{l}\\}_{t=0}^{1000000}$ at location $l=435$ using four MH
algorithms. These plots align with Table 2. The ACF plots in Figure 13
indicate that the MH algorithm with $q_{3}$ as the proposal distribution
produces the least correlated samples. It shows samples approximately $15,000$
time iterations apart, indicating insignificant correlation. Compared to the
MH algorithm with $q_{1}$ as the proposal distribution, the ACF plots for
$q_{3}$ decrease much faster. Additionally, ACF plots for the MH algorithm
with proposal distribution $q_{2}$ also decrease rapidly compared to $q_{1}$
and $q_{4}$. While autocorrelation decreases similarly for $q_{1}$ and $q_{4}$
in gas plots, MALA decreases noticeably faster in oil and clay plots compared
to $q_{1}$.

The trace plots depicted in Figure 14 are consistent with the findings
presented in Table 2. They illustrate that the MH algorithm employing $q_{3}$
as the proposal distribution exhibits optimal mixing, as evidenced by the
trajectory frequently traversing around a central value. Similarly, the trace
plots suggest that the MH algorithm utilizing proposal distribution $q_{2}$
achieves good mixing, albeit with slightly less coverage compared to the most
efficient algorithm. Moreover, the trace plots highlight that the MALA mixes
notably better than the standard random walk but not as effectively as the MH
algorithms employing $q_{2}$ and $q_{3}$ as proposal distributions.

However, the trace plots indicate that the MH algorithm with $q_{1}$ as the
proposal distribution has not yet converged, as evidenced by a trend
suggesting ongoing exploration of the posterior. Comparison of the trace plots
for the standard random walk with those of other algorithms reveals that the
former has not reached convergence; for instance, the trace plot in Figure
14(a) has not approached values close to $-10$, an area visited by all three
other algorithms.

At this particular location, the ACF and trace plots corroborate the results
in Table 2. It’s noteworthy that mixing behaviors, as reflected in trace and
ACF plots, can vary across different locations. Nonetheless, the table’s
Effective Sample Size (ESS) represents the mean ESS across all locations,
indicating that the MH algorithm employing $q_{3}$ as the proposal
distribution exhibited the most favorable overall mixing.

## References

  * Abrahamsen et al., [2018] Abrahamsen, P., Kvernelv, V., and Barker, D. (2018).  Simulation of gaussian random fields using the fast fourier transform (fft).  In ECMOR 16th European Conference on the Mathematics of Oil Recovery, volume 2018 (1), pages 1–14. European Association of Geoscientists & Engineers.
  * Alquier, [2024] Alquier, P. (2024).  User-friendly introduction to pac-bayes bounds.  Foundations and Trends® in Machine Learning, 17(2):174–303.
  * Alquier and Ridgway, [2020] Alquier, P. and Ridgway, J. (2020).  Concentration of tempered posteriors and of their variational approximations.  The Annals of Statistics, 48(3):1475–1497.
  * Avseth et al., [2008] Avseth, P., Dræge, A., van Wijngaarden, A.-J., Johansen, T. A., and Jørstad, A. (2008).  Shale rock physics and implications for avo analysis: A north sea demonstration.  The Leading Edge, 27(6):788–797.
  * Avseth and Lehocki, [2016] Avseth, P. and Lehocki, I. (2016).  Combining burial history and rock-physics modeling to constrain avo analysis during exploration.  The Leading Edge, 35(6):528–534.
  * Avseth et al., [2010] Avseth, P., Mukerji, T., Mavko, G., and Dvorkin, J. (2010).  Rock-physics diagnostics of depositional texture, diagenetic alterations, and. reservoir heterogeneity in high-porosity siliciclastic sediments and rocks - a review of selected models and suggested work flows.  Geophysics, 75:A31–A47.
  * Bissiri et al., [2016] Bissiri, P. G., Holmes, C. C., and Walker, S. G. (2016).  A general framework for updating belief distributions.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 78(5):1103–1130.
  * Bosch et al., [2010] Bosch, M., Mukerji, T., and Gonzalez, E. F. (2010).  Seismic inversion for reservoir properties combining statistical rock physics and geostatistics: A review.  Geophysics, 75:A165–A176.
  * Buland et al., [2003] Buland, A., Kolbjørnsen, O., and Omre, H. (2003).  Rapid spatially coupled avo inversion in the fourier domain.  Geophysics, 68(3):824–836.
  * Chan and Wood, [1997] Chan, G. and Wood, A. T. (1997).  Algorithm as 312: An algorithm for simulating stationary gaussian random fields.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 46(1):171–181.
  * Chen et al., [2013] Chen, M., Sun, Y., Fu, P., Carrigan, C. R., Lu, Z., Tong, C. H., and Buscheck, T. A. (2013).  Surrogate-based optimization of hydraulic fracturing in pre-existing fracture networks.  Computers & Geosciences, 58:69–79.
  * Chen et al., [2014] Chen, M., Tompson, A. F., Mellors, R. J., Ramirez, A. L., Dyer, K. M., Yang, X., and Wagoner, J. L. (2014).  An efficient bayesian inversion of a geothermal prospect using a multivariate adaptive regression spline method.  Applied Energy, 136:619–627.
  * Cotter et al., [2013] Cotter, S., Roberts, G., Stuart, A., and White, D. (2013).  Mcmc methods for functions: Modifying old algorithms to make them faster.  Statistical Science, 28(3):424.
  * Davies and Bryant, [2013] Davies, T. M. and Bryant, D. (2013).  On circulant embedding for gaussian random fields in r.  Journal of Statistical Software, 55:1–21.
  * Efron and Hastie, [2016] Efron, B. and Hastie, T. (2016).  Computer Age Statistical Inference Algorithms, Evidence and Data Science.  Cambridge University Press.
  * Eidsvik et al., [2004] Eidsvik, J., Avseth, P., Omre, H., Mukerji, T., and Mavko, G. (2004).  Stochastic reservoir characterization using prestack seismic data.  Geophysics, 69(4):978–993.
  * Forberg et al., [2021] Forberg, O. B., Kjøsnes, Ø., and Omre, H. (2021).  Bayesian seismic avo inversion using a laterally coupled multimodal prior model.  IEEE Transactions on Geoscience and Remote Sensing, 60:1–16.
  * Friedman, [1991] Friedman, J. H. (1991).  Multivariate adaptive regression splines.  The annals of statistics, 19(1):1–67.
  * Gelman et al., [1997] Gelman, A., Gilks, W. R., and Roberts, G. O. (1997).  Weak convergence and optimal scaling of random walk metropolis algorithms.  The annals of applied probability, 7(1):110–120.
  * Givens and Hoeting, [2012] Givens, G. H. and Hoeting, J. A. (2012).  Computational Statistics, volume 703.  John Wiley & Sons.
  * Grana, [2022] Grana, D. (2022).  Bayesian rock-physics inversion with kumaraswamy prior models.  Geophysics, 87(3):M87–M97.
  * [22] Grana, D., Azevedo, L., De Figueiredo, L., Connolly, P., and Mukerji, T. (2022a).  Probabilistic inversion of seismic data for reservoir petrophysical characterization: Review and examples.  Geophysics, 87:M199–M216.
  * [23] Grana, D., de Figueiredo, L., and Mosegaard, K. (2022b).  Markov chain monte carlo for petrophysical inversion.  Geophysics, 87(1):M13–M24.
  * Grana et al., [2021] Grana, D., Mukerji, T., and Doyen, P. (2021).  Seismic reservoir modeling: Theory, examples, and algorithms.  John Wiley & Sons.
  * Grünwald and Van Ommen, [2017] Grünwald, P. and Van Ommen, T. (2017).  Inconsistency of bayesian inference for misspecified linear models, and a proposal for repairing it.  Bayesian Analysis, 12(4):1069–1103.
  * Hairer et al., [2014] Hairer, M., Stuart, A., and Vollmer, S. (2014).  Spectral gaps for a metropolis–hastings algorithm in infinite dimensions.  The Annals of Applied Probability, 24(6):2455–2490.
  * Hammer et al., [2023] Hammer, H. L., Riegler, M. A., and Tjelmeland, H. (2023).  Approximate bayesian inference based on expected evaluation.  Bayesian Analysis, pages 1–22.
  * Hastie et al., [2009] Hastie, T., Tibshirani, R., and Friedman, J. (2009).  The Elements of Statistical Learning. Data Mining, Inference, and Prediction.  Springer, 2 edition.
  * Hayfield and Racine, [2008] Hayfield, T. and Racine, J. S. (2008).  Nonparametric econometrics: The np package.  Journal of Statistical Software, 27(5):1–32.
  * Hong and Martin, [2020] Hong, L. and Martin, R. (2020).  Model misspecification, bayesian versus credibility estimation, and gibbs posteriors.  Scandinavian Actuarial Journal, 2020(7):634–649.
  * Khoshkholgh et al., [2021] Khoshkholgh, S., Zunino, A., and Mosegaard, K. (2021).  Informed proposal monte carlo.  Geophysical Journal International, 226(2):1239–1248.
  * Khoshkholgh et al., [2022] Khoshkholgh, S., Zunino, A., and Mosegaard, K. (2022).  Full-waveform inversion by informed-proposal monte carlo.  Geophysical Journal International, 230(3):1824–1833.
  * Knoblauch et al., [2022] Knoblauch, J., Jewson, J., and Damoulas, T. (2022).  An optimization-centric view on bayes’ rule: Reviewing and generalizing variational inference.  Journal of Machine Learning Research, 23(132):1–109.
  * Mai, [2023] Mai, T. T. (2023).  A reduced-rank approach to predicting multiple binary responses through machine learning.  Statistics and Computing, 33:136.
  * Mai and Alquier, [2017] Mai, T. T. and Alquier, P. (2017).  Pseudo-bayesian quantum tomography with rank-adaptation.  Journal of Statistical Planning and Inference, 184:62–76.
  * Malinverno and Briggs, [2004] Malinverno, A. and Briggs, V. A. (2004).  Expanded uncertainty quantification in inverse problems: Hierarchical bayes and empirical bayes.  Geophysics, 69(4):1005–1016.
  * Matsubara et al., [2022] Matsubara, T., Knoblauch, J., Briol, F.-X., and Oates, C. J. (2022).  Robust generalised bayesian inference for intractable likelihoods.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 84(3):997–1022.
  * Milborrow, [2023] Milborrow, S. (2023).  earth: Multivariate Adaptive Regression Splines.  R package version 5.3.2 Milborrow. Derived from mda:mars by Trevor Hastie and Rob Tibshirani. Uses Alan Miller’s Fortran utilities with Thomas Lumley’s leaps wrapper.
  * Mosegaard and Tarantola, [1995] Mosegaard, K. and Tarantola, A. (1995).  Monte carlo sampling of solutions to inverse problems.  Journal of Geophysical Research: Solid Earth, 100(B7):12431–12447.
  * Pinski et al., [2015] Pinski, F. J., Simpson, G., Stuart, A. M., and Weber, H. (2015).  Algorithms for kullback–leibler approximation of probability measures in infinite dimensions.  SIAM Journal on Scientific Computing, 37(6):A2733–A2757.
  * R Core Team, [2022] R Core Team (2022).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Racine and Li, [2004] Racine, J. and Li, Q. (2004).  Nonparametric estimation of regression functions with both categorical and continuous data.  Journal of Econometrics, 119(1):99–130.
  * Rimstad et al., [2012] Rimstad, K., Avseth, P., and Omre, H. (2012).  Hierarchical bayesian lithology/fluid prediction: A north sea case study.  Geophysics, 77:B69–B85.
  * Rudolf and Sprungk, [2018] Rudolf, D. and Sprungk, B. (2018).  On a generalization of the preconditioned crank–nicolson metropolis algorithm.  Foundations of Computational Mathematics, 18:309–343.
  * Rue and Held, [2005] Rue, H. and Held, L. (2005).  Gaussian Markov random fields: theory and applications.  CRC press.
  * Spremić et al., [2024] Spremić, M., Eidsvik, J., and Avseth, P. (2024).  Bayesian rock-physics inversion using a localized ensemble-based approach—with an application to the alvheim field.  Geophysics, 89(2):R95–R108.
  * Ulrych et al., [2001] Ulrych, T. J., Sacchi, M. D., and Woodbury, A. (2001).  A bayes tour of inversion: A tutorial.  Geophysics, 66(1):55–69.