# Multiscale Sampling for the Inverse Modeling of Partial Differential
Equations

Alsadig Ali  Abdullah Al-Mamun  Felipe Pereira  Arunasalam Rahunanthan
Department of Mathematical Sciences, The University of Texas at Dallas,
Richardson, TX, 75080, USA  Institute of Natural Sciences, United
International University, Dhaka-1212, Bangladesh  Department of Mathematics
and Computer Science, Central State University, Wilberforce, OH, 45384, USA
Corresponding author: Felipe Pereira<EMAIL_ADDRESS>
###### Abstract

We are concerned with a novel Bayesian statistical framework for the
characterization of natural subsurface formations, a very challenging task.
Because of the large dimension of the stochastic space of the prior
distribution in the framework, typically a dimensional reduction method, such
as a Karhunen-Leove expansion (KLE), needs to be applied to the prior
distribution to make the characterization computationally tractable. Due to
the large variability of properties of subsurface formations (such as
permeability and porosity) it may be of value to localize the sampling
strategy so that it can better adapt to large local variability of rock
properties.

In this paper, we introduce the concept of multiscale sampling to localize the
search in the stochastic space. We combine the simplicity of a preconditioned
Markov Chain Monte Carlo method with a new algorithm to decompose the
stochastic space into orthogonal subspaces, through a one-to-one mapping of
the subspaces to subdomains of a non-overlapping domain decomposition of the
region of interest. The localization of the search is performed by a
multiscale blocking strategy within Gibbs sampling: we apply a KL expansion
locally, at the subdomain level. Within each subdomain, blocking is applied
again, for the sampling of the KLE random coefficients.

The effectiveness of the proposed framework is tested in the solution of
inverse problems related to elliptic partial differential equations arising in
porous media flows. We use multi-chain studies in a multi-GPU cluster to show
that the new algorithm clearly improves the convergence rate of the
preconditioned MCMC method. Moreover, we illustrate the importance of a few
conditioning points to further improve the convergence of the proposed method.

###### keywords:

Preconditioned MCMC , MCMC convergence , Inverse modeling , Multiscale
Sampling

##  1 Introduction

Markov chain Monte Carlo (MCMC) methods have important applications and have
experienced enormous developments since the original work of Metropolis-
Hastings [1, 2, 3]. These methods differ on the strategy used in their
sampling stage. They can be divided in two large classes, depending on whether
they use or not the gradient of the likelihood function. Methods that take
advantage of the gradient information typically show improvement in the
convergence rate to the equilibrium distribution. These methods include, among
many others, the MALA-Gibbs [4], Hamiltonian Monte Carlo [5], and the active
subspace method [6]. There are also Hessian-based procedures [7]. In the case
of porous media flows gradient calculations can be computationally very
expensive (see [8] and references therein).

Our motivation for this work is the gradient-free uncertainty quantification
in inverse problems associated with porous media flow problems in the field
scale. For this class of problems the likelihood calculation typically
requires the numerical solution of very large problems in fine computational
grids [9]. There are methods known as upscaling [10], that aim at solving an
approximate problem (with effective coefficients) on coarse grids, thus
reducing drastically the cost of the simulations. The idea of upscaling was
used in [11] to define a two-stage or preconditioned MCMC (see also [12] where
a procedure of this type was proposed). In the preconditioned MCMC a sample is
first tested through a coarse grid numerical simulation with upscaled
coefficients. If it passes this filter, then a full fine grid simulation has
to be performed to determine if the sample is accepted. Computationally this
is a quite effective procedure because they do not require gradient
calculations and samples can be discarded with a coarse grid simulation. The
methods of [11, 12] have been further investigated over the years. Some
developments of these methods include their application to flows in fractured
porous media [13], their multi-physics version [14], and their parallelization
in multi-core devices [15]. This procedure has also been successfully applied
in the solution of inverse problems in geophysics [16] and more recently a
multi-level version of [12] has been introduced [17]. We remark that although
the preconditioned MCMC is computationally more competitive than the
Metropolis-Hastings algorithm, it still shows slow convergence for large
dimensional problems [16]. Thus, the development of MCMC methods that show
good convergence properties and do not require gradient calculations remains
as an important area for research.

In this work, the new Multiscale Sampling Method (MSM) is proposed. It is
motivated by multiscale methods for the solution of second order elliptic
equations that are based on a domain decomposition (see [18, 19] and
references therein) that can produce solution for large problems taking
advantage of the solution of a family of smaller boundary value problems. In
these methods the domain of the equation is decomposed into non-overlapping
subdomains, local multiscale basis functions are computed for each subdomain,
and a global interface problem is constructed and solved to couple the local
solutions and produce the global solution. In the Multiscale Sampling Method
the domain of the partial differential equation is also decomposed into non-
overlapping subdomains and a local truncated Karhunen-Loève expansion (KLE)
[20] is used for each subdomain. The final stage in the construction of one
sample consists in applying a local averaging procedure to remove
discontinuities between adjacent subdomains. The localized sampling is
performed by Gibbs sampling [21].

We perform several multi-chain MCMC studies in a multi-GPU cluster to compare
the convergence of the preconditioned MCMC with and without multiscale
sampling for high-dimensional problems. We find that multiscale sampling has a
huge impact in improving convergence rates for all problems considered.
Moreover, acceptance rates also increase when the multiscale sampling is used.
One example is included to illustrate the improved convergence that results
from combining the proposed method with a few points for conditioning the
field of interest.

This work is organized as follows. We begin by describing the governing
equations for contaminant transport problems in Section 2. In Section 3 we
present a Bayesian framework for subsurface characterization and the KLE for
the dimensional reduction. We also recall methods for convergence assessment
of MCMC methods. In Section 4 we describe the proposed method. Kriging and
conditioning that will be used in our numerical studies are discussed in
Section 5. Numerical results from our experiments appear in Section 6. Our
conclusions appear in Section 7.

##  2 The Model Problem

###  2.1 Motivation

We consider a model for contaminant transport problems (or single-phase flow
problems) in a subsurface aquifer $\Omega$ with a heterogeneous permeability
field. In this model, the first equation of a system of governing equations is
an elliptic equation

$\displaystyle\begin{cases}\boldsymbol{u}~{}~{}~{}~{}~{}&=-k(\boldsymbol{x})\nabla
p~{}~{}~{}~{}~{}\text{ in }\Omega\\\
\nabla\cdot\boldsymbol{u}&=f\qquad\qquad~{}\text{ in }\Omega,\end{cases}$ (1)

where $\boldsymbol{u}$ and $p$ represent the Darcy velocity and the fluid
pressure, respectively, $k(\boldsymbol{x})$ is known as the absolute
permeability field of the rock (a positive definite tensor), and $f$
represents sources and sinks. The elliptic equation is coupled to a hyperbolic
equation

$\displaystyle\phi(\boldsymbol{x})\frac{\partial s(\boldsymbol{x},t)}{\partial
t}+\nabla\cdot[s(\boldsymbol{x},t)\boldsymbol{u}(\boldsymbol{x})]=0,$ (2)

where $s(\boldsymbol{x},t)$ is the contaminant concentration in the water and
$\phi(\boldsymbol{x})$ is the porosity of the rock.

The aquifer may contain many monitoring and injection wells. Our goal is to
characterize the permeability field of the domain of interest by using
available fractional flow data defined by

$\displaystyle
F(t)=1-\frac{\int_{{\partial\Omega}_{\text{out}}}u_{n}(\boldsymbol{x})s(\boldsymbol{x},t)~{}dy}{\int_{{\partial\Omega}_{\text{out}}}u_{n}(\boldsymbol{x})~{}dy},$

where ${\partial\Omega}_{\text{out}}$ and $u_{n}(\boldsymbol{x})$ are the well
outflow boundary and the normal components of the velocity field,
respectively. More details about single-phase flow problems can be found in
[22, 23]. In this paper we illustrate the proposed method in terms of the
elliptic equation (1).

###  2.2 Variational Formulation of the Pressure Equation

In this work we consider $\Omega\subset\mathbb{R}^{2}$, a bounded domain with
a Lipschitz boundary $\partial\Omega$. For problems in $\mathbb{R}^{3}$ a
formulation similar to the one we describe here is also applicable. The
velocity-pressure system is given by Eq. (1). In porous media flow
applications typical boundary conditions that occur are Dirichlet (the
pressure is given) and Neumann (the normal component of the velocity is
specified), which are expressed as

$p=g_{p}\in
H^{\frac{1}{2}}(\partial\Omega_{p}),\qquad\boldsymbol{u}\cdot\boldsymbol{\hat{n}}=g_{\boldsymbol{u}}\in
H^{-\frac{1}{2}}(\partial\Omega_{\boldsymbol{u}}),$

where
$\partial\Omega=\overline{\partial\Omega_{p}}\cup\overline{\partial\Omega_{\boldsymbol{u}}}$,
$\overline{\partial\Omega_{p}}\cap\overline{\partial\Omega_{\boldsymbol{u}}}=\varnothing$
and $\boldsymbol{\hat{n}}$ is the outward unit normal vector. Moreover, we
assume $f\in L^{2}(\Omega)$.

Our numerical approximation is derived from the weak formulation of the above
velocity-pressure problem. In order to introduce the weak formulation, we
first define the following spaces

$\displaystyle W(\Omega)$ $\displaystyle=L^{2}(\Omega),$ $\displaystyle
H(\div;\Omega)$ $\displaystyle=\\{\boldsymbol{v}\in(L^{2}(\Omega))^{2}\ |\
\nabla\cdot{{\boldsymbol{v}}}\in L^{2}(\Omega)\\},$

and the set

$\displaystyle V_{g_{\boldsymbol{u}}}(\Omega)$
$\displaystyle=\\{\boldsymbol{v}\in H(\div;\Omega)\ |\
\boldsymbol{v}\cdot\boldsymbol{\hat{n}}=g_{\boldsymbol{u}}\ \hbox{on}\
\partial\Omega_{\boldsymbol{u}}\\},$

for some function $g_{\boldsymbol{u}}$. The global weak form of the pressure-
velocity system (1) is given by finding $\\{p,\boldsymbol{u}\\}\in W\times
V_{g_{\boldsymbol{u}}}$ such that

$\displaystyle(\nabla\cdot\boldsymbol{u},w)_{\Omega}=(f,w)_{\Omega},~{}~{}\forall
w\in W,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$ (3)
$\displaystyle(k^{-1}({\boldsymbol{x}})\boldsymbol{u},{\boldsymbol{v}})_{\Omega}-(p,\nabla\cdot{\boldsymbol{v}})_{\Omega}+\langle
g_{b},{\boldsymbol{v}}\cdot\boldsymbol{\hat{n}}\rangle_{\partial\Omega_{p}}=0,~{}~{}\forall{\boldsymbol{v}}\in
V_{0},$ (4)

where $(\cdot,\cdot)_{\Omega}$ is the $L^{2}(\Omega)$ inner product and
$\langle\cdot,\cdot\rangle_{\partial\Omega}$ is the $L^{2}(\partial\Omega)$
inner product involving line integration over $\partial\Omega$.

The system (3-4) is approximated by the lowest order Raviart-Thomas space [24,
25] that is equivalent to cell-centered finite differences for a uniform
partition of $\Omega$. The resulting problem for the pressure variable is
symmetric positive definite. Thus, it can be efficiently solved by a
preconditioned gradient method [26]. We use the algebraic multigrid method as
a preconditioner and our elliptic solver has been developed to run on GPUs
[15, 27].

##  3 Subsurface Characterization

###  3.1 The Bayesian Framework

Our focus in this work is in the characterization of the permeability field
conditioned on pressure data. The available pressure data comes in the form of
a red-black chessboard pattern: we assume that pressure measurements are
available at all black cells (this type of problem has been investigated in
[4] and references therein). The (log of the) permeability field is denoted by
$\boldsymbol{\eta}$ and $R_{p}$ refers to the reference pressure data. A
Bayesian statistical approach consisting of a preconditioned MCMC method
combined with a novel multiscale sampling strategy is used to solve the
inverse problem for the permeability field. The posterior probability
conditioned on the pressure data $R_{p}$ is given by Bayes’ rule:

$P(\boldsymbol{\eta}|R_{p})\propto
P(R_{p}|\boldsymbol{\eta})P(\boldsymbol{\eta}),$ (5)

where $P(\boldsymbol{\eta})$ denotes a prior distribution. The normalizing
constant is not needed in an iterative search within MCMC methods. The (log
of) permeability field $\boldsymbol{\eta}(\boldsymbol{\theta})$ is built by
using a local KLE strategy and the chain $\boldsymbol{\theta}$ is evolved by
an MCMC method. A Gaussian likelihood function is assumed (as in [11]), and it
is given by

$P(R_{p}|\boldsymbol{\eta})\propto\exp\Big{(}-(R_{p}-R_{\boldsymbol{\eta}})^{\top}\Sigma(R_{p}-R_{\boldsymbol{\eta}})\Big{)},$
(6)

where $R_{\boldsymbol{\eta}}$ refers to the simulated pressure data. We set
the covariance matrix $\Sigma$ to be $\Sigma=\boldsymbol{I}/2\sigma_{R}^{2}$,
where $\boldsymbol{I}$ and $\sigma_{R}^{2}$ refer to the identity matrix and
the precision parameter, respectively.

An MCMC algorithm is used to sample from the posterior distribution (5). In
the MCMC algorithm an instrumental distribution
$I(\boldsymbol{\eta}_{p}|\boldsymbol{\eta})$ is used to propose a sample
$\boldsymbol{\eta}_{p}=\boldsymbol{\eta}(\boldsymbol{\theta}_{p})$ at each
iteration, where $\boldsymbol{\eta}$ denotes the previously accepted sample.
For a given permeability field, the system (3)-(4) is solved numerically to
give $R_{\boldsymbol{\eta}}$, and the original Metropolis-Hastings [1, 2, 3]
acceptance probability of a proposed sample is given by

${\alpha}(\boldsymbol{\eta},\boldsymbol{\eta}_{p})=\text{min}\left(1,\frac{I(\boldsymbol{\eta}|\boldsymbol{\eta}_{p})P(\boldsymbol{\eta}_{p}|R_{p})}{I(\boldsymbol{\eta}_{p}|\boldsymbol{\eta})P(\boldsymbol{\eta}|R_{p})}\right).$
(7)

In this work we consider a preconditioned MCMC method that will be discussed
in Section 4.2.

###  3.2 Dimensional Reduction

We use a Bayesian statistical framework along with MCMC methods where our
numerical simulator requires a permeability value in each cell of a partition
of the domain of interest. Therefore, we need to generate a large number of
random permeability values (based on the grid size) in each iteration that is
infeasible from a practical point of view. Thus, we need to reduce the
dimension of the uncertainty parameter space describing the permeability
field. We use KLE [20, 28] to achieve the desired dimensional reduction of the
parameter space. Next, we briefly discuss the KLE.

We consider $\log\left[k(\boldsymbol{x})\right]=Y^{k}(\boldsymbol{x})$ to be a
sample of a Gaussian field, where $k(\boldsymbol{x})$ represents the
permeability field and $\boldsymbol{x}$ is a point in the domain $\Omega$. We
also consider $Y^{k}(\boldsymbol{x})\in L^{2}(\Omega)$ with unit probability,
i.e., $Y^{k}(\boldsymbol{x})$ is a second-order stochastic process. If we
assume $E[(Y^{k})^{2}]=0$, then, the permeability field
$Y^{k}(\boldsymbol{x})$ can be written for a given orthonormal basis
$\left\\{\varphi_{i}\right\\}$ of $L^{2}(\Omega)$ as follows:

$Y^{k}(\boldsymbol{x})=\sum_{i=1}^{\infty}Y_{i}^{k},\varphi_{i}(\boldsymbol{x}),\hskip
8.5359pt$ (8)

where $Y_{i}^{k}=\int_{\Omega}Y^{k}(\boldsymbol{x})$ are random coefficients,
and $\varphi_{i}(\boldsymbol{x})$ are eigenfunctions with the corresponding
eigenvalues $\lambda_{i}=E[(Y_{i}^{k})^{2}]>0$. The pairs
$(\lambda_{i},\varphi_{i}(\boldsymbol{x}))$ satisfy the integral equation

$\int_{\Omega}R(\boldsymbol{x}_{1},\boldsymbol{x}_{2})\varphi_{i}(\boldsymbol{x}_{2})d\boldsymbol{x}_{2}=\lambda_{i}\varphi_{i}(\boldsymbol{x}_{1}),~{}~{}i=1,2,...$
(9)

for a given covariance function $R(\boldsymbol{x}_{1},\boldsymbol{x}_{2})$.
Setting $\theta_{i}^{k}=Y_{i}^{k}/\sqrt{\lambda_{i}}$ in Eq. (8) we have

$Y^{k}(\boldsymbol{x})=\sum_{i=1}^{\infty}\sqrt{\lambda_{i}}{\theta}_{i}^{k}\varphi_{i}(\boldsymbol{x}),$
(10)

where the eigenfunctions $\varphi_{i}$ and the corresponding eigenvalues
$\lambda_{i}$ satisfy Eq. (9). The eigenvalues are assumed to be arranged in
descending order. The series shown in the Eq. (10) is known as the Karhunen-
Loève expansion. The first $N$ dominating eigenvalues are considered in the
KLE so that the energy $E$ is above $95\%$ [29], i.e.,

$E=\dfrac{\displaystyle\sum_{i=1}^{N}\lambda_{i}}{\displaystyle\sum_{i=1}^{\infty}\lambda_{i}}\geq
95\%.$ (11)

We, thus, can define the truncated KLE by

$Y^{k}_{N}(\boldsymbol{x})=\displaystyle\sum_{i=1}^{N}\sqrt{\lambda_{i}}\theta_{i}^{k}\varphi_{i}(\boldsymbol{x}).$
(12)

###  3.3 Convergence Assessment of MCMCs

We consider a problem that consists of sampling the permeability field
conditioned on pressure measurements. We use a Bayesian statistical approach
(discussed in Section 3.1) along with a preconditioned MCMC method [11, 12] to
characterize the permeability field of our domain of interest. Two critical
issues, namely, where to begin (burn-in) and when to terminate (convergence),
need to be addressed when MCMC methods are used. We now discuss the
convergence of MCMC methods that we use in our investigation to construct one
of the rock properties (permeability field).

A number of convergence criteria [30, 31, 32] for MCMCs have been developed
with a solid theoretical foundation. Several review papers, where authors used
MCMC convergence diagnostics, are available in the literature [33, 34, 35,
36]. Note that in [34] the authors discussed thirteen MCMC convergence
diagnostics. The convergence diagnostics described in [36] are now widely
used. In this work we use two popular diagnostic tools, namely, the Potential
Scale Reduction Factor (PSRF) and the multivariate PSRF (MPSRF), to diagnose
the convergence of MCMC algorithms. Between these two, the MPSRF method takes
all the parameters into account for accessing convergence of the MCMC methods.
Thus, the MPSRF is more restrictive than the PSRF.

The PSRF and MPSRF measures rely on multiple chains. Thus we are required to
run $m>1$ independent chains in parallel with different initial points drawn
from an overdispersed distribution. The effect of starting at different
initial points is made minimal by discarding the first few iterations as burn-
in. Let us denote by $\boldsymbol{\theta}$ an $N$-dimensional parameter
vector, and let $l$ represents the number of posterior draws for each of the
$m$ chains. Furthermore, assume that $\boldsymbol{\theta}_{j}^{c}$ denote the
value of the parameter vector $\boldsymbol{\theta}$ generated at iteration $c$
in $j$th chain of the MCMC algorithm. The posterior variance-covariance matrix
is then estimated as

$\mathbf{\widehat{V}}=\frac{l-1}{l}\mathbf{W}+\left(1+\frac{1}{m}\right)\frac{\mathbf{B}}{l}.$
(13)

The within- and between-sequence (chain) covariance matrix $\mathbf{W}$ and
$\mathbf{B}$ are calculated as

$\mathbf{W}=\frac{1}{m(l-1)}\sum\limits_{j=1}^{m}\sum\limits_{c=1}^{l}\left(\boldsymbol{\theta}_{j}^{c}-\boldsymbol{\bar{\theta}}_{j.}\right)\left(\boldsymbol{\theta}_{j}^{c}-\boldsymbol{\bar{\theta}}_{j.}\right)^{T},$
(14)

and

$\mathbf{B}=\frac{l}{m-1}\sum\limits_{j=1}^{m}\left(\boldsymbol{\bar{\theta}}_{j.}-\boldsymbol{\bar{\theta}}_{..}\right)\left(\boldsymbol{\bar{\theta}}_{j.}-\boldsymbol{\bar{\theta}}_{..}\right)^{T},$
(15)

respectively. $\boldsymbol{\bar{\theta}}_{j.}$ denote within chain mean and
$\boldsymbol{\bar{\theta}}_{..}$ represent the mean between $m$ combined
chains, respectively. $T$ denotes the transpose of a matrix. The PSRFs are
calculated using the two estimators $\mathbf{\widehat{V}}$ and $\mathbf{W}$
defined by

$\displaystyle\text{PSRF}_{\text{i}}=\sqrt{\frac{\text{diag}(\mathbf{\widehat{V}})_{i}}{\text{diag}(\mathbf{W})_{i}}},~{}~{}~{}~{}\text{where}\,\,\,i=1,2,...,N.$
(16)

A large PSRFi suggests that either the estimate of the between variance can be
decreased by taking more samples into account or by taking further samples one
could increase the within variance. It indicates that the simulated sequences
have not yet traversed the parameter space completely. On the other hand, if
the maximum of PSRF values is close to $1$, we can draw the conclusion that
each of the $m$ chains of $l$ simulated samples is close to the target
distribution. The MPSRF is estimated by using the maximum root statistic. As
in [36] it is defined by

MPSRF $\displaystyle=\sqrt{\max_{\mathbb{a}}\vskip 3.0pt plus 1.0pt minus
1.0pt\frac{\mathbb{a}^{T}\mathbb{\widehat{V}}\mathbb{a}}{\mathbb{a}^{T}\mathbb{W}\mathbb{a}}}$
$\displaystyle=\sqrt{\max_{\mathbb{a}}\vskip 3.0pt plus 1.0pt minus
1.0pt\frac{\mathbb{a}^{T}\left[\frac{l-1}{l}\mathbb{W}+\left(1+\frac{1}{m}\right)\frac{\mathbb{B}}{l}\right]\mathbb{a}}{\mathbb{a}^{T}\mathbb{W}\mathbb{a}}}$
$\displaystyle=\sqrt{\frac{l-1}{l}+\left(\frac{m+1}{m}\right)\max_{\mathbb{a}}\frac{\mathbb{a}^{T}\frac{\mathbb{B}}{l}\mathbb{a}}{\mathbb{a}^{T}\mathbb{W}\mathbb{a}}}$
$\displaystyle=\sqrt{\frac{l-1}{l}+\left(\frac{m+1}{m}\right)\lambda},$

where $\mathbb{a}\in\mathbb{R}^{N}$ is an arbitrary vector, and $\lambda$ is
the greatest eigenvalue of the positive definite matrix
$\mathbb{W}^{-1}\mathbb{B}/l$. If the means of between chains are equal, the
between chain covariance matrix $\mathbf{B}$ becomes zero. In this case, the
chains mix well and $\lambda\rightarrow 0$. Thus, as the MPSRF approaches to
$1$, it guarantees a convergence for sufficiently large sample size.

##  4 Multiscale Sampling

###  4.1 The Multiscale Prior Distribution

We begin with the description of a decomposition of the domain $\Omega$. Our
multiscale sampling strategy is based on two non-overlapping partitions of the
domain $\Omega$: the first is a uniform fine Cartesian mesh $\Omega^{f}$ where
the values of the absolute permeability field are piecewise constant. This is
also the mesh used for the numerical solution of the system (3)-(4). The
second is a coarse Cartesian mesh $\Omega^{c}$ constructed as sets of elements
in $\Omega^{f}$ (see Figure 1) where a KLE will be applied for local
dimensional reduction. The proposed method is based on partitions
$\Omega^{\gamma}$ into rectangles $\\{\Omega_{i}^{\gamma},\
i=1,\dots,M_{\gamma}\\}$ (see Figure 1), such that

$\bar{\Omega^{\gamma}}=\bigcup^{M_{\gamma}}_{i=1}\bar{\Omega_{i}^{\gamma}};\quad\Omega_{i}^{\gamma}\cap\Omega_{k}^{\gamma}=\varnothing,\quad
i\neq k,\quad\gamma=c,f.$

Define $\Gamma=\partial\Omega$ and, for $i=1,\dots,M_{\gamma}$:

$\quad\Gamma_{ik}^{\gamma}=\Gamma_{ki}^{\gamma}=\partial\Omega_{i}^{\gamma}\cap\partial\Omega_{k}^{\gamma},\quad\gamma=c,f.$

For each element of the coarse partition $\\{\Omega_{i}^{\gamma},\
i=1,\dots,M_{c}\\}$, we define the set

$\mathcal{S}_{i}=\\{j:\Omega_{j}^{f}\subset\Omega_{i}^{c}\\}.$

As indicated in Figure 1, we refer to two length scales in the description of
the new multiscale procedure: $H$, the mesh size for the coarse partition and
$h$, the mesh size of an underlying fine grid.

Figure 1: The fine $\Omega^{f}$ (dashed lines) and coarse $\Omega^{c}$ (solid
lines) partitions of $\O$ along with the three spatial scales used in the
definition of the new multiscale procedure.

We will consider a blocking strategy [36] for Gibbs sampling within a
Metropolis-Hastings algorithm. In order to define it we decompose the
$\boldsymbol{\theta}$ vector in Eq. (12) into orthogonal subspaces
corresponding to blocks with the same number of components, that are denoted
by $\boldsymbol{\theta}^{i}$, for $i=1,\dots,M_{c}$. Each block of thetas is
used to generate a local Gaussian field within its corresponding subdomain, as
illustrated in Fig. 2. The update of each $\boldsymbol{\theta}^{i}$ block is
based on the random walk sampler (RWS) of [37]. It is given, for
$i=1,\dots,M_{c}$, by

${\boldsymbol{\theta}_{p}}^{i}=\sqrt{1-\beta^{2}}\,{\boldsymbol{\theta}}^{i}+\beta\,{\boldsymbol{\epsilon}}^{i},$
(17)

where the current sample is denoted by ${\boldsymbol{\theta}_{p}}^{i}$ and the
previously accepted sample by $\boldsymbol{\theta}^{i}$. The algorithmic
parameter $\beta$ is used for tuning the sampler and
${\boldsymbol{\epsilon}}^{i}$ represents a $\mathcal{N}(0,1)$-random vector.
Not all components of ${\boldsymbol{\theta}_{p}}^{i}$ are updated
simultaneously. Blocking is used again so that only a subset of the components
is modified in one MCMC iteration. We view ${\boldsymbol{\theta}_{p}}^{i}$ as
a column vector and we are going to refer to the local blocking number as the
number of contiguous components of this column vector that are updated
simultaneously. For the purpose of sampling, we consider the
${\boldsymbol{\theta}_{p}}^{i}$ components ordered by their corresponding
eigenvalues in the local KLE, from the largest towards the smallest one.

The samples produced by the local sampling strategy discussed above produces
Gaussian samples that show discontinuities in $\Gamma_{ik}^{c}$. Motivated by
downscaling strategies developed for multiscale methods (that aim at removing
flux discontinuities at subdomain boundaries [23]) in order to complete the
construction of one sample from our multiscale prior distribution an averaging
method is used to condition each sample on the available data at nearest
neighbor subdomains. The averaging procedure is illustrated in Figure 2.

Figure 2: The mapping between blocks of theta variables and elements of the
coarse $\Omega^{c}$ (solid lines) partition of $\Omega$. Each block of thetas
is used to generate a local Gaussian field within its corresponding subdomain.

A length scale $\overline{H}$ is set (a fraction of the correlation length
that enters in the construction of the prior distribution) and, for each $i$,
all cells of $\mathcal{S}_{i}$ that are at a distance of $\overline{H}$ (or
less) to $\Gamma_{ik}^{c}$ have their current values replaced by local
averages (that preserve both their mean value and variance - if they were
uncorrelated). Figure 3 illustrates a sample before and after this averaging
procedure. In this figure, $H=0.25$ and the averaging is applied on the
boundary of the top right subdomain. Note that if the correlation lengths are
not equal, the circle for the averaging in Figure 2 should be replaced by an
ellipse.

In conclusion, the multiscale prior distribution requires three user-specified
parameters:

  * 1.

The value of $H$: the subdomain size;

  * 2.

The value of $\overline{H}$: the length scale local averages are taken;

  * 3.

The local blocking number.

Figure 3: Sample of a permeability field before (left) and after (right) a
local averaging is applied at the top right subdomain of a domain
decomposition with $H=0.25$.

We remark that the number of blocks for the localized Gibbs sampling in each
subdomain is given by $N_{\text{local}}=N/(M_{c}*N_{lb})$, where $N_{lb}$
denotes the local blocking number. We refer to these blocks as $B_{k}$,
$k=1,\dots,N_{\text{local}}$, and we also define the local stochastic
dimension to be $N_{c}=N/M_{c}$. We conjecture that improved convergence for
MCMCs should be observed for $H>L$, where $L$ denotes the largest correlation
length in the definition of the prior distribution. Further studies are needed
to check the validity of this conjecture.

###  4.2 The Multiscale Sampling Method

We now provide a detailed algorithm of the new Multiscale Sampling Method
(MSM) that consists of a preconditioned MCMC with a multiscale prior
distribution. If the number of subdomains $M_{c}=1$, then the proposed
algorithm reduces to the classical preconditioned MCMC [11, 12] with Gibbs
sampling associated with the local blocking number.

We first discuss the algorithm of the preconditioned MCMC method. The
filtering step of this method is based on a coarse-scale model approximation
of the governing system (3)-(4). The coarse-scale discretization is similar to
the fine-scale discretization and the permeability field
$\boldsymbol{\eta}(\boldsymbol{\theta})$ is projected on the coarse-scale. An
upscaling procedure [38] is used to set an effective permeability field that
provides a similar average response as that of the underlying fine-scale
problem. The numerical simulator is run on the coarse-scale model and produces
the coarse-grid pressure field $R_{c}$. The coarse-scale and fine-scale
acceptance probabilities are estimated as

$\displaystyle{\alpha}_{c}(\boldsymbol{\eta},\boldsymbol{\eta}_{p})$
$\displaystyle=\text{min}\left(1,\dfrac{I(\boldsymbol{\eta}|\boldsymbol{\eta}_{p})P_{c}(\boldsymbol{\eta}_{p}|R_{p})}{I(\boldsymbol{\eta}_{p}|\boldsymbol{\eta})P_{c}(\boldsymbol{\eta}|R_{p})}\right),\text{and}$
(18) $\displaystyle{\alpha}_{f}(\boldsymbol{\eta},\boldsymbol{\eta}_{p})$
$\displaystyle=\text{min}\left(1,\dfrac{P_{f}(\boldsymbol{\eta}_{p}|R_{p})P_{c}(\boldsymbol{\eta}|R_{p})}{P_{f}(\boldsymbol{\eta}|R_{p})P_{c}(\boldsymbol{\eta}_{p}|R_{p})}\right),$

where $P_{c}$ and $P_{f}$ are the posterior probabilities calculated at
coarse- and fine-scale, respectively. In MSM we first construct a local
permeability field $\boldsymbol{\eta}(\boldsymbol{\theta}^{i})$ for each
subdomain $\Omega^{c}_{i},i=1,\dots M_{c}$ using Eq. (12). To do so, we
generate local KLE data in subdomains $\Omega^{c}_{i},i=1,\dots M_{c}$ with
size $\dfrac{1}{2^{n}}$, $n=1,2,\dots$. Then we construct the global
permeability field by taking local averages. The MSM algorithm is presented in
Algorithm 1.

Algorithm 1 The Multiscale Sampling Method (MSM)

1:For a given covariance function $R$ solve Eq. (9) to get a KLE in Eq. (12),
which is used in all the subdomains, $\Omega_{i}^{c},i=1,\dots,M_{c}$.

2:for $j=1$ to $M_{\text{mcmc}}$ do

3: for $i=1$ to $M_{c}$ do

4: for $k=1$ to $N_{\text{local}}$ do

5: Generate i.i.d., $\mathcal{N}(0,1)$ Gaussian variables to construct
$\boldsymbol{\theta}_{p}$ using Eq. (17) for block $B_{k}$ in
$\Omega_{i}^{c}$.

6: Construct a local Gaussian sample (in each subdomain) using the KLE to set
a preliminary value for the Gaussian sample at the $\mathcal{S}_{i}$ cells.

7: Run the local averaging algorithm to remove discontinuities.

8: Compute the upscaled permeability on the coarse-scale using
$\boldsymbol{\eta}_{p}$.

9: Solve the forward problem on the coarse-scale to get $R_{c}$.

10: Compute the coarse-scale acceptance probability
${\alpha}_{c}(\boldsymbol{\eta},\boldsymbol{\eta}_{p})$.

11: if $\boldsymbol{\eta}_{p}$ is accepted then

12: Use $\boldsymbol{\eta}_{p}$ in the fine-scale simulation to get $R_{f}$.

13: Compute the fine-scale acceptance probability
${\alpha}_{f}(\boldsymbol{\eta},\boldsymbol{\eta}_{p})$.

14: if $\boldsymbol{\eta}_{p}$ is accepted then
$\boldsymbol{\eta}=\boldsymbol{\eta}_{p}$.

15: end if

16: end if

17: $j=j+1.$

18: end for

19: end for

20:end for

##  5 Kriging and Conditioning

In this section, we combine the multiscale sampling method with conditioning
by projection for the sampling. The conditioning by projection method has been
discussed in detail in [39]. This method consists of two steps. In the first
step, for given permeability values at sparse locations in the domain,
generate a kriged field for the domain. In the second step, project the i.i.d
$\mathcal{N}(0,1)$-random vector onto the nullspace of a data matrix defined
in terms of KLE to calculate the linear combination in Eq. (12). The final
permeability field is obtained by adding the fields defined in both steps
(after taking exponential). In the following subsections we briefly discuss
the kriging interpolation and projection method for conditioning.

###  5.1 Kriging Interpolation

Kriging is an interpolation method that is derived from a regionalized
variable theory [40, 41]. It employs a limited set of sampled data points to
compute the value of a variable throughout a continuous spatial field. Kriging
uses the spatial correlation between sampled points to interpolate the values
in the spatial field. This interpolation gives the exact values of the field
at the known locations.

###  5.2 The Projection Method for Conditioning

Following the discussion in [39], we need to extract a data matrix that is
defined in terms of KLE. We assume the Gaussian field $Y(\boldsymbol{x})$
defined in Eq. (12) is a Gaussian perturbation on top of a kriged field
$\widehat{Y}(\boldsymbol{x})$, thus we can write

$\displaystyle
Y(\boldsymbol{x})-\widehat{Y}(\boldsymbol{x})=\sum_{i=1}^{N}\sqrt{\lambda_{i}}\varphi_{i}(\boldsymbol{x})\theta_{i}=\boldsymbol{\phi}^{T}(\boldsymbol{x})\sqrt{D}\boldsymbol{\theta},$
(19)

where, for each
$\boldsymbol{x},\boldsymbol{\phi}(\boldsymbol{x})\in\mathbb{R}^{N}$, and $D$
is a diagonal matrix containing $N$ dominating eigenvalues. If we have $M$
measured data values at sparse locations, then we can define the following
homogeneous linear system of equations

$\displaystyle A\boldsymbol{\theta}=\boldsymbol{0},$

where
$A=\boldsymbol{\phi}^{T}(\boldsymbol{\hat{x}})\sqrt{D}\in\mathbb{R}^{M\times
N}$ is the desired data matrix. Finally, we project the vector
$\boldsymbol{\theta}$ onto the nullspace of the matrix $A$ to get the closest
vector to $\boldsymbol{\theta}$ in the nullspace of the data matrix $A$. i.e.

$\boldsymbol{\hat{\theta}}=P\boldsymbol{\theta},$

where $P$ is a projection matrix [42]. Therefore, we can write

$\displaystyle
Y(\boldsymbol{x})=\widehat{Y}(\boldsymbol{x})+\sum_{i=1}^{N}\sqrt{\lambda_{i}}\varphi_{i}(\boldsymbol{x})\hat{\theta_{i}},\quad\boldsymbol{\hat{\theta}}=(\hat{\theta}_{1},\dots,\hat{\theta}_{N}).$

##  6 Numerical Results

In this section, we describe the simulation study for the problem of interest.
We test the proposed multiscale sampling method in four examples. In each
example, we numerically solve the system containing Eqs. (3)-(4) on the domain
$\Omega=[0,1]\times[0,1]$. We present a comparative study between the
preconditioned MCMC method with and without multiscale sampling in the first
three examples. In the last example, we analyze MSM with and without
conditioning for a problem of higher dimensional stochastic space. In MSM, we
apply KLE to construct a permeability field for each subdomain, and then
construct the global permeability field. In KLE, we use the following
covariance function:

$\displaystyle
R(\boldsymbol{x}_{1},\boldsymbol{x}_{2})=\sigma_{Y}^{2}\,\text{exp}\left(-\frac{|x_{1}-x_{2}|^{2}}{2L_{x}^{2}}-\frac{|y_{1}-y_{2}|^{2}}{2L_{y}^{2}}\right),$
(20)

where $L_{x}$ and $L_{y}$ are the correlation lengths and
$\sigma_{Y}^{2}=\text{Var}[(Y^{k})^{2}]$. We take $\sigma_{Y}^{2}=1$ in all
the four examples. Moreover, we set the source term $f=0$ and impose Dirichlet
boundary conditions, $p=1$ and $p=0$, on the left and right boundaries,
respectively. We also set a no-flow (Neumann-type boundary condition)
condition on the other two boundaries. We run four MCMC chains for each
method. In order to remove the discontinuities between subdomains in our
numerical studies, we set the length scale $\overline{H}$ to be
$\overline{H}=\min\\{\frac{L_{x}}{2},\frac{L_{y}}{2}\\}$. Below we discuss the
numerical results.

###  6.1 Example 1

In the first example, we consider $L_{x}=L_{y}=0.2$ in Eq. (20). We then
generate KLEs for the global and MSM $2\times 2$ samplings. In MSM $2\times 2$
sampling, we use $H=0.5$. Figure 4 illustrates the decay of the eigenvalues
(in log scale) for both samplings. Note that the relationship between the
eigenvalues in the global sampling and the eigenvalues in multiscale sampling
can be obtained directly by a change of variables in Eq. (9). We take the
first $20$ eigenvalues that preserve more than $97\%$ of the total energy for
the global sampling. Five eigenvalues are used for each subdomain in the
multiscale sampling. We generate a reference synthetic permeability field on a
computational fine mesh of size $16\times 16$, and then run our numerical
simulator to generate the corresponding reference pressure field. We run the
MCMC algorithms conditioned on this pressure field. Figure 5 shows these
reference fields. Furthermore, we use a coarse mesh of size $8\times 8$ as a
filtering step in the preconditioned MCMC. We let the local blocking number
$N_{lb}=1$ and $\beta=0.5$ in Eq. (17).

Figure 4: Decay of eigenvalues for the global and multiscale samplings in the
first example.

Figure 5: Reference log permeability field (left) and the corresponding
reference pressure field (right) for the first example.

As we discussed in subsection 3.3, we analyze the convergence of the MCMCs
using PSRFs and MPSRF. An MCMC method converges to the stationary distribution
if both MPSRF and the maximum of PSRFs get closer to 1. In [43] the author
considered a value of $1.2$ for these parameters to confirm the convergence of
the chains. In line with that, we decide to stop the simulation once these
parameters reach $1.2$. Figure 6 shows that the preconditioned MCMC methods
with and without multiscale sampling converge. However, the plots at the
bottom in Figure 6 show that the preconditioned MCMC method with multiscale
sampling converges to the stationary distribution earlier than the method
without multiscale sampling. Table 1 shows the acceptance rates for both
methods as well as the precision parameters for coarse- and fine-grid
simulations. The acceptance rate increases slightly when we use MSM. The
errors between the reference and simulated pressure data, which are used in
the likelihood function, for both methods are shown in Figure 7. Both methods
produce similar error curves.

Figure 6: Top: The maximum of PSRFs and MPSRF for the MCMC method with and without multiscale sampling in the first example. Bottom: Tails of the maximum of PSRFs and MPSRF curves. Table 1: A comparison of acceptance rates for the MCMC with and without MSM for the first example. |  MCMC with global sampling |  MCMC with multiscale sampling |
---|---|---|---
$\sigma_{F}^{2}$ | $10^{-3}$ | $10^{-3}$ |
$\sigma_{C}^{2}$ | $5\times 10^{-3}$ | $5\times 10^{-3}$ |
acc. rate | $53\%$ | $55\%$ |
Figure 7: Error curves of the preconditioned MCMC with and without multiscale
sampling for the first example.

After the convergence of both MCMC methods, we take 10000 log permeability
values from each chain and draw the posterior histograms for three cells with
high, medium and low permeability values in the computational domain. See
Figure 5 for those three cells. Figure 8 shows the posterior histograms with
the true (red vertical line) and mean (green vertical line) values of the log
permeability for the cells.

---
Cell 1

---
Cell 2

---
Cell 3

Figure 8: Posterior histograms for three selected cells.

Table 2 shows these values and the corresponding standard deviation. We
observe in Figure 8 that when we use the multiscale sampling method, the mean
values are within one standard deviation (green horizontal line) for all the
three cells. Also, in MSM, the mean of the posterior histogram is almost the
same as the true value for the cell 3. We do not observe a similar behavior in
the posterior histograms in the global sampling method.

Table 2: True and mean values of log permeability with the corresponding
standard deviation for three cells.

|  Cell 1 |  Cell 2 |  Cell 3
---|---|---|---
|  Global |  MSM |  Global |  MSM |  Global |  MSM
True |  $-1.54$ |  $-$ |  $-0.3$ |  $-$ |  $1.3$ |  $-$
Mean |  $-1.08$ |  $-1.16$ |  $-0.69$ |  $-0.07$ |  $0.78$ |  $1.35$
SD |  $0.36$ |  $0.34$ |  $0.39$ |  $0.27$ |  $0.42$ |  $0.26$

We now compare the reference field with some of the simulated permeability
fields from two selected chains. See Figures 9 and 10. Other chains also show
a similar behavior. Although both MCMC methods converged, we observe that at
iteration 60000 the permeability field obtained from the multiscale sampling
is closer to the reference permeability field than the field obtained without
the multiscale sampling.

Figure 9: First row: Reference log permeability filed. Second row: Accepted
permeability fields in the global sampling method. Third row: Accepted
permeability fields in MSM $2\times 2$. From left to right, log permeability
fields at 20000, 40000 and 60000 iterations, respectively, from chain 1 in the
first example.

Figure 10: First row: Reference log permeability filed. Second row: Accepted
permeability fields in the global sampling method. Third row: Accepted
permeability fields in MSM $2\times 2$. From left to right, log permeability
fields at 20000, 40000 and 60000 iterations, respectively, from chain 2 in the
first example.

###  6.2 Example 2

In the second example we consider the case where the correlation lengths are
not equal, i.e., $L_{x}=0.2$ and $L_{y}=0.06$ in Eq. (20). Figure 11 shows the
decay of the eigenvalues (in log scale) for the global sampling, MSM $2\times
2$, and MSM $4\times 4$ methods.

Figure 11: Decay of eigenvalues for the global and multiscale sampling in the
second example.

In KLE, we consider the first $64$ eigenvalues, which preserve $97.87\%$ of
the total energy, in the global sampling method. For MSM $2\times 2$ and MSM
$4\times 4$, we take $16$ and $4$ eigenvalues, respectively. We generate a
reference synthetic permeability field on a computational grid of size
$32\times 32$ and then run the numerical simulator to generate the
corresponding reference pressure field. See Figure 12 for the reference
fields.

Figure 12: Reference log permeability field (left) and the corresponding
reference pressure field (right) for the second example.

We use the same coarse mesh of size $8\times 8$ as in the first example. We
also use the same local blocking number $N_{lb}=1$. We set the tuning
parameter $\beta=0.75$ in Eq. (17).

Let us consider the convergence analysis of these methods. We take $170000$
proposals from each chain to compute the MPSRF and the maximum of PSRFs.
Figure 13 shows the maximum of PSRFs and MPSRF curves. At the tails of the
maximum of PSRFs and MPSRF curves, we have the values $1.2$ and $1.4$,
respectively, for MSM $4\times 4$. These values are slightly higher for MSM
$2\times 2$. Thus, we can conclude that MSM $4\times 4$ converges to the
stationary distribution faster than MSM $2\times 2$. On the other hand, the
PSRF and MPSRF curves in the global sampling method do not show any sign of
converging at the same number of iterations. Moreover, the acceptance rates
are better for the multiscale sampling methods. See Table 3. The error curves
for this study are shown in Figure 14. They are very comparable. Figures 15
and 16 present simulated permeability fields from two selected chains. From
these figures, we observe that both MSM $2\times 2$ and MSM $4\times 4$
recover the permeability fields better than the global sampling method.

Figure 13: The maximum of PSRFs and MPSRF for the MCMC method with and without multiscale sampling for the second example. Table 3: A comparison of acceptance rates for the MCMC with and without MSM in the second example. |  MCMC global |  MCMC with MSM $2\times 2$ |  MCMC with MSM $4\times 4$
---|---|---|---
$\sigma_{F}^{2}$ | $10^{-3}$ | $10^{-3}$ | $10^{-3}$
$\sigma_{C}^{2}$ | $5\times 10^{-3}$ | $5\times 10^{-3}$ | $5\times 10^{-3}$
acc. rate | $50\%$ | $54\%$ | $55\%$
Figure 14: Error curves of the preconditioned MCMC with and without multiscale
sampling for the second example.

Figure 15: First row: Reference log permeability filed. Second row: Accepted
permeability fields in the global sampling method. Third row: Accepted
permeability fields in MSM $2\times 2$. Fourth row: Accepted permeability
fields in MSM $4\times 4$. From left to right, log permeability fields at
20000, 50000 and 100000 iterations, respectively, from chain 1 in the second
example.

Figure 16: First row: Reference log permeability filed. Second row: Accepted
permeability fields in the global sampling method. Third row: Accepted
permeability fields in MSM $2\times 2$. Fourth row: Accepted permeability
fields in MSM $4\times 4$. From left to right, log permeability fields at
20000, 50000 and 100000 iterations, respectively, from chain 2 in the second
example.

###  6.3 Example 3

In this example, we test the proposed method on a large grid size of $64\times
64$ with the correlation lengths, $L_{x}=L_{y}=0.1$ in Eq. (20). Figure 4
shows the decay of the eigenvalues for the methods, the global sampling, MSM
$2\times 2$, and MSM $4\times 4$. We consider $64$ eigenvalues, which preserve
$95.8\%$ of the total energy in KLE, in the global sampling. The numbers of
eigenvalues for MSM $2\times 2$ and MSM $4\times 4$ are 16 and 4,
respectively. The synthetic reference permeability field is generated on a
computational mesh of size $64\times 64$. Then, the numerical simulator is
used to generate the corresponding reference pressure field. Figure 18 shows
the reference permeability field and the corresponding pressure distribution
on the grid. We use a coarse mesh of size $16\times 16$ in the filtering step
in the preconditioned MCMC. Let $N_{lb}=2$. We set $\beta=0.2$ in Eq. (17).

Figure 17: Decay of eigenvalues for the global and multiscale sampling for the
third example.

Figure 18: Reference log permeability field (left) and the corresponding
reference pressure field (right) for the third example.

Let us consider PSRFs and MPSRF curves for these methods. We take $240000$
proposals from each chain in constructing the PSRF and MPSRF curves. We show
the maximum of PSRFs and MPSRF curves in Figure 19. For MSM $4\times 4$, the
values at the tails of the PSRF and MPSRF curves are $1.2$ and $1.6$,
respectively. These values indicate that the curves in MSM $4\times 4$ are
closer to the convergence. However, the curves in MSM $2\times 2$ and global
sampling method are very far from reaching a convergence. Table 4 shows that
the acceptance rates for the multiscale sampling methods are also slightly
better than that of the global sampling method. The error curves are
comparable for these methods. See Figure 20. Figures 21 and 22 compare the
accepted permeability fields for two selected chains in the MCMC simulation.
Both MSM $2\times 2$ and MSM $4\times 4$ recover the fields better than the
global sampling method.

Figure 19: The maximum of PSRFs and MPSRF for the MCMC method with and without multiscale sampling for the third example. Table 4: A comparison of acceptance rates for the MCMC with and without MSM for the third example. |  MCMC global |  MCMC with MSM $2\times 2$ |  MCMC with MSM $4\times 4$
---|---|---|---
$\sigma_{F}^{2}$ | $10^{-3}$ | $10^{-3}$ | $10^{-3}$
$\sigma_{C}^{2}$ | $5\times 10^{-3}$ | $5\times 10^{-3}$ | $5\times 10^{-3}$
acc. rate | $41\%$ | $43\%$ | $43\%$
Figure 20: Error curves of the preconditioned MCMC with and without multiscale
sampling for the third example.

Figure 21: First row: Reference log permeability filed. Second row: Accepted
permeability fields in the global sampling method. Third row: Accepted
permeability fields in MSM $2\times 2$. Fourth row: Accepted permeability
fields in MSM $4\times 4$. From left to right, log permeability fields at
80000, 160000 and 240000 iterations, respectively, from chain 1 in the third
example.

Figure 22: First row: Reference log permeability filed. Second row: Accepted
permeability fields in the global sampling method. Third row: Accepted
permeability fields in MSM $2\times 2$. Fourth row: Accepted permeability
fields in MSM $4\times 4$. From left to right, log permeability fields at
80000, 160000 and 240000 iterations, respectively, from chain 2 in the third
example.

###  6.4 Example 4

In this example we compare MSM $4\times 4$ with and without conditioning for
the problem in example 3. However, in the conditioning approach, we
incorporate the permeability measurements at eight sparse locations in the
field.

Four MCMCs are simulated in each case. We compute the maximum of the PSRFs and
MPSRF by taking $100,000$ samples from each chain (total $400,000$ samples).
In Figure 23 we present the maximum of PSRFs and MPSRF curves. The values in
the tails of the maximum of the PSRFs and the MPSRF are $1.09$ and $1.17$,
respectively, for the method with conditioning. For the method without
conditioning, these values are $2.5$ and $3.5$, respectively. Therefore, we
can say that MSM $4\times 4$ with conditioning reaches convergence earlier
than MSM without conditioning.

Figure 23: The maximum of PSRFs and MPSRF for the multiscale sampling MCMC
method with and without conditioning for the fourth example.

Fig. 24 shows the accepted fields in the method with conditioning for two
chains. The permeability fields were not fully recovered, however, we see a
considerable improvement in the fields in comparison to the fields in Figures
21 and 22, which were obtained using MSM $4\times 4$ without conditioning. We
thus conclude that the conditioning speeds-up the convergence and improves the
characterization in this example

Figure 24:  First row: Reference log permeability filed. Second row: Accepted
permeability fields from chain 1. Third row: Accepted permeability fields from
chain 2. From left to right, log permeability fields at 40000, 60000 and 80000
iterations, respectively, in MSM $4\times 4$ with conditioning.

##  7 Conclusions

We have presented a novel multiscale sampling method aiming at subsurface
characterization. The proposed method is based on a non-overlapping partition
of the domain of the governing partial differential equation that leads to the
localization of the search in the underlying stochastic space. The novel
method is implemented in the framework of a preconditioned Markov Chain Monte
Carlo algorithm.

Through several multi-chain MCMC examples, motivated by subsurface flow
problems, we compare the usual preconditioned Markov Chain Monte Carlo
algorithm with the proposed procedure. Our results show that the new
multiscale sampling method considerably improves the convergence rate of the
preconditioned Markov Chain Monte Carlo algorithm. We also incorporated sparse
measurements of the permeability field in the multiscale sampling method, and
showed that conditioning on this data further improves the convergence of the
proposed method.

The authors and their collaborators are currently applying the method
introduced here to solve the inverse problems associated with single and
multiphase flows in porous media. In these studies the Multiscale Perturbation
Method [44] will be used to speed-up the numerical solution of elliptic
equations in the forward solution of the governing system of equations.
Multiscale sampling procedures based on overlapping domain decompositions are
also being considered.

## Acknowledgments

A. Rahunanthan research was supported by NIFA/USDA through Central State
University Evans-Allen Research Program.

## References

  * Metropolis and Ulam [1949] N. Metropolis, S. Ulam,  The Monte Carlo method,  Journal of the American statistical association 44 (1949) 335–341.
  * Metropolis et al. [1953] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, E. Teller,  Equation of state calculations by fast computing machines,  The journal of chemical physics 21 (1953) 1087–1092.
  * Hastings [1970] W. K. Hastings,  Monte Carlo sampling methods using Markov chains and their applications,  Biometrika 57 (1970) 97–109.
  * Tong et al. [2020] X. T. Tong, M. Morzfeld, Y. M. Marzouk,  MALA-within-Gibbs samplers for high-dimensional distributions with sparse conditional structure,  SIAM Journal on Scientific Computing 42 (2020) A1765–A1788.
  * Brooks et al. [2011] S. Brooks, A. Gelman, G. Jones, X.-L. Meng, Handbook of Markov chain Monte Carlo, CRC press, 2011.
  * Constantine et al. [2016] P. G. Constantine, C. Kent, T. Bui-Thanh,  Accelerating Markov chain Monte Carlo with active subspaces,  SIAM Journal on Scientific Computing 38 (2016) A2779–A2805.
  * Martin et al. [2012] J. Martin, L. C. Wilcox, C. Burstedde, O. Ghattas,  A stochastic Newton MCMC method for large-scale statistical inverse problems with application to seismic inversion,  SIAM Journal on Scientific Computing 34 (2012) A1460–A1487.
  * Moraes et al. [2017] R. Moraes, J. Rodrigues, H. Hajibeygi, J. Jansen, et al.,  Multiscale gradient computation for multiphase flow in porous media,  in: SPE Reservoir Simulation Conference, Society of Petroleum Engineers, 2017, pp. 1168–1193.
  * Jaramillo et al. [2021] A. Jaramillo, R. T. Guiraldello, S. Paz, R. F. Ausas, F. S. Sousa, F. Pereira, G. C. Buscaglia,  Towards HPC simulations of billion-cell reservoirs by multiscale mixed methods,  arXiv preprint arXiv:2103.08513 (2021).
  * Durlofsky [2005] L. J. Durlofsky,  Upscaling and gridding of fine scale geological models for flow simulation,  in: 8th International Forum on Reservoir Simulation Iles Borromees, Stresa, Italy, volume 2024, Citeseer, 2005, pp. 1–59.
  * Efendiev et al. [2006] Y. Efendiev, T. Hou, W. Luo,  Preconditioning Markov Chain Monte Carlo Simulations Using Coarse-Scale Models,  SIAM J. Sci. Comput. 28 (2006) 776–803.
  * Christen and Fox [2005] J. A. Christen, C. Fox,  Markov chain Monte Carlo using an approximation,  Journal of Computational and Graphical Statistics 14 (2005) 795–810.
  * Ginting et al. [2011] V. Ginting, F. Pereira, M. Presho, S. Wo,  Application of the two-stage Markov chain Monte Carlo method for characterization of fractured reservoirs using a surrogate flow model,  Comput. Geosci. 15 (4) (2011) 691–707.
  * Ginting et al. [2015] V. Ginting, F. Pereira, A. Rahunanthan,  Multi-physics Markov chain Monte Carlo methods for subsurface flows,  Mathematics and Computers in Simulation 118 (2015) 224–238.
  * Ginting et al. [2014] V. Ginting, F. Pereira, A. Rahunanthan,  A prefetching technique for prediction of porous media flows,  Computational Geosciences 18 (2014) 661–675.
  * Stuart et al. [2019] G. K. Stuart, S. E. Minkoff, F. Pereira,  A two-stage Markov chain Monte Carlo method for seismic inversion and uncertainty quantification,  Geophysics 84 (2019) R1003–R1020.
  * Lykkegaard et al. [2020] M. B. Lykkegaard, G. Mingas, R. Scheichl, C. Fox, T. J. Dodwell,  Multilevel delayed acceptance MCMC with an adaptive error model in PyMC3,  arXiv preprint arXiv:2012.05668 (2020).
  * Francisco et al. [2014] A. Francisco, V. Ginting, F. Pereira, J. Rigelo,  Design and implementation of a multiscale mixed method based on a nonoverlapping domain decomposition procedure,  Mathematics and Computers in Simulation 99 (2014) 125 – 138.
  * Guiraldello et al. [2018] R. T. Guiraldello, R. F. Ausas, F. S. Sousa, F. Pereira, G. C. Buscaglia,  The multiscale Robin coupled method for flows in porous media,  Journal of Computational Physics 355 (2018) 1–21.
  * Loève [1997] M. Loève, Probability theory, Springer, Berlin, 1997\.
  * Smith and Roberts [1993] A. F. Smith, G. O. Roberts,  Bayesian computation via the Gibbs sampler and related Markov chain Monte Carlo methods,  Journal of the Royal Statistical Society: Series B (Methodological) 55 (1993) 3–23.
  * Chen et al. [2006] Z. Chen, G. Huan, Y. Ma, Computational methods for multiphase flows in porous media, SIAM, Philadelphia, PA, 2006\.
  * Guiraldello et al. [2020] R. T. Guiraldello, R. F. Ausas, F. S. Sousa, F. Pereira, G. C. Buscaglia,  Velocity postprocessing schemes for multiscale mixed methods applied to contaminant transport in subsurface flows,  Computational Geosciences (2020) 1–21.
  * Raviart and Thomas [1977] P.-A. Raviart, J.-M. Thomas,  A mixed finite element method for second order elliptic problems,  in: Mathematical aspects of finite element methods, Springer, 1977, pp. 292–315.
  * Douglas et al. [1997] J. Douglas, F. Furtado, F. Pereira,  On the numerical simulation of waterflooding of heterogeneous petroleum reservoirs,  Computational Geosciences 1 (1997) 155–190.
  * Saad [2003] Y. Saad, Iterative methods for sparse linear systems, SIAM, 2003.
  * Pereira and Rahunanthan [2011] F. Pereira, A. Rahunanthan,  A semi-discrete central scheme for the approximation of two-phase flows in three space dimensions,  Mathematics and Computers in Simulation 81 (2011) 2296–2306.
  * Ginting et al. [2013] V. Ginting, F. Pereira, A. Rahunanthan,  A multi-stage Bayesian prediction framework for subsurface flows,  International J. Uncertain. Quantif. 3 (2013) 499–522.
  * Laloy et al. [2013] E. Laloy, B. Rogiers, J. Vrugt, D. Mallants, D. Jacques,  Efficient posterior exploration of a high-dimensional groundwater model from two-stage Markov chain Monte Marlo simulation and polynomial chaos expansion,  Water Resour. 49 (2013) 2664–2682.
  * Brooks et al. [1998] S. P. Brooks, et al.,  MCMC convergence diagnosis via multivariate bounds on log-concave densities,  Annals of statistics 26 (1998) 398–433.
  * Polson [1996] N. G. Polson,  Convergence of Markov chain Monte Carlo algorithms,  Bayesian statistics 5 (1996) 297–321.
  * Rosenthal [1995] J. S. Rosenthal,  Minorization conditions and convergence rates for Markov chain Monte Carlo,  Journal of the American Statistical Association 90 (1995) 558–566.
  * Roy [2019] V. Roy,  Convergence diagnostics for Markov chain Monte Carlo,  Annual Review of Statistics and Its Application 7 (2019) 387–412.
  * M. K. Cowles [1996] B. P. C. M. K. Cowles,  Markov chain Monte Carlo convergence diagnostics: A comparative review,  J. of American Statistical Association 91 (1996) 883–904.
  * Mengersen et al. [1999] K. L. Mengersen, C. P. Robert, C. Guihenneuc-Jouyaux,  MCMC convergence diagnostics: A review,  in: M. Bernardo, J. O. Berger, A. P. Dawid, A. F. M. Smtith (Eds.), Bayesian Statistics, volume 6, Oxford University Press, 1999, pp. 415–440.
  * Brooks and Gelman [1998] S. Brooks, A. Gelman,  General methods for monitoring convergence of iterative simulations,  Journal of Computational and Graphical Statistics 7 (1998) 434–455.
  * Cotter et al. [2013] S. L. Cotter, G. O. Roberts, A. M. Stuart, D. White,  MCMC methods for functions: Modifying old algorithms to make them faster,  Statistical Science 28 (2013).
  * Durlofsky [1991] L. Durlofsky,  Numerical calculation of equivalent grid block permeability tensors for heterogeneous porous media,  Water Resour. Res. 27(5) (1991) 699–708.
  * Ali et al. [2021] A. Ali, A. Al-Mamun, F. Pereira, A. Rahunanthan,  Conditioning by projection for the sampling from prior Gaussian distributions,  in: International Conference on Computational Science and Its Applications, Springer, 2021, pp. 653–669.
  * Delhomme [1978] J. P. Delhomme,  Kriging in the hydrosciences,  Advances in water resources 1 (1978) 251–266.
  * Journel and Huijbregts [1978] A. G. Journel, Huijbregts, Mining Geostatistics, Technical Report, Academic Press, 1978.
  * Strang [2019] G. Strang, Linear algebra and learning from data, Wellesley-Cambridge Press Cambridge, 2019\.
  * Smith [2007] B. J. Smith,  boa: An R package for MCMC output convergence assessment and posterior inference,  Journal of Statistical Software 21 (2007) 1–37.
  * Ali et al. [2020] A. Ali, H. Mankad, F. Pereira, F. S. Sousa,  The multiscale perturbation method for second order elliptic equations,  Applied Mathematics and Computation 387 (2020) 125023.