Multilevel approximate Bayesian approaches for flows in highly heterogeneous porous media and their applications

Journal of Computational and Applied Mathematics (

)

–

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Multilevel approximate Bayesian approaches for flows in highly heterogeneous porous media and their applications Nilabja Guha a,b , Xiaosi Tan c,∗ a

Department of Statistics, Texas A & M University, College Station, TX 77843, USA

b

ISC, Texas A&M University, College Station, TX 77843, USA

c

Department of Mathematics, Texas A&M University, College Station, TX 77843, USA

article

info

Article history: Received 4 July 2015 Received in revised form 10 October 2016 Keywords: Multiscale MCMC Porous media Two-phase flow MLMC Approximate Bayesian

abstract Estimation of quantities related to high-contrast flow problems such as permeability field plays an important role in porous media characterization. A Generalized Multiscale Finite Element Method (GMsFEM) can be used for solving parameter-dependent (or stochastic) flow problems with multiscale nature. A hierarchy of approximations of different resolution can be provided by GMsFEM. Hence, it can be coupled with Multilevel Markov Chain Monte Carlo (MLMCMC) to generate samples in different levels and form the multilevel estimator. Karhunen–Lo’eve Expansion (KLE) is used to parameterize the underlying random field by a function of Gaussian random field. Instead of MCMC, an Approximate Bayesian Computation (ABC) method can be used within the Multilevel Monte Carlo framework. ABC can be incorporated in different levels to reduce the computational cost and to produce an approximate solution by ensembling different levels. © 2016 Elsevier B.V. All rights reserved.

1. Introduction Reservoir performance forecasting is greatly influenced by the uncertainties in the description of the underlying permeability field. The permeability field is modeled as a spatial random process and the uncertainty is reduced by incorporating observed data (pressure/fractional flow). A Bayesian framework gives a natural structure for updating the random field based on the observed data. The process is computationally intensive as it requires a large number of forward simulations. In the multiscale problems of high contrast field estimation, a Multilevel Markov Chain Monte Carlo (MLMCMC) method integrated with the Generalized Multiscale Finite Element Method (GMsFEM) can be used [1], in which the underlying field is estimated by samples from different levels. In this paper, we implement a multilevel ABC algorithm to estimate the underlying field. A new mixed GMsFEM is implemented for multiscale estimations, which provides a hierarchy of approximations to the solution that can be utilized in the multilevel framework. The basic idea of our Multilevel Monte Carlo (MLMC) method is first introduced in [2] for finite- and infinite-dimensional integration. Later on, it was applied to stochastic ODEs and PDEs in [3–6]. And in [1], the idea of MLMC is integrated with the Metropolis–Hastings algorithm and the GMsFEM to form a multilevel MCMC framework. In MLMC, a respective number of samples are used at different levels and the estimate can be formed in terms of a telescoping sum. More realizations are used at the coarser levels with inexpensive forward computations, and fewer samples are needed at the finer and more expensive levels due to the smaller variances.

∗

Corresponding author. E-mail addresses: [email protected] (N. Guha), [email protected], [email protected] (X. Tan).

http://dx.doi.org/10.1016/j.cam.2016.10.008 0377-0427/© 2016 Elsevier B.V. All rights reserved.

2

N. Guha, X. Tan / Journal of Computational and Applied Mathematics (

)

–

The underlying fine-scale problem considered in this paper is a coupled flow and transport equations. The flow equations require mass conservative discretization to avoid suspicious saturation values in the convection equation. For this reason, we use mixed finite element discretization and its multiscale version. A mixed GMsFEM is described for solving the parameterdependent two-phase flow equation with transport. In this method, the pressure field is approximated by piecewise constant basis functions, while the velocity field is approximated by multiscale basis which is constructed with an offline–online procedure similar to the GMsFEM introduced in [1]. A hierarchy of approximations to the solution is provided by the resulting mixed GMsFEM, which can be integrated into the multilevel Monte Carlo framework. Approximate Bayesian Computation [7] is a novel way to approximate Bayesian computations where the likelihood does not have a closed form and is expensive to compute. Using the known data generating mechanisms, given the parameter, a candidate data set is generated. The parameter value is feasible when the generated data set is close enough to the observed data. The closeness in that acceptance step controls the degree of approximation. Because of no burden of likelihood computation and flexible nature, it is adaptable in many different scenarios. Some of the recent developments include adaptive approximation [8] and application in linear regression model [9]. In this article we implement ABC in a multiscale scenario. The underlying permeability field is estimated in a multilevel set up. We use a novel mixed GMsFEM for the multiscale procedure. A multilevel ABC (MLABC) algorithm is proposed and the mean field is estimated from the posterior sample generated by the multilevel ABC. Some relevant results and justification related to the MLABC algorithm are also proved for some theoretical insight. Results about the stationary property of the chain at each level generated by the ABC algorithm are shown and the approximation property is discussed. We implement the proposed methodology in two-phase and single-phase flow problems. The proposed ABC based method shows improvement of computational cost. Also, from simulation it shows that a higher number of observations are accepted from the higher level. Thus, it suggests the scope of a multilevel ABC algorithm in the context of flow estimation. The article is arranged as follows. In Section 2, we briefly introduce the model and the probabilistic framework. Then, we describe mixed GMsFEM and MLMC in Section 3. In Section 4, we explicitly write down the permeability characterization in terms of KLE. Later, we discuss ABC thoroughly and write down the stationary distribution in our multilevel case in Section 5. In Section 6, we compare the approximation and computational gain in both single- and two-phase flow examples. 2. Model We consider two-phase flow in a subsurface formation (denoted by D) under the assumption that the displacement is dominated by viscous effects. The effects of gravity, compressibility and capillary pressure are neglected. The two phases are referred to as water and oil, designated by subscripts w and o, respectively. Darcy’s law can be written as follows for each phase, with all quantities dimensionless:

κrj (S ) κ · ∇ u, (1) µj where vj , j = w, o, is the phase velocity, κ is the permeability tensor, κrj is the relative permeability to phase j, S is the water vj = −

saturation (volume fraction) and u is the pressure. Combining Darcy’s law with a statement of conservation of mass allows us to express the governing equations in terms of pressure and saturation equations:

∇ · (λ(S )κ∇ u) = Qs , ∂S + v · ∇ fs (S ) = 0, (2) ∂t where λ(S ) is the total mobility, Qs is the source term, v is the total velocity and fs (S ) is the fraction flux of water, which are respectively given by:

κr w (S ) κro (S ) + , µw µo v = vw + vo = −λ(S )κ · ∇ u, κr w (S )/µw fs (S ) = . κr w (S )/µw + κro (S )/µo λ(S ) =

(3) (4) (5)

The above description is referred to as the fine model of the two-phase flow problem. For single-phase flow, we have λ(S ) = 1 and fs (S ) = S. Throughout, the porosity is assumed to be constant. In this article, the permeability field will be sampled based on fractional flow, specifically, the fraction of water produced in relation to the total production rate, denoted by F (t ) (or F in further discussion). F for a two-phase water–oil flow is defined as the fraction of water in the produced fluid and is given by qw /qt , where qt = qo + qw , with qo and qw the flow rates of oil and water at the production edge of the model, i.e., out vn fs (S )dl , F (t ) = ∂ D ∂ Dout vn dl



where ∂ Dout is outflow boundaries and vn is the normal velocity field.

(6)

N. Guha, X. Tan / Journal of Computational and Applied Mathematics (

)

–

3

For the high contrast single flow equation on D in Rd for d = 2, 3, we have a simpler representation,

−∇.(κ(x, µ)∇ u) = f in D,

(7)

with suitable boundary condition, where f is the source term and κ(x, µ) is the heterogeneous permeability field. Here, u is the pressure inside the medium. The underlying field is assumed to be high contrast and µ is the multi-dimensional parameter determining the field. In this case, we will sample the permeability field conditioned on the observed pressure data. Given the observed data Fobs (fraction flow or pressure), our aim is to sample from and estimate the permeability field κ . Because of the PDE structure and the error in observation, a direct solution is not feasible. A Bayesian framework is used for this purpose. The main idea is to model the field of interest by some appropriate spatial process, and exploit the PDE structure and the error distribution to find the distribution of the permeability field given the observed points. As mentioned earlier, we do not observe the true data of quantity of interest. The observed data points are convoluted with some random error. We assume the true data F (κ) is observed with some random error ϵ and the error terms are assumed to be independent normal. That is, Fobs = F (κ) + ϵ,

(8)

ϵ ∼ N (0, σf2 ). Here, κ is related to the flow equation through the PDE. It is parametrized by Karhunen–Loéve expansion (KLE) method with parameter θ , determining κ . Under prior p(.) on θ (or equivalently κ ), the posterior is

π (κ) = p(κ|Fobs ) ∝ p(Fobs |κ)p(θ ).

(9)

Because of the underlying PDE, a closed form posterior sampling is not possible. Simple Metropolis–Hastings method amounts to solving the PDE at each draw in fine scale which is computationally intensive. A multiscale method is therefore appropriate to save the computational cost, using KLE for parameterization of the random field. A similar parameterization can be adapted for the ABC approximation. Instead of direct sampling using MCMC, in a ABC method we generate data for the candidate κ and check if it is close to the true observed data in some appropriate metric. A more detailed account is given in latter sections. Using ABC samples, we implement MLABC and estimate the κ . 3. Mixed GMsFEM as a multilevel solver We introduce mixed GMsFEM which achieves efficient forward model simulation and also provides a hierarchy of approximation which can be used for constructing multilevel Monte Carlo estimates. More details regarding GMsFEM can be found in [10–12]. In [1], GMsFEM for single-phase flow is introduced as a multilevel solver, which is coupled with the MLMC methods to arrive at a general framework for the uncertainty quantification of the quantities of interest in multiscale flow problems. In this section, a mixed GMsFEM for the two-phase flow model problem (2) will be introduced. 3.1. Preliminaries We consider the following two-phase flow problem with transport in mixed formulation

κ(x, µ)−1 v + ∇ u = 0 in D, div (v) = f in D,

(10)

with non-homogeneous Neumann boundary condition v · n = g on ∂ D, where κ(x, µ) is the permeability field, D is the computational domain, and n is the outward unit-normal vector on ∂ D. In the mixed GMsFEM considered in this chapter, we construct the basis functions for the velocity field v with an offline–online procedure. For the pressure field u, we will use piecewise constant approximations. First, we introduce the following definitions of coarse and fine grids to discretize the model problem (10). Let T H be a conforming partition of the computational domain D into finite elements (triangles, quadrilaterals, tetrahedra, etc.). This partition is referred to as the coarse grid, and each coarse subregion is further partitioned into a connected union of fine grid N blocks. We denote the fine grid partition by T h . We use {Ei }i=e 1 (where Ne is the number of coarse edges) to denote the set H of all edges of the coarse mesh T , and define the neighborhood ωi corresponding to the coarse edge Ei by

ωi =

 {Kj ∈ T H ; Ei ∈ ∂ Kj }.

(11)

See Fig. 1 for an illustration of neighborhoods and elements subordinated to the coarse and fine discretization. Let QH be the space of piecewise constant functions with respect to the coarse grid T H , the approximation of the pressure u will be obtained in this space. We define the multiscale space for velocity field v as the linear span of all local multiscale

4

N. Guha, X. Tan / Journal of Computational and Applied Mathematics (

)

–

Fig. 1. Illustration of a coarse edge Ei and the associated coarse neighborhood ωi . 0 basis functions corresponding to edge Ei , which is denoted as VH = ξ H {Ψi }. We also define VH = VH ∩ {v ∈ VH : v · n = 0 on ∂ D}. The mixed GMsFEM is to find (vH , uH ) ∈ VH × QH such that





κ

−1

vH · wH −

D



div(wH )uH

=

div(vH )qH

=

0,

∀wH ∈ VH0 ,

D



(12)



D

fqH ,

∀qH ∈ QH ,

D

where vH · n = gH on ∂ D, and for each coarse edge Ei ∈ ∂ D, we have E (gH − g )Ψj · n = 0 for all basis functions Ψj i corresponding to the edge Ei . With the above definition, we proceed to describe the offline–online computational procedure which will construct the multiscale space VH for the velocity field v .



3.2. The construction of multiscale basis functions At the offline stage, we first generate the snapshot space, and then a lower-dimensional offline space. This offline space Ei is computed once and used repeatedly at the online stage. We first construct a snapshot space Vsnap on each coarse edge Ei , (i)

(i)

which involves solving a set of local problems for a number of selected parameters. Specifically, we will find (vl,j , ul,j ) by solving the following problem on each coarse neighborhood ωi corresponding to the edge Ei :

κ(µj )−1 vl(,ij) + ∇ ul(,ij) = 0 in ωi , (i)

(i)

(13)

in ωi ,

div (vl,j ) = αl,j

(i)

subject to the boundary condition vl,j · ni = 0 on ∂ωi , where ni denotes the outward unit-normal vector on ∂ωi . Here, {µj }Jj=1 is a set of selected parameters. The above problem (13) will be solved separately in the coarse-grid blocks forming ωi . Therefore, we will need an extra boundary condition on Ei . Notice that Ei can be written as a union of fine-grid edges, Li Ei = l=1 el , where Li is the total number of fine-grid edges on Ei and el denotes a find-grid edge. We define a piecewise (i)

constant function ∆l on Ei as (i)

∆l =

1, 0,



on el , on other fine grid edges on Ei ,

for l = 1, 2, . . . , Li . The remaining boundary condition on the coarse edge Ei for the local problem (13) is then taken as

vl(,ij) · mi = ∆(l i) ,

(14)

where mi is a fixed unit-normal vector on Ei . The constant

α (i) Kn l,j



=

 Ei

(i)

∆l for all Kn ⊂ ωi .

αl(,ij)

in (13) is chosen to satisfy the compatibility condition i,snap

The collection of the solutions of the above local problems generates the snapshot space. Let Ψl,j snapshot fields. The snapshot space Vsnap for each ωi is defined by i,snap

Vsnap = span{Ψl,j

:= vl(,ij) be the

: 1 ≤ j ≤ J , 1 ≤ l ≤ Li }.

snap Ψi

snap

snap

Notice that each can be represented by a vector ψi containing the coefficient in the expansion of Ψi in the finegrid basis functions. After the construction of the snapshot space, we proceed to construct the offline space by performing a dimension reduction within the space of snapshots using some local spectral problems following the general framework of [10]. Let (i) Vsnap be the snapshot space corresponding to the coarse edge Ei , which is defined by i,snap

(i) Vsnap = span{Ψj

: 1 ≤ j ≤ Msnap }.

N. Guha, X. Tan / Journal of Computational and Applied Mathematics (

)

–

5

(i)

We consider the local spectral problem: find a real number λ ≥ 0 and a function v ∈ Vsnap such that ai (v, w) = λsi (v, w),

(i) ∀w ∈ Vsnap .

(15)

Here, we take ai (v, w) =



κ −1 v · w,

ωi

si (v, w) =



  uv

uw ,



(16)

Ei

where κ is the domain-based parameter-averaged coefficient, (v, uv ), (w, uw ) are solutions of the local problem (13), and   p denotes the jump of the function p. To generate the offline space, we choose the smallest Moff eigenvalues to (15) and i,off

take the corresponding eigenvectors, ψk

i,off

, in the space of snapshots by setting Ψk i,off

to form the reduced snapshot space, where ψkj i,off

Voff = span{Ψk

=

i,off

are the coordinates of the vector ψk



j

ψkji,off Ψji,snap for k = 1, . . . , Moff

. Then the offline space is

: 1 ≤ k ≤ Moff , 1 ≤ i ≤ Ne },

where Ne is the total number of edges in the domain. One can solve the model problem using the offline space at the online stage. But here we will further construct an online space Von for each given input parameter µ. We want it to be a low-dimensional subspace of the offline space for computational efficiency. In particular, we seek a subspace of the offline space by solving a same eigenvalue problem to (15), but solved in the offline space Voff , and in (16) we use the input κ(µ) instead of κ . Similarly, eigenvectors corresponding to the smallest Mon eigenvalues will be chosen. We note that, at the online stage, the bilinear forms as in (16) will be formed with the input parameter µ, which means that they are parameter-dependent. Mon can be chosen such that different size of the online multiscale space will be formed, which provides the hierarchy approximations. This online space Von is used to solve the model problem. Remark 1. A postprocessing technique is applied in ourexperiments to obtain conservative velocity fields on the fine-grid.  We notice that the solution of (12) satisfies ∂ K vH · n = K f for every coarse-grid block K . When f has fine-scale oscillation in some coarse blocks, the velocity field needs to be postprocessed in these coarse blocks. Apostprocessed velocity vh⋆ will  be constructed such that conservation on the fine grid is obtained, as follows: ∂τ (vh⋆ · n) = ∂τ f , ∀τ ∈ T h . Above, the offline–online procedure of the mixed GMsFEM is introduced for the model problem (10). To solve the twophase flow and transport problem, we solve the flow equation for each time step and use the fine-scale velocity to advance the saturation front. More details and discussions about mixed GMsFEM can be found in [13]. 4. Permeability parameterization To obtain a permeability field in terms of an optimal L2 basis, we use the KLE [14]. For our numerical tests, we truncate the expansion and represent the permeability matrix by a finite number of random parameters. We consider the random field Y (x, ω) ˜ = log[k(x, ω)] ˜ , where ω˜ represents randomness. We assume a zero mean E[Y (x, ω)] ˜ = 0, with a known covariance operator R(x, y) = E [Y (x)Y (y)]. Then, we expand the random field Y (x, ω) ˜ as Y (x, ω) ˜ =

∞ 

Yk (ω) ˜ Φk (x),

k=1

with Yk (ω) ˜ =

 Ω

Y (x, ω) ˜ Φk (x)dx.

The functions {Φk (x)} are eigenvectors of the covariance operator R(x, y), and form a complete orthonormal basis in L2 (Ω ), i.e.,

 Ω

R(x, y)Φk (y)dy = λk Φk (x),

k = 1, 2, . . . ,

(17)

√

where λk = E[Yk2 ] > 0. We note that E[Yi Yj ] = 0 for all i ̸= j. By denoting γk = Yk / λk (whence E[γk ] = 0 and E[γi γj ] = δij ), we have Y (x, ω) ˜ =

∞   k=1

λk γk (ω) ˜ Φk (x),

(18)

6

N. Guha, X. Tan / Journal of Computational and Applied Mathematics (

)

–

where Φk and λk satisfy (17). The randomness is represented by the scalar random variables γk . After discretizing the domain Ω into a rectangular mesh, we truncate the KLE (18) to a finite number of terms. In other words, we keep only the leadingorder terms (quantified by the magnitude of λk ), and capture most of the energy of the stochastic process Y (x, ω) ˜ . For an N-term KLE approximation YN =

N  

λk γk Φk ,

k =1

the energy ratio of the approximation is defined by N 

e(N ) :=

N 2

E ∥Y ∥

E ∥Y ∥2

=

λk

k=1

∞ 

. λk

k=1

If the eigenvalues {λk } decay very fast, then the truncated KLE with the first few terms would be a good approximation of the stochastic process Y (x, ω) in the L2 sense. Prior on the random field Gaussian random field is a natural choice for prior on a spatial field. Here, we assume that the log of the permeability field follows a Gaussian random field with a known covariance kernel. Here, the KLE expansion of the process is given in (18). Thus, we have the following prior: Y (x, ω) ˜ =

∞   λk γk (ω) ˜ Φk (x), k=1

γk ∼ N (0, 1),

(19)

where γk ’s are independent. 5. MCMC and ABC approximation In MCMC, the stationary distribution is achieved by a Metropolis–Hastings or a similar algorithm, where the balanced equation corresponding to the stationary distribution is satisfied. Consider a simple MCMC, starting with proposal density q(θ1 |θ2 ) which moves θ2 to θ1 . For observed data Y, where the likelihood and prior is given by p(Y|θ ) and p(θ ), respectively and target posterior density is given by π (θ |Y) ∝ p(θ , Y), MCMC steps can be written down as follows. We start with initial θ0 .

• By the transition kernel q(.) we move θk → θ ∗ . • The move from θk → θ ∗ is accepted with probability   q(θk |θ ∗ )p(θ ∗ , Y) . α(θk , θk+1 ) = min 1, q(θ ∗ |θk )p(θk , Y) • If accepted, θk+1 = θ ∗ , otherwise θk+1 = θk . Though this sampling is very useful for many practical purposes for obtaining the posterior distribution as the stationary distribution, it may encounter difficulty in case of inverse problem where we need to solve an underlying PDE at each draw. To reduce the computational cost, an alternative route can be employed through an approximate Bayesian approach. The main idea behind ABC is to generate a pseudo observation Z at each step given the proposed θ ∗ . If the distance between Y and Z, d(Y, Z) is less than some chosen η > 0, then we accept θ ∗ in the first step and move to the second step. The acceptance criterion for the second step is given in the algorithm below. If Y has a discrete distribution, then η = 0 gives the exact posterior density. The choice of tolerance limit η is subjective and it determines the order of the approximation. Hence, the procedure can be given as follows. We start with initial θ0 .

• By the transition kernel q(.) we move θk → θ ∗ . • Generate Z from the p(.|θ ) with θ = θ ∗ . • If d(Y, Z) < η, the move from θk → θ ∗ is accepted with probability   q(θk |θ ∗ )p(θ ∗ ) . α(θk , θk+1 ) = min 1, q(θ ∗ |θk )p(θk ) • If accepted, θk+1 = θ ∗ , otherwise θk+1 = θk .

N. Guha, X. Tan / Journal of Computational and Applied Mathematics (

)

–

7

The resulting stationary distribution is not exactly p(θ|Y). Instead, it is given by

πη (θ , Z|Y) = 

p(θ , Z)Id(Y,Z)<η AηY ×θ

p(θ , Z)dZdθ

,

where AYη = {Z : d(Y, Z) < η}. Note that if η approaches zero, the stationary distribution approaches the posterior distribution. Therefore, we obtain an approximate posterior and the choice of η determines the trade-off between computational cost and accuracy. 5.1. ABC for multilevel MCMC In the multilevel scenario, we use a similar approach to ordinary ABC. Let η1 , . . . , ηl be the cut-off of the affinity between true Y, Fobs , and the pseudo data points FlZ (κ) = Z’s generated by solving the coarse level problem at lth level. The ηl ’s from different levels may be different and produce an approximation of the true posterior at each level. The algorithm is formally given in Algorithm 1. Here, κ parametrized by γk ’s is our parameter and we will replace θ by κ for notational convenience. The prior for level l is given by the prior on γk ’s and is denoted by pl (κ) = p(κ). Let ξ l be the set such that the target density at level l has density value greater than zero. Also, let Dl be the set where the proposal density is greater than zero. Under the normality of (8), and the proposal density given by random perturbation around the current value q0 (κ|κn ) = κn + δϵn , the condition ξ = ξ l = Dm = D is satisfied. Here D is the set where the proposal density is greater than zero for all m and D equals the whole domain, as normal distribution has density greater than zero at all values. Similarly for the target density ξ l = ξ = whole domain. Also, the resulting chain is aperiodic at each level l [1]. Finally, we can state the following result about the resulting stationary distribution. Lemma 1. Under the model (8) and prior (18), the stationary distribution corresponding to the ABC algorithm in Algorithm 1 is given by p (κ, Z)Id(Fobs ,Z )<η

l η πl (κ, Z|Fobs ) = 

F

obs Aη,κ

pl (κ, Z)dZdθ

, F

obs where Z is a draw from the field value κ at level l, and Aη,κ denotes the set of observations y = Fl (κ) + ϵ which have a distance less than η = ηl from the observed Fobs .

Proof. We show the result for l + 1th level. Let c be the value of the integral in the denominator at l + 1th level. Without m m ′ loss of generality, we assume ρl+1 (κlm +1 , κ) ≤ 1. For a transition (κl+1 , Z) → (κ, Z ), the transition probability for κl+1 ̸= κ is m m Ql+1 (κlm +1 , κ) = ql (κ|κl+1 )ρl+1 (κl+1 , κ).

Let I denote the indicator function that Z and Z′ is in A for level l from Algorithm 1, we have for level l + 1

Fobs

η,κlm+1

F

obs , respectively. Writing down the stationary condition and Aη,κ

η

πl+1 (κlm+1 , Z|Fobs )Ql+1 (κlm+1 , κ)pl+1 (Z′ |κ)I = c −1 pl+1 (Z|κlm+1 )pl+1 (κlm+1 )ql (κ|κlm+1 )ρl+1 (κlm+1 , κ)pl+1 (Z′ |κ)I ql (κlm +1 |κ)pl+1 (κ) pl+1 (Z′ |κ)I = c −1 pl+1 (Z|κlm+1 )pl+1 (κlm+1 )ql (κ|κlm+1 ) m ql (κ|κlm +1 )pl+1 (κl+1 ) = c −1 pl+1 (Z′ |κ)pl+1 (κ)ql (κlm+1 |κ)pl+1 (Z|κlm+1 )I = c −1 pl+1 (Z′ |κ)pl+1 (κ)ρl+1 (κ, κlm+1 )pl+1 (Z|κ)I η

= πl+1 (κ, Z′ |Fobs )Ql+1 (κ, κlm+1 )pl+1 (Z|κlm+1 )I. For κlm +1 = κ , the condition is satisfied trivially.



Remark 2. The marginal posterior of k would be therefore p(κ)tA (ηl , κ)

η πl (κ|Fobs ) = 

F

Aηobs

where tA (ηl , θ ) =



F

Aηobs ,κ l

pl (κ, Z)dZdκ

pl (Z|κ)dZ.

,

(20)

8

N. Guha, X. Tan / Journal of Computational and Applied Mathematics (

)

–

The choice of ηl ’s determines the order of approximation of the posterior distribution at the lth level. With ηl approaching η zero, the true posterior is achieved. If some regularity conditions are satisfied, the ABC estimate πl m (κ|Fobs ) approaches the true posterior πl (κ|Fobs ) at the lth level in some appropriate topology, with ηl converging to zero. A simple result for weak topology can be stated in the form of the following lemma. η

Lemma 2. Suppose p(Z|κ) is a continuous function in Z and supZ,κ p(Z|κ) < ∞. Then, the ABC estimate πl m (κ|Fobs ) converges to the true posterior πl (κ|Fobs ) at the lth level in terms of weak topology, as ηl converges to zero. Proof. For notational convenience, we drop the subscript l. Let p(Z|k) < co < ∞ for all κ and Z. Let Vn (η) be the volume of F the n dimensional set Aηobs = {Z : d(Fobs , Z) < η}. Here, n be the length of the observed vector. Let M be the prior dimension used in KL expansion. Given δ > 0, there exists a compact set Kc such that p(Kc′ ) < δ . F Given δ ′ > 0, there exists η′ > 0 such that for all η < η′ and Z ∈ Aηobs , we have |p(Fobs |κ)− p(Z|κ)| < δ ′ for κ ∈ Kc . Let Nη and Dη be the numerator and denominator of (20). Hence, we have |Dη − m(Fobs )Vn (η)| < 2δ ′ Vn (η) + 2Vn (η)co δ . Here, m() denotes the marginal density. Also, |Nη (κ)−p(Fobs , κ)Vn (η)| ≤ δ ′ Vn (η) on Kc and on Kc′ , |Nη (κ)−p(Fobs , κ)Vn (η)| < 2co Vn (η). η Thus, on Kc , πl (κ|Fobs ) = πlκ (κ|Fobs ) + O(δ) + O(δ ′ ). Now for any bounded continuous function of θ , g (θ ), we have Eπ η (κ|Fobs ) (g (κ)) = Eπ η (κ|Fobs ) (g (κ)IKc ) + c δ ′ , l

l

and hence, Eπ η (κ|Fobs ) (g (κ)) = Eπl (κ|Fobs ) (g (κ)) + O(δ) + O(δ ′ ). l

As δ and δ ′ can be arbitrarily small, it concludes the proof.



Remark 3. For practical purposes, the choice of ηl can be data driven. One approach may be using some small quantile of the posterior distribution of the proposed metric. This can be done by solving the system in the coarse grid, which is computationally cheaper. 5.2. ML-ABC implementation The idea of multilevel MCMC is useful in context of drawing posterior samples from different levels. Functions of the parameter such as posterior mean may be of interest. In the MLMCMC, higher levels associate with higher cost and accuracy, while lower levels may be cheaper but less accurate. Thus, using different levels we can estimate the subsequent improvements of adding a new level, and estimate the target function as a telescoping sum. We start with the telescopic sum

 EπL [FL ] =

 =

FL (x)πL (x)dx F0 (x)π0 (x)dx +

L  

(Fl (x)πl (x) − Fl−1 (x)πl−1 (x))dx,

l=1

where πl denotes the approximated target distribution at level l, and π0 is our initial level. We note that after the initial level each expectation involves two measures, πl and πl−1 . We can rewrite the integration using a product measure and get,

Eπ L [FL ] = Eπ0 [F0 ] +

L 

Eπl ,πl−1 [Fl − Fl−1 ].

(21)

l=1

The idea of our multilevel method is to estimate each term of the right hand side of Eq. (21) independently. The proposed ABC estimator can be used to estimate each term in (21). To create different levels, we adapt the proposal distribution q(k|km ) to the target distribution π (κ) using the GMsFEM with different sizes of the online space. The process modifies the proposal distribution q(κ|κm ) by incorporating the online coarse-scale information. Let Fl (κ) be the pressure/production computed by solving the online coarse problem at level l for a given κ . The target distribution π (κ) is approximated on level l by πl (κ), with π (κ) ≡ πL (κ). Here, we have

  ∥Fobs − Fl (κ)∥2 × p(κ). πl (κ) ∝ exp − 2σl2

(22)

In the algorithm, we still keep the same offline space for each level. From level 0 to level L, we increase the size of the online space as we go to a higher level, which means for any levels l, l + 1 ≤ L, samples of level l are cheaper to generate than that of level l + 1. This idea underlies the cost reduction using the multilevel estimator. Hence, the posterior distribution for coarser levels πl , l = 0, . . . , L − 1, does not have to model the measured data as faithfully as πL , which in particular implies that by choosing a suitable value of σl2 it is easier to match the result Fl (κ) with the observed data. We denote the number of samples at level l by Ml , where we will have M0 ≥ · · · ≥ ML . As was discussed above, our quantity of interest can be approximated by the telescopic sum (21).

N. Guha, X. Tan / Journal of Computational and Applied Mathematics (

)

–

9

Algorithm 1 Multilevel ABC 1: 2: 3:

Pick η1 , η2 , . . . , ηL Given κ m , draw a trial proposal κ from distribution q(κ|κ1m ) = q0 (κ|κ1m ) Compute the acceptance probability

  q0 (κ1m |κ)p1 (κ) ρ1 (κ1m , κ) = min 1, q0 (κ|κ1m )p1 (κ1m ) 4: 5: 6: 7: 8: 9: 10:

Compute d(Fobs , F1Z (κ)) if d(Fobs , F1Z (κ)) ≤ η1 then Accept κ with probability ρ1 end if for l = 1 : L − 1 do if κ is accepted at level l then Form the proposal distribution ql (on the l + 1th level) by m m m ql (κ|κlm +1 ) = ρl (κl+1 , κ)ql−1 (κ|κl+1 ) + δκ (1 − l+1

11:

13: 14: 15: 16: 17: 18: 19: 20: 21: 22:

ρl (κlm+1 , κ)ql−1 (κ|κlm+1 )dκlm+1 )

Compute the acceptance probability

ρl+1 (κ 12:



m l +1



, κ) = min 1,

ql (κlm +1 |κ)pl+1 (κ)



m ql (κ|κlm +1 )pl+1 (κl+1 )

Compute d(Fobs , FlZ (κ)) if d(Fobs , FlZ (κ)) ≤ ηl then Accept κ with probability ρl+1 end if if κ is accepted then κlm++1 1 = κ and go to next level (if l = L − 1, accept κ and set κLm+1 = κ ) else κlm++1 1 = κlm+1 , and break end if end if end for

6. Numerical results In this section, we present some numerical examples of the multilevel ABC method. As mentioned earlier, the permeability field is given by κ(x, µ), where x is defined on the unit square Ω = [0, 1]2 . We assume that the permeability field κ(x, µ) is a log-normal process and its covariance function is known. We use the Karhunen–Loéve expansion (KLE) to parameterize the permeability field. Then, we apply the multilevel ABC method described earlier and compare the simulation results with multilevel MCMC method. Both single-phase and two-phase flow cases will be considered. In our examples, the permeability field κ is assumed to follow a log-normal distribution with a known spatial covariance, with the correlation function R(x, y) given by

 |x − y |2 1 1

R(x, y) = σ 2 exp −

2l21

−

|x2 − y2 |2  2l22

,

(23)

where l1 and l2 are the correlation lengths in x1 - and x2 -direction, respectively, and σ 2 = E[Y 2 ] is a constant that determines the variation of the permeability field. Choice of correlation length is important. Small correlation length can be helpful in modeling non-smooth field. For more about the choice of l1 and l2 and prior specifications, see [15]. In the following simulations, we use known l1 and l2 . Later, we consider an example in which we put prior on σ 2 and use it as an additional parameter in the ABC setup. We use the telescoping sum as mentioned in Section 5.2 to compute the mean permeability in MLABC algorithm. 6.1. Single-phase flow For the first set of the numerical tests, we consider the stationary, single-phase flow model (7) with f ≡ 1 and linear boundary conditions. We apply GMsFEM to solve the forward problem, with a 5 × 5 coarse grid and a 50 × 50 fine grid. The ‘‘true’’ data Fobs is obtained by generating a reference permeability field, solving the forward problem with the GMsFEM, and evaluating the pressure at nine selected points away from the boundary.

10

N. Guha, X. Tan / Journal of Computational and Applied Mathematics (

)

–

Table 1 Number of accepted samples at each level in MCMC, 3-level MCMC and 3-level ABC for the single-phase flow test case.

Total Level 1 Level 2 Level 3

MCMC

3-level MCMC

3-level ABC

5534

4723 2299 1573 1000

4197 1754 1176 1000

1000

1

1 3

3

2.5

0.8

2

0.6

1.5

2

0.4

1 0.5

0

0

0.2 0.4 0.6 0.8

1

0

MCMC

1

1

0.2 0

0

0.2 0.4 0.6 0.8

1

3-level MCMC

1 3 2.5

3

0.8

2

0.6

1.5

0.4

2

1 0.5

0

1

0

0.2 0.4 0.6 0.8

1

0

1

0.2 0

0

0.2 0.4 0.6 0.8

1

0

3-level ABC 3.5 3 2.5 2 1.5 1 0.5

0

0

Fig. 2. Numerical results for the single-phase flow test case. Top left: The true log-permeability field. Top right: Initial log-permeability field. Middle Left/Right: The mean of the sampled log-permeability field from MCMC and 3-level MCMC. Bottom: The mean of the sampled log-permeability field from 3-level ABC.

The prior permeability distribution p(κ) is parametrized by KL expansion as introduced in Section 4. We keep N = 5 terms in the KL expansion. The Gaussian field is of correlation length l1 = l2 = 0.1. Our proposal distribution is a random walk sampler in which the proposal distribution depends on the previous value of the permeability field and is given by q(κ|κn ) = κn + δϵn , where ϵn is a random perturbation with zero mean and unit variance, and δ is a step size. The random perturbations are imposed on the γk coefficients in the KL expansion. The distance d(Fobs , FlZ (κ)) in Algorithm 1 is taken as the L2 -norm. The acceptance criterion ηl is decided by 25% quantile of the L2 -norm of pseudo data value Z’s from the draws of a coarse level MCMC sampling. To compare the performance of different methods, we test the same problem with MCMC, 3-level MCMC and 3-level ABC, respectively. For the multilevel methods, the numbers of eigenvalues to generate the online space at each level are taken to be 4, 8 and 16. For the 3-level ABC, we set the acceptance criterion η = 0.06 for all the levels. We run the algorithms until 1000 total samples pass the final level of acceptance. The number of accepted samples is listed in Table 1. As in ABC, we use distance metric as the acceptance criterion, unlikely values may have higher distance in even coarser levels and get rejected in the coarser level. In MCMC the procedure depends on the likelihood instead of the pseudo data generation. Therefore, rejecting more ‘bad’ samples in coarser level, we can have a higher acceptance rate in higher level in an ABC approach. In Fig. 2, we plot the reference permeability field, initial permeability field, and the result posterior mean of the 3 methods. We notice that all 3 sample means are very close to the reference field. Hence, the proposed multi-level ABC method is of higher efficiency.

N. Guha, X. Tan / Journal of Computational and Applied Mathematics (

)

–

11

Table 2 Number of accepted samples at each level in MCMC, 2-level MCMC and 2-level ABC for the isotropic case in Example 1 of the two-phase flow test case.

Total Level 1 Level 2

MCMC

2-level MCMC

2-level ABC

5436

16 524 5 345 1 000

4240 2185 1000

1000

Table 3 Number of accepted samples at each level in MCMC, 2-level MCMC and 2-level ABC for the anisotropic case in Example 1 of the two-phase flow test case.

Total Level 1 Level 2

MCMC

2-level MCMC

2-level ABC

4142

6623 3002 1000

2995 1477 1000

1000

6.2. Two-phase flow and transport 6.2.1. Example 1 In these simulations, we present the performance of the multilevel ABC method with two-phase flow and transport problems. We consider the model problem (2) on Ω = (0, 1)2 with zero Neumann boundary condition. The flow equation is solved by mixed GMsFEM [13] with a 50 × 50 fine grid and a 5 × 5 coarse grid, and the saturation equation is solved on the fine grid by the finite volume method. The prior permeability distribution p(k) is also parametrized by KLE as above. The quantity of interest is taken to be the water-cut function F . One injecter at (0, 0) and one producer at (1, 1) are considered when we run the forward model in the reference permeability field to get the fractional flow. Two sets of Gaussian covariance functions are considered in this simulation.

• Isotropic Gaussian field with correlation lengths l1 = l2 = 0.1, and a KL dimension 8. • Anisotropic Gaussian field with correlation lengths l1 = 0.05, l2 = 0.1, and a KL dimension 8. Our proposal distribution is a random walk sampler in which the proposal distribution depends on the previous value of the permeability field and is given by q(κ|κn ) = κn + δϵn , where ϵn is a random perturbation with mean zero and unit variance, and δ is a step size. The random perturbations are imposed on the coefficients in the KL expansion. The distance d(Fobs , FlZ (κ)) in Algorithm 1 is taken as the L2 -norm. The acceptance criterion ηl is decided by 20% quantile of L2 -norm of pseudo data value Z’s from draws of a MCMC sampling. Again, to compare the performance of different methods we run MCMC, 2-level MCMC and 2-level ABC on the same problem. For the multilevel methods, the numbers of eigenvalues to generate the online space at each level are taken to be 2 and 8. For the 2-level ABC test, we set different acceptance criterion ηl for each level. In the isotropic case, η1 = 0.38 and η2 = 0.1. In the anisotropic case η1 = 0.4 and η2 = 0.15. We run the algorithms until 1000 samples are accepted on the last level. For the isotropic case, the number of accepted samples is listed in Table 2. We can see that the 2-level ABC method has a higher acceptance rate than the other two methods on the more expensive level. In Fig. 3, we plot the reference fraction flow, the initial sample and the sample mean of each method. In this example and all the following examples, the water cut is plotted versus pore volume injected (PVI). We observe that the mean estimate of the fractional flow is very close to the observed data. In Fig. 4, we plot the reference log-permeability field, the initial sample, and the posterior mean of the 3 methods. For the anisotropic case, the number of accepted samples is listed in Table 3. We plot the reference fraction flow, the initial sample and the sample mean of each method in Fig. 5. The reference log-permeability field, the initial sample, and the posterior means from the 3 methods are shown in Fig. 6. Similarly, we observe that the 2-level ABC method has a higher acceptance rate than the other two methods on the more expensive level and the mean estimate of the fractional flow is very close to the observed data. We also notice that the run time of the 2-level ABC of both the tests is around 15% faster than the 2-level MCMC simulations. By these results, we can see that we achieve higher computational efficiency with the proposed multi-level ABC method. 6.2.2. Example 2 In this example, we use different parameterization for generating reference field and prior specification in earlier setup. We use two different classes of bases for this purpose. We use the first 4 bases of the anisotropic KLE with correlation lengths l1 = 0.07, l2 = 0.1 to generate the reference field and the isotropic KLE with correlation lengths l1 = l2 = 0.1 is used to estimate the field with 8 terms in the KL expansion. Therefore, the reference and the prior have different basis. For the multilevel methods, the numbers of eigenvalues to generate the online space at each level are taken to be 2 and 8. For the 2-level ABC test, the acceptance criterion for level 1 and level 2 is η1 = 0.28 and η2 = 0.08, respectively. The eigenvalue plots for isotropic and anisotropic cases are given in Fig. 7.

12

N. Guha, X. Tan / Journal of Computational and Applied Mathematics (

)

–

Fig. 3. Water cut results for the isotropic case in Example 1 of the two-phase flow test case. Red line designates the reference water cut, the dash line designates the initial water cut and the green, the blue and the cyan lines designate water cut corresponding to mean of the sampled water cut from MCMC, 2-level MCMC and 2-level ABC, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4. Numerical results for the isotropic case in Example 1 of the two-phase flow test case. Top left: The true log-permeability field. Top right: Initial log-permeability field. Middle Left/Right: The mean of the sampled log-permeability field from MCMC and 2-level MCMC. Bottom: The mean of the sampled log-permeability field from 2-level ABC.

The relative performances of MCMC, 2-level MCMC and 2-level ABC are given. Table 4 shows the number of accepted samples of each method. The reference and fitted log permeability fields are given in Fig. 9, while the true and fitted water cuts are given in Fig. 8. It can be seen that MLABC and MLMCMC can capture the structure of κ in this case. Also, like earlier cases, the ABC based method offers higher efficiency.

N. Guha, X. Tan / Journal of Computational and Applied Mathematics (

)

–

13

Fig. 5. Water cut results for the anisotropic case in Example 1 of the two-phase flow test case. Red line designates the reference water cut, the dash line designates the initial water cut and the green, the blue and the cyan lines designate water cut corresponding to mean of the sampled water cut from MCMC, 2-level MCMC and 2-level ABC, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6. Numerical results for the anisotropic case in Example 1 of the two-phase flow test case. Top left: The true log-permeability field. Top right: Initial log-permeability field. Middle Left/Right: The mean of the sampled log-permeability field from MCMC and 2-level MCMC. Bottom: The mean of the sampled log-permeability field from 2-level ABC.

6.2.3. Example 3 A reference field which is not generated by KLE will be considered in this example, and it will be estimated by MCMC, 2-level MCMC and 2-level ABC with KLE for parameterization. Specifically, we consider the exponential of the reference field to be a smooth field generated with the equation log(κ(x, y)) = x(1 − x)y(1 − y). The prior permeability will be similarly

14

N. Guha, X. Tan / Journal of Computational and Applied Mathematics (

)

–

Fig. 7. Left column shows the eigenvalues in isotropic and the right column in the anisotropic cases, Example 2, two-phase flow. Eigenvalues decrease very fast and less than 0.1 after roughly 5 and 8 terms respectively. Table 4 Number of accepted samples at each level in MCMC, 2-level MCMC and 2-level ABC for Example 2 of the two-phase flow test case.

Total Level 1 Level 2

MCMC

2-level MCMC

2-level ABC

9968

10 720 5 906 1 000

6366 3713 1000

1000

Fig. 8. Water cut results for Example 2 of the two-phase flow test case. Red line designates the reference water cut, the dash line designates the initial water cut and the green, the blue and the cyan lines designate water cut corresponding to mean of the sampled water cut from MCMC, 2-level MCMC and 2-level ABC, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

parametrized by KLE as the previous examples, and here we take the isotropic field with correlation length l1 = l2 = 0.1, and a KL dimension 5. In the multilevel methods, the numbers of eigenvalues to generate the online space at each level are also 2 and 8. The acceptance criterion is set as η1 = 0.36 and η2 = 0.07 for level 1 and level 2 in the 2-level ABC, respectively. The relative performances of MCMC, 2-level MCMC and 2-level ABC are presented. Table 5 shows the number of accepted samples of each method. The reference and fitted log permeability fields are illustrated in Fig. 11, while the true and fitted water cuts are presented in Fig. 10. From the figures we observe that the feature of the reference smooth field can also be captured by MLMCMC and MLABC. This shows the capability of these methods to recover a reference field not generated with the same parameterization. We may raise the number of KL dimension to achieve a better match. Also, the ABC based method offers higher efficiency. 6.2.4. Example 4 In the following example, we consider our model with a different prior on σ 2 in (23). Similar to the isotropic case in Example 1, we set l1 = l2 = 0.1 in (23), and KL dimension is taken to be 8. But instead of taking a constant σ 2 , the prior of σ 2 is assumed to be an inverse Gamma distribution as σ 2 ∼ InvGamma(3, 1/2). For the multilevel methods, the numbers of

N. Guha, X. Tan / Journal of Computational and Applied Mathematics (

)

–

15

Fig. 9. Numerical results for Example 2 of the two-phase flow test case. Top left: The true log-permeability field. Top right: Initial log-permeability field. Middle Left/Right: The mean of the sampled log-permeability field from MCMC and 2-level MCMC. Bottom: The mean of the sampled log-permeability field from 2-level ABC.

Table 5 Number of accepted samples at each level in MCMC, 2-level MCMC and 2-level ABC for Example 3 of the two-phase flow test case.

Total Level 1 Level 2

MCMC

2-level MCMC

2-level ABC

14 758

16 861 8 546 2 000

7735 3456 2000

2 000

Table 6 Number of accepted samples at each level in MCMC, 2-level MCMC and 2-level ABC for Example 4 of the two-phase flow test case.

Total Level 1 Level 2

MCMC

2-level MCMC

2-level ABC

12 024

13 598 8 113 2 000

11 569 4 860 2 000

2 000

eigenvalues to generate the online space at each level are 2 and 8. The acceptance criterion is set as η1 = 0.31 and η2 = 0.12 for level 1 and level 2 in the 2-level ABC, respectively. Table 6 shows the number of accepted samples of each method. The reference and fitted log permeability fields are illustrated in Fig. 13, while the true and fitted water cuts are presented in Fig. 12. From this case we can observe that MLABC and MLMCMC are capable to capture the structure of κ with an inverse gamma prior on σ 2 , with MLABC offering higher efficiency.

16

N. Guha, X. Tan / Journal of Computational and Applied Mathematics (

)

–

Fig. 10. Water cut results for Example 3 of the two-phase flow test case. Red line designates the reference water cut, the dash line designates the initial water cut and the green, the blue and the cyan lines designate water cut corresponding to mean of the sampled water cut from MCMC, 2-level MCMC and 2-level ABC, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 11. Numerical results for Example 3 of the two-phase flow test case. Top left: The true log-permeability field. Top right: Initial log-permeability field. Middle Left/Right: The mean of the sampled log-permeability field from MCMC and 2-level MCMC. Bottom: The mean of the sampled log-permeability field from 2-level ABC.

7. Discussion The proposed ABC based MLMCMC has been used successfully in single-phase and two-phase flow problems. The goal of using ABC is to gain computational efficiency without sacrificing the accuracy of the estimation. The estimation of the mean field κ is comparable to direct 2-level MCMC as seen in Figs. 2, 3 and 5. The computational efficiency is evident from Tables 1–3, where the number of acceptance in the higher level has increased significantly for the higher level which is

N. Guha, X. Tan / Journal of Computational and Applied Mathematics (

)

–

17

Fig. 12. Water cut results for Example 4 of the two-phase flow test case. Red line designates the reference water cut, the dash line designates the initial water cut and the green, the blue and the cyan lines designate water cut corresponding to mean of the sampled water cut from MCMC, 2-level MCMC and 2-level ABC, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 13. Numerical results for Example 4 of the two-phase flow test case. Top left: The true log-permeability field. Top right: Initial log-permeability field. Middle Left/Right: The mean of the sampled log-permeability field from MCMC and 2-level MCMC. Bottom: The mean of the sampled log-permeability field from 2-level ABC.

more expensive. This gives a scope for more accurate posterior analysis with more posterior samples from the finer grid. From limited simulation, the proposed ABC is faster than the MCMC-based method. The results shown in Section 5.1 and Remark 2 give a theoretical justification of the approximation. A more extensive study can be conducted about the choice of approximation parameter η and the theoretical properties can be explored further. The ABC based method has been implemented with a mixed GMsFEM based multiscale procedure. From this exercise, ABC shows promise for multiscale

18

N. Guha, X. Tan / Journal of Computational and Applied Mathematics (

)

–

procedures as an alternative to MCMC and applying ABC for general flow characterization can be an interesting area to study. References [1] Yalchin Efendiev, Bangti Jin, Presho Michael, Xiaosi Tan, Multilevel Markov Chain Monte Carlo method for high-contrast single-phase flow problems, Commun. Comput. Phys. 17 (01) (2015) 259–286. [2] Stefan Heinrich, Multilevel Monte Carlo methods, in: Svetozar Margenov, Jerzy Waśniewski, Plamen Yalamov (Eds.), Lecture Notes in Computer Science, Vol. 2179, Springer-Verlag, Berlin, 2001, pp. 58–67. [3] Michael B. Giles, Multilevel Monte Carlo path simulation, Oper. Res. 56 (3) (2008) 607–617. [4] Mike Giles, Improved multilevel Monte Carlo convergence using the Milstein scheme, in: Monte Carlo and Quasi-Monte Carlo Methods 2006, Springer, Berlin, 2008, pp. 343–358. [5] Andrea Barth, Christoph Schwab, Nathaniel Zollinger, Multi-level Monte Carlo finite element method for elliptic PDEs with stochastic coefficients, Numer. Math. 119 (1) (2011) 123–161. [6] K.A. Cliffe, M.B. Giles, R. Scheichl, A.L. Teckentrup, Multilevel Monte Carlo methods and applications to elliptic PDEs with random coefficients, Comput. Vis. Sci. 14 (1) (2011) 3–15. [7] Katalin Csilléry, Michael G.B. Blum, Oscar E. Gaggiotti, Olivier François, Approximate Bayesian computation (abc) in practice, Trends Ecol. Evol. 25 (7) (2010) 410–418. [8] Mark A. Beaumont, Jean-Marie Cornuet, Jean-Michel Marin, Christian P. Robert, Adaptive approximate Bayesian computation, Biometrika (2009) asp052. [9] Michael G.B. Blum, Olivier François, Non-linear regression models for approximate Bayesian computation, Stat. Comput. 20 (1) (2010) 63–73. [10] Yalchin Efendiev, Juan Galvis, Thomas Y. Hou, Generalized multiscale finite element methods (GMsFEM), J. Comput. Phys. (2013). [11] Yalchin Efendiev, Juan Galvis, Guanglian Li, Michael Presho, Generalized multiscale finite element methods: Nonlinear elliptic equations, Commun. Comput. Phys. 15 (03) (2014) 733–755. [12] Eric T. Chung, Yalchin Efendiev, Guanglian Li, An adaptive GMsFEM for high-contrast flow problems, J. Comput. Phys. 273 (2014) 54–76. [13] Eric T. Chung, Yalchin Efendiev, Chak Shing Lee, Mixed generalized multiscale finite element methods and applications, Multiscale Model. Simul. 13 (1) (2015) 338–366. [14] Michel Loève, Probability Theory. II, fourth ed., Springer-Verlag, New York, 1978. [15] Anirban Mondal, Bani Mallick, Yalchin Efendiev, Akhil Datta-Gupta, Bayesian uncertainty quantification for subsurface inversion using a multiscale hierarchical model, Technometrics 56 (3) (2014) 381–392.