Rapid quantification of uncertainty in permeability and porosity of oil reservoirs for enabling predictive simulation

+Model MATCOM-3939; No. of Pages 14

ARTICLE IN PRESS Available online at www.sciencedirect.com

Mathematics and Computers in Simulation xxx (2013) xxx–xxx

Original article

Rapid quantification of uncertainty in permeability and porosity of oil reservoirs for enabling predictive simulation夽 V. Ginting a , F. Pereira b , A. Rahunanthan c,∗,1 b

a Department of Mathematics, University of Wyoming, Laramie, WY 82071, United States School of Energy Resources, Department of Chemical and Petroleum Engineering, Department of Mathematics, University of Wyoming, Laramie, WY 82071, United States c Department of Mathematics and Statistics, University of Toledo, Toledo, OH 43606, United States

Received 17 August 2011; received in revised form 8 October 2012; accepted 5 April 2013

Abstract One of the most difficult tasks in subsurface flow simulations is the reliable characterization of properties of the subsurface. A typical situation employs dynamic data integration such as sparse (in space and time) measurements to be matched with simulated responses associated with a set of permeability and porosity fields. Among the challenges found in practice are proper mathematical modeling of the flow, persisting heterogeneity in the porosity and permeability, and the uncertainties inherent in them. In this paper we propose a Bayesian framework Monte Carlo Markov Chain (MCMC) simulation to sample a set of characteristics of the subsurface from the posterior distribution that are conditioned to the production data. This process requires obtaining the simulated responses over many realizations. In reality, this can be a prohibitively expensive endeavor with possibly many proposals rejection, and thus wasting the computational resources. To alleviate it, we employ a two-stage MCMC that includes a screening step of a proposal whose simulated response is obtained via an inexpensive coarse-scale model. A set of numerical examples using a two-phase flow problem in an oil reservoir as a benchmark application is given to illustrate the procedure and its use in predictive simulation. © 2013 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Dynamic data integration; GPU; Two-phase flows; Uncertainty quantification

1. Introduction Reliable characterization of subsurfaces is one of the most challenging tasks in the flow through porous media community. A common practice is to use the available dynamic data to guide the characterization of a particular reservoir. This practice, which is coined “history matching,” amounts to adjusting the reservoir model until it closely reproduces the recorded data. In reality, history matching requires sampling subsurface characteristics, such as porosity and permeability fields, submitting them as input parameters to the reservoir model (usually expressed in terms of a 夽 This work is partially supported by the grants from DOE (DE-FE0004832 and DE-SC0004982), the Center for Fundamentals of Subsurface Flow of the School of Energy Resources of the University of Wyoming (WYDEQ49811GNTG, WYDEQ49811PER, WYDEQ49811FRTD), 2011 Clean Coal Technologies Research Program of the School of Energy Resources of the University of Wyoming (1100 20352 2012), and from NSF (DMS-1016283). ∗ Corresponding author. Tel.: +1 419 530 2425. E-mail addresses: [email protected] (V. Ginting), [email protected] (F. Pereira), [email protected] (A. Rahunanthan). 1 This work was done while the author was in the Department of Mathematics at the University of Wyoming.

0378-4754/$36.00 © 2013 IMACS. Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.matcom.2013.04.015

Please cite this article in press as: V. Ginting, et al. Rapid quantification of uncertainty in permeability and porosity of oil reservoirs for enabling predictive simulation, Math. Comput. Simul. (2013), http://dx.doi.org/10.1016/j.matcom.2013.04.015

+Model MATCOM-3939; No. of Pages 14 2

ARTICLE IN PRESS

V. Ginting et al. / Mathematics and Computers in Simulation xxx (2013) xxx–xxx

set of governing mathematical principles), obtaining the output of simulated responses, and then comparing them with the measured data. However, direct simulations of flow and transport problems may be computationally expensive or infeasible to compute. Ideally, the level of mesh resolution on which the flow and transport are sought has to be in a comparable scale with the level of subsurface characteristics resolution. Often this would require designing a highly unstructured grid along with the possibility of local refinements in various subregions. Putting this into the scenario of sampling large numbers of subsurface characteristics, one can see the tremendous computational load that has to be committed in history matching. On a related topic, it is not immediately obvious how to rigorously represent the subsurface characteristics before sending them to the reservoir model. In this regard, the predicament lies on the uncertain nature of these characteristics. In our context, the values of these characteristics have to be projected to the underlying computational mesh on which the simulated responses are to be evaluated. This will translate into the dimensional immensity of the uncertainty space which greatly multiplies the already huge computational cost. To recapitulate, establishing a complete statistical description of the subsurface that is consistent with measurement data necessitates focusing the effort on addressing three fundamental problems, namely, (1) proper modeling of the flow and transport in a subsurface, (2) appropriate parametrization of the subsurface characteristics in hopes of representing the uncertainty correctly, and (3) accurate numerical procedures that make the overall computational work tractable. Our thesis in this paper is that these fundamental problems can be resolved by implementing careful modeling reduction techniques. The reduced system is then placed as a black-box in the framework of Bayesian statistical inference in combination with the Markov chain Monte Carlo (MCMC) method. This statistical approach aims at generating a Markov chain from which a stationary, posterior distribution of the subsurface characteristics may be constructed. In our case, the goal is to describe the posterior distribution of the permeability and porosity conditioned to known measurement data. Furthermore, as the uncertainty space of a permeability and porosity field may be exceptionally large, a reduction of the space dimension is a step which must be performed for computationally feasible simulations. The well known Karhunen–Loève expansion [21] allows for parametrization of the random space and thereby leads to the desired reduction. The parametrization is achieved by solving an eigenvalue problem involving a given covariance structure and the assumption of uncorrelated random coefficients. This technique has been used previously in flow through porous media applications (see, e.g. [10,11]), where a random permeability field defined on a large number of underlying grid points is expressed using the expansion. Many authors studied Bayesian methods for inverse problems in reservoir modeling. We refer the reader to [19,23] for recent works in inverse problem in applications of flow through porous media. Lee et al. [19] used an intrinsically stationary Markov random field, which compares favorably to Gaussian process models and offers some additional flexibility and computational advantages, for the choice of prior for the unknown permeability field. Through a Bayesian approach, using MCMC methods to explore the high-dimensional posterior distribution, they demonstrated how to characterize an aquifer based on flow data. In contrast, Efendiev et al. [6,10,11] used the Karhunen–Loève expansion to parametrize the permeability field. In practice however, MCMC can be computationally expensive, particularly for solving forward problems in real flow problems. Higdon et al. [17] presented a methodology for improving the speed and efficiency of an MCMC by combining runs on different scales. They introduced a coarse-scale to make the MCMC chain run faster and better explore the posterior, and linked the coarse chain back to the original fine-scale chain of interest. The proposed coupled MCMC runs more efficiently without sacrificing the accuracy achieved at the fine-scale. In [10], a two-stage MCMC, that utilizes inexpensive coarse-scale models to screen out detailed flow and transport simulations, was used to explore the posterior distribution of permeability field. In the first stage, a new proposal is first tested at the coarse-scale model. If the proposal passes the testing at the coarse-scale model, then at the second stage the fine-scale simulation will be run and this fine-scale run is very expensive compared to the coarse-scale run. We use the Metropolis–Hasting MCMC [16] with a random walk sampler to explore the posterior distribution. To achieve high efficiency in this process, we implement a tuning mechanism so that the sampler has an acceptance rate that is neither too high nor too low. Since we explore a high-dimensional posterior distribution, we implement a component-wise auto-tuning mechanism [14] at the coarse-scale stage in a two-stage MCMC. It has been generally known that permeability and porosity are two of the most important reservoir properties governing the movement and storage of the fluids. To the best of the authors’ knowledge, the current paper is the first Please cite this article in press as: V. Ginting, et al. Rapid quantification of uncertainty in permeability and porosity of oil reservoirs for enabling predictive simulation, Math. Comput. Simul. (2013), http://dx.doi.org/10.1016/j.matcom.2013.04.015

+Model MATCOM-3939; No. of Pages 14

ARTICLE IN PRESS

V. Ginting et al. / Mathematics and Computers in Simulation xxx (2013) xxx–xxx

3

attempt to simultaneously quantify uncertainty associated with permeability and porosity. Mainly because the variation in porosity is typically at least one order of magnitude smaller than the variation in permeability, it has been a routine assumption to take the porosity as a deterministic constant in predictive simulations. Still, others (see for example [22,29]) have investigated possible correlations between permeability and porosity which allows for making them to vary along the reservoir, and these have been largely empirical in nature and were on the case by case basis. Instead of assuming ad-hoc correlations between permeability and porosity fields, we model these fields as independent random fields; possible correlations between them should be identified within our MCMC framework. In this paper we consider two-phase immiscible, incompressible displacement in petroleum reservoirs, and the resulting subsurface characterization from the MCMC step is exploited in two ways. One is on predicting the future production. In this case, we use partial curves of the production wells to characterize the permeability and porosity fields. A set of the accepted proposals is then sent to the numerical model to run the forward problem and predict the rest of the production curves by aggregating the solutions of the forward problem. The other is on devising an informed decision making strategy for placing the wells. This is even more crucial, for example, in the context of toxic geological sequestration. Owing to the computational burden in the MCMC sampling, we parallelize the numerical algorithms for solving the forward problem. Recent investigation on the use of a GPU for solving multiphase flow problems in a three dimensional setting indicates a speed-up of up to 60 in relative comparison to the traditional serial calculation on a CPU [24]. Our strategy in the current work is to exploit the advances in GPU computing by using a heterogeneous CPU-GPU system to solve the forward problem. This paper is organized as follows. We discuss the physical modeling of the problem at hand in Section 2. In Section 3, numerical schemes to solve the flow problem on a fine-grid, and an upscaling procedure to solve the same problem on a coarse-grid are discussed. The parametrization of uncertainty using Karhunen–Loève expansion for unknown permeability and porosity fields is discussed in Section 4. In Section 5, we discuss a Bayesian approach for quantifying uncertainty in both permeability and porosity fields, using two-stage MCMC methods. In Section 6, numerical experiments are presented to demonstrate the power of the approach in reservoir modeling. Section 7 contains conclusions and future work.

2. Physical modeling We consider a heterogeneous oil reservoir which is confined in a domain . The reservoir is equipped with an injection well, from which water is discharged to displace the trapped oil towards the production wells, situated on the perimeter of the domain. See Fig. 1. The dynamics of the movement of the fluids in the reservoir is categorized as an immiscible two-phase system with water and oil (denoted by w and o, respectively) that is incompressible. Capillary pressure is not included in the model. Further simplifying assumptions that we use are gravity-free environment with no sources or sinks (injection

Fig. 1. Description of the problem. (a) Both the production wells are on the perimeter. (b) A production well in the domain and another on the perimeter and (c) Discretization of the domain for configuration A and configuration B.

Please cite this article in press as: V. Ginting, et al. Rapid quantification of uncertainty in permeability and porosity of oil reservoirs for enabling predictive simulation, Math. Comput. Simul. (2013), http://dx.doi.org/10.1016/j.matcom.2013.04.015

+Model MATCOM-3939; No. of Pages 14 4

ARTICLE IN PRESS

V. Ginting et al. / Mathematics and Computers in Simulation xxx (2013) xxx–xxx

and production are modeled through boundary conditions), and the two fluids fill the pore space. Then, the Darcy’s law with a statement of conservation of mass allow us to write the governing equations of the flow as ∇ · v = 0,

where v = −λ(s)k(x)∇p,

(2.1)

and ∂s + ∇ · (f (s)v) = 0, (2.2) ∂t where v is the Darcy velocity, s is the water saturation, k is the permeability and φ is the porosity. The total mobility λ(s) and the flux function f(s) are respectively given by φ(x)

λ(s) =

krw kro + , μw μo

f (s) =

krw (s)/μw , λ(s)

(2.3)

where krj , j = w, o, is the relative permeability  of the phase j [4]. In the simulation we have used boundary conditions v · n = gN on N and p = 0 on D , with N D = ∂, and for a given flux gN . We consider sampling of two important physical quantities, permeability and porosity, of the oil reservoir using partial fractional curves as governing measurements. For each production boundary the fractional flow F(t) is defined as the fraction of oil in the produced fluid, i.e.,  ∂ vn f (s)dl F (t) = 1 −  out , (2.4) ∂out vn dl where ∂out denotes outflow boundary, vn is normal velocity field and t in dimensionless time PVI, which is computed as   T PVI = (2.5) Vp −1 vn dl dτ, 0

∂out

where Vp is the total pore-volume of the reservoir and T denotes the time taken for injection of water. 3. Numerical modeling At the core of the numerical approximation is the partition of  into non-overlapping rectangular elements of Ii,j . With xi = (i − 1/2)hx , yj = (j − 1/2)hy , xi±1/2 = xi ± 0.5hx and yj±1/2 = yj ± 0.5hy , where i = 1, . . ., N and j = 1, . . ., M with N and M denote the number of elements/cells in x- and y-direction, respectively, we define Ii,j : = [xi−1/2 , xi+1/2 ] × [yj−1/2 , yj+1/2 ]. See the right most picture of Fig. 1. We employ an operator splitting technique to compute the solutions of the pressure (2.1) and saturation Eq. (2.2), which are coupled together by the Darcy velocity v and water saturation s. Typically, for computational efficiency, larger time steps are used to compute the pressure Eq. (2.1). The splitting technique allows for using larger time steps in the pressure calculation than the steps allowed under an appropriate CFL condition in the saturation calculation. We thus introduce a variable time step ts for the saturation calculation, and a longer time step tp for the pressure calculation. The pressure and thus the seepage velocity are approximated at times tm = mtp , where m = 0, 1, . . .; the  saturation is approximated at times t m,k = t m + ki=1 (ts )i for tm < tm,k ≤ tm+1 . We remark that we must specify the water saturation at t = 0. For further detail on the operator splitting procedure, we refer the reader to [1,2,8]. In the simulation of (2.1) and (2.2), maintaining the divergence-free property of the Darcy’s velocity is very important. This is precisely the reason for employing a mixed finite element method to solve (2.1). To this end, let U g = {v ∈ H(div, ) : v · n|N = g} and we set U h,g ⊂ U g and W h ⊂ L2 () to be the lowest-order Raviart–Thomas space [26] over rectangles. The mixed finite element approximation consists of finding (vh , ph ) ∈ U h,gN × W h such that (λ−1 (s)k−1 (x)vh , u) − (ph , ∇ · u) = 0, (∇ · vh,w ) = 0,

∀ w ∈ W h.

∀ u ∈ U h,0 ,

(3.1)

The resulting algebraic system is then efficiently solved in a GPU by the conjugate gradient method preconditioned by the algebraic multi-grid scheme [20]. The approximate pressure ph is computed at the center of each cell Ii,j while the normal velocity vh · n is computed at the four edges of the cell Ii,j . Please cite this article in press as: V. Ginting, et al. Rapid quantification of uncertainty in permeability and porosity of oil reservoirs for enabling predictive simulation, Math. Comput. Simul. (2013), http://dx.doi.org/10.1016/j.matcom.2013.04.015

+Model MATCOM-3939; No. of Pages 14

ARTICLE IN PRESS

V. Ginting et al. / Mathematics and Computers in Simulation xxx (2013) xxx–xxx

5

To solve (2.2), we use a conservative second-order, high-resolution central scheme, which can be viewed as a generalization of Kurganov–Tadmor central scheme [18]. The scheme is based on the assumption that an approximation of the solution at time level t = tn is already available in the form of  n s˜ i,j sh (x, t n ) := (x)χi,j (x), with s˜ i,j (x) = si,j + ∇si,j · ri,j , (3.2) i,j

where ∇si,j is the gradient of sh approximated by MinMod limiter, ri,j = [(x − xi ) (y − yj )] , χi,j is the characteristic function of Ii,j , and si,j is the average of sh over Ii,j . With the above definitions, in the semi-discrete formulation, the evolution step is performed by integrating over non-equal control volumes, whose size are determined by a local speed of propagation. The evolved solution is then projected back to the original grid with a piecewise polynomial reconstruction. The resulting non-staggered fully discrete scheme is passed through t → 0, while keeping hx and hy fixed to yield y

φi,j

y

x x (t) − Hi−1/2,j (t) Hi,j+1/2 (t) − Hi,j−1/2 (t) Hi+1/2,j d si,j (t) = − − , dt hx hy

(3.3)

with the numerical flux x Hi+1/2,j = (vh · n)i+1/2,j {f (sh )}i+1/2,j −

 1  ai+1/2,j  φi+1/2,j [sh ]i+1/2,j , 2

and the rest of the terms are expressed similarly. Here the index pair (i + 1/2, j) refers to the center (xi+1/2 , yj ) of the right edge of Ii,j , so the subscript signifies the evaluation of the quantities at this location: (vh · n)i+1/2,j is the normal component of vh , {f(sh )}i+1/2,j is the average of f(sh ) over Ii,j and Ii+1,j cross-weighted by their respective φ, ai+1/2,j is the local maximum speed, φi+1/2,j is the harmonic average of φ over Ii,j and Ii+1,j , and [sh ]i+1/2,j is the jump of sh . The resulting semi-discrete scheme can then be solved by a SSP Runge–Kutta scheme [28]. We refer the reader to [24,25] for further discussion on the numerical schemes used in solving the forward problem that takes advantage of the GPU computing. 4. Parametrization of uncertainty As the uncertainty space describing the permeability and porosity can be exceptionally large, a reduction of the space dimension is a step that must be performed for computationally feasible simulations. The Karhunen–Loève (KL) expansion [21,30] allows for the parametrization of the uncertainty space and thereby leads to the desired reduction. We note that this reduction technique has been applied within Bayesian MCMC in [6,10,11,15,23]. The KL expansion relies on a basic prior information on the structure of the permeability and porosity fields that reasonably match the actual reservoir. In practice however, this is hindered by the very limited information available to use. A standard assumption in geostatistics is to model the permeability to follow a log-normal distribution [5], i.e., log [k(x, ω)] = Y k (x, ω), where x = (x, y) ∈  ∈ R2 , and ω is a random element in a probability space, with Yk (x, ω) a field possessing a Gaussian distribution and a covariance function  2 2 − x | − y | |x |y 1 2 1 2 . R(x1 , x2 ) = σY2 exp − − (4.1) 2L2x 2L2y 2

Here Lx and Ly are correlation lengths in x- and y-direction, respectively, and σY2 = E[(Y k ) ]. We now briefly describe the essentials of the KL expansion. Suppose Yk (x, ω) is a second-order stochastic process, that is, Yk (x, ω) ∈ L2 () with a probability of one. We will assume that E[Yk (x, ω)] = 0. Given an arbitrary orthonormal basis {ϕi } in L2 , we can expand Yk (x, ω) as Y (x, ω) = k

∞  i=1

 Yik (ω)ϕi (x),

withYik (ω)

=

Y k (x, ω)ϕi (x)dx

(4.2)



Please cite this article in press as: V. Ginting, et al. Rapid quantification of uncertainty in permeability and porosity of oil reservoirs for enabling predictive simulation, Math. Comput. Simul. (2013), http://dx.doi.org/10.1016/j.matcom.2013.04.015

+Model MATCOM-3939; No. of Pages 14 6

ARTICLE IN PRESS

V. Ginting et al. / Mathematics and Computers in Simulation xxx (2013) xxx–xxx

functions of a random variable. We are interested in the special basis in L2 () which allows Yik to be uncorrelated, E[Yik Yjk ] = 0 for all i = / j. Because by definition R(x1 , x2 ) = E[Yk (x1 )Yk (x2 )], such basis functions {ϕi } satisfy   E[Yik Yjk ] = ϕi (x1 )dx1 R(x1 , x2 )ϕj (x2 )dx2 = 0, i = / j. (4.3) 



Completeness of L2 () implies that ϕi (x) is an eigenfunction of integral equation involving R(x1 , x2 ) expressed as  R(x1 , x2 )ϕi (x2 )dx2 = λi ϕi (x1 ), i = 1, 2, . . . , (4.4) 

√ 2 where λi = E[(Yik ) ] > 0. Denoting θik = Yik / λi , then θik satisfies E(θik ) = 0 and E(θik θjk ) = δij , and thus Y (x, ω) = k

∞



λi θik (ω)ϕi (x),

(4.5)

i=1

where ϕi and λi satisfy (4.4). We assume that eigenvalues λi are ordered so that λ1 ≥ λ2 ≥ . . .. The expansion (4.5) is called the Karhunen–Loève expansion. In (4.4), the L2 basis functions ϕi (x) are deterministic and resolve the spatial dependence of the permeability field. The uncertainty is represented by the scalar random variables θik . In general, we only need to keep the leading order terms (quantified by the magnitude of λi ) and still capture most of the energy of the stochastic process Yk (x, ω). For a Nk -term KL expansion YNk k =

Nk 

λi θik ϕi .

(4.6)

i=1

If λi decays very fast, then the truncated KL expansion (4.6) is a good approximation of the stochastic process in the L2 sense. We remark that one could use a different covariance function, such as for example,  |x1 − x2 | |y1 − y2 | 2 . − (4.7) R(x1 , x2 ) = σY exp − 2Lx 2Ly However, the eigenvalues of KL expansion may decay algebraically [12,10]. To have a good approximation for (4.5), we must add more terms in the truncated KL expansion (4.6). This will increase the parameter space dimension, and thus make the sampling of permeability and porosity more expensive. We first solve the eigenvalue problem (4.4) numerically (typically in coarser grids than those used for the approximation of the (2.1 – 2.2)) and obtain the eigenpairs (λi , ϕi ). The truncated KL expansion approximates the stochastic process Yk (x, ω), and sampling Yk (x, ω) is proceeded by generating Gaussian random variables θik . With respect to the porosity field, we model it by assuming that the permeability and porosity fields share the same spatial distribution (we do not use ad-hoc correlations between these fields). In turn, this allows us to employ the Karhunen–Loève basis functions (4.4). As in the KL expansion for the permeability field, we use a Nφ -term KL expansion φ YNφ (x, ω)

=

Nφ 

φ

λi θi (ω)ϕi (x).

(4.8)

i=1

Then we define the porosity as φ(x) =

φmin + φmax e 1+e

φ YN φ

φ

YN

φ

,

φmin and φmax ∈ (0, 1),

(4.9)

where φmin and φmax are the lower and upper limits of the porosity of the reservoir. Please cite this article in press as: V. Ginting, et al. Rapid quantification of uncertainty in permeability and porosity of oil reservoirs for enabling predictive simulation, Math. Comput. Simul. (2013), http://dx.doi.org/10.1016/j.matcom.2013.04.015

+Model MATCOM-3939; No. of Pages 14

ARTICLE IN PRESS

V. Ginting et al. / Mathematics and Computers in Simulation xxx (2013) xxx–xxx

7

5. A Bayesian MCMC framework 5.1. Posterior exploration Our problem of interest is to sample the permeability and porosity fields conditioned on the available fractional flow data with measurement errors. This is translated into sampling from the conditional distribution P(ψ|F), where ψ = [θ k , θ φ ] with θ k and θ φ vectors containing the random coefficients in the KL expansions (4.6) and (4.8), respectively. According to Bayes’ theorem this distribution satisfies the proportionality relation π(ψ) = P(ψ|Fm ) ∝ P(Fm |ψ)P(ψ),

(5.1)

where P(Fm |ψ) represents the likelihood function (that requires the forward solution of the two-phase flow) and P(ψ) is the prior distribution of ψ. The normalizing constant in this expression is not important because we use an iterative updating procedure. We assume that the likelihood function follows a Gaussian distribution. i.e., P(Fm |ψ) ∝ exp( − Fm − Fψ 2 ),

(5.2)

where Fm are the reference fractional flow curves, F␺ are the simulated fractional flow curves that is obtained by solving the forward problem with known permeability k and porosity φ, in other words with the proposed ψ, and  T Fm − Fψ 2 = Σ(t)(Fm (t) − Fψ (t))2 dt, (5.3) 0

with  the covariance function representing the measurement errors. We take (t) = 1/2σF2 , where σF2 is the precision associated with the measurement Fm and numerical solution F␺ . In practice, both porosity and permeability might be dependent on each other, and that is the reason we take the same KL expansion structure for both fields. However, when we make proposal, we assume that θ k and θ φ are independent of each other. We use the Metropolis–Hasting MCMC to sample from the posterior distribution. The main idea of MCMC is to generate a Markov chain whose stationary distribution is given by π(ψ). At each iteration, θ k and θ φ are proposed using an instrumental distribution q(ψp |ψ), where ψ and ψp represent the previously accepted proposal and the current proposal, respectively. The forward problem is then solved to determine the acceptance probability,  q(ψ|ψp )P(ψp |Fm ) α(ψ, ψp ) = min 1, (5.4) , q(ψp |ψ)P(ψ|Fm ) i.e., ψp is accepted with the probability α(ψ, ψp ). This single-stage MCMC chain may burn within a few thousand cycles after initialization, and an additional several thousand cycles are needed to explore the posterior in real flow problems. In most reservoir simulations, solving the fine-scale forward problem is computationally expensive. Thus, the single-stage MCMC limits the exploration of the posterior in a practical time. Therefore, we use a coarse-scale solution based on upscaling before running the fine-scale flow simulator. Here, we run the flow simulator on the coarse-scale model and compare the fractional flow curves to determine if we need to run on the fine-scale model. 5.2. A two-stage MCMC method As alluded to earlier, we propose to develop a two-stage MCMC in which a screening process is implemented using a coarse-scale model approximating the original Eqs. (2.1) and (2.2). First, we proceed by giving a brief description of this coarse-scale model. It is solved on a coarse-grid whose discretization is done in a similar fashion as in the fine-scale calculation described in Section 3, but with discretization parameters significantly larger than their finescale counterpart. The key lies on a rigorous representation of k and φ on the coarse-grid obtained from the fine-grid resolution. In this case, we use an upscaling procedure such that effective permeability and porosity fields on the coarse-grid deliver the same average response as that of the underlying fine-scale problem locally. Please cite this article in press as: V. Ginting, et al. Rapid quantification of uncertainty in permeability and porosity of oil reservoirs for enabling predictive simulation, Math. Comput. Simul. (2013), http://dx.doi.org/10.1016/j.matcom.2013.04.015

+Model MATCOM-3939; No. of Pages 14 8

ARTICLE IN PRESS

V. Ginting et al. / Mathematics and Computers in Simulation xxx (2013) xxx–xxx

To compute the effective permeability on the coarse-grid, we solve for ζ governed by equations similar to (2.1) subject to the constant pressure and no flux boundary conditions on each coarse-grid cell I. After solving these equations on the coarse-grid, we compute the effective permeability tensor as [9]  1 (k(x)ei , ej ) = (5.5) (k(x)∇ζi (x), ej )dx, |I| I where ei is the standard Euclidean basis vector. Once we know the effective permeability on the coarse-grid cells, we can solve (2.1) on the coarse-grid with k replacing k to get the upscaled Darcy’s velocity. In a similar fashion, we solve (2.2) on the coarse-grid with an upscaled porosity computed as the arithmetic average over each I. It is not difficult to see that the numerical solution of the coarse-scale model is relatively inexpensive compared to the finescale calculation and yet it manages to capture the general trend of the process. Hence it is suitable for the screening purpose mentioned above. We remark that the upscaling procedure may produce a systematic bias in the numerical results [3,27]. The coarse-grid calculations should be computationally inexpensive, but not necessarily very accurate. However, for improved approximations on coarse-grids one could use multi-scale methods (see for example [13]). Algorithm 5.1. Two-stage MCMC Given covariance function R, generate KLE parametrization {λn , ϕn }N n=1 for p = 1 to Mmcmc do φ (1) At ψ = (θ k , θ φ ), generate ψp = (θ kp , θ p ) from q(ψp |ψ) (2) With ψp , compute k and φ using KLE and solve

the forward problem  on the coarse-scale to get Fc (3) Accept ψp with probability αc (ψ, ψp ) = min 1,

q(ψ|ψp )Pc (ψp |Fm ) q(ψp |ψ)Pc (ψ|Fm )

(4) If ψp is rejected, go back to (1), else use ψp in the fine-scale simulation  to get Ff (5) Take ψp with probability αf (ψ, ψp ) = min 1,

Pf (ψp |Fm )Pc (ψ|Fm ) Pf (ψ|Fm )Pc (ψp |Fm )

end for

The basic procedure of the two-stage MCMC is shown in Algorithm 5.1 [11]. We consider the random walk sampler, L L L L L θL p = θ + γ , where L represents k or φ, θ is the previously accepted proposal, θ p is the current proposal, γ is the k φ tuning parameter and  is a set of Gaussian random variables. Here, we treat θ and θ separately. Since we explore the high-dimensional posterior distribution, we use a component-wise auto-tuning mechanism in MCMC with a target acceptance rate of 0.44. The auto-tuning mechanism is implemented to select the tuning parameter γ automatically in the MCMC chain, where γ = [γ 1 , γ 2 , . . ., γ n ] for n dimensional posterior distribution. In this mechanism, we first assign an initial value for each tuning parameter γ i and run the random walk sampler for a specified number of proposals (say 100 proposals). We then calculate the acceptance rate for each component i. If there are any components whose acceptance rates fall outside 0.44, the corresponding γ i s are adjusted, and the MCMC chain is continuously run without interrupting the chain for another 100 proposals. We continue in this manner throughout the posterior exploration aiming at achieving optimum acceptance rates for all components [14]. For the single-stage MCMC, implementing an auto-tuning mechanism is straightforward. However, for the two-stage MCMC, we can implement an auto-tuning mechanism at the coarse-scale or fine-scale stage. The more (or less) the number of acceptance at the coarse-scale stage, the more (or less) the number of acceptance at the fine-scale stage. In other words, we can control the acceptance at the fine-scale stage by controlling the acceptance at the coarse-scale stage. One could also tune by implementing the auto-tuning at the fine-scale stage. This, in turn, will control the acceptance at the coarse-scale. Since the former naturally follows the flow of two-stage MCMC concept, we adapt the tuning at the coarse-scale stage. 6. Numerical results We now discuss the simulations of the two-phase flow problem in an oil reservoir as illustrated in Fig. 1 and present the associated numerical results. We note that in the numerical implementation all the wells are modeled through boundary conditions. At the injection wells, we specify the flux of the water phase. We have two production wells: one well is situated along the diagonal, opposite to the injection well, and the other is situated at the center of a side which is one of the two sides that enclose the production well at the corner. Let us call this configuration as configuration A. In this configuration, we refer to them as corner and perimeter wells, respectively (see Fig. 1a). In configuration B, one well is situated along the diagonal, opposite to the injection well, and the other is situated at the center of the domain. Please cite this article in press as: V. Ginting, et al. Rapid quantification of uncertainty in permeability and porosity of oil reservoirs for enabling predictive simulation, Math. Comput. Simul. (2013), http://dx.doi.org/10.1016/j.matcom.2013.04.015

+Model MATCOM-3939; No. of Pages 14

ARTICLE IN PRESS

V. Ginting et al. / Mathematics and Computers in Simulation xxx (2013) xxx–xxx

9

Fig. 2. Top: Water saturation plots for configuration A. Bottom: Water saturation plots for configuration B. Left: At 0.4 PVI. Right: 1.3 PVI.

Here, we refer to them as corner and center wells, respectively (see Fig. 1b). Fig. 2 shows the water saturation plots at t = 0.4 PVI and t = 1.3 PVI for configuration A and configuration B. Spatial profiles of permeability and porosity fields for the reservoir of interest are given in Fig. 4. The relative permeability functions of water and oil take the form of s2 and (1 − s)2 , respectively, and the viscosity ratio between water and oil is 1:20. We assume that at t = 0, the reservoir is saturated by oil without any water, i.e., s = 0. The water is then injected at the injection well at the rate of one pore-volume every five years. We assume that we have observed production curves for the wells until 0.4 PVI. Therefore, we run the MCMC until 0.4 PVI. For KLE in (4.1), we select the correlation length Lx = Ly = 0.2 and variance σ Y 2 = 4. Fig. 3 shows that the eigenvalues decay very fast for these values, and it is enough to consider the first twenty eigenvalues in the KLE. Since we assume that the permeability and porosity share the same spatial structure, we share the same KLE structure for the permeability and porosity fields, with Nk = Nφ = 20. Next, we turn our attention to σF2 in the likelihood function (5.2). The σF2 must be fixed a priori [19], because treating σF2 as an unknown parameter results in an unacceptably large estimate. The smaller the value of σF2 , the better the sampled production curves are. If σF2 is chosen to be too large, then the posterior will ignore the likelihood and simply draw from the prior; if σF2 is too small, then the posterior probability exists only in extreme modes, which unrealistically restricts possible permeability and porosity configurations. It is ideal to specify σF2 from knowledge of the error in measuring permeabilities and porosities. In our application, we specify σF2 = 10−4 . For the simulation, we use a fine grid of size 128 × 128, while for the two-stage MCMC, we use a coarse-grid of 32 × 32. In the random walk Metropolis–Hastings sampler, we implement the component-wise tuning for (θ k , θ φ ) based on the coarse-scale acceptance rate in the two-stage MCMC. Table 1 shows comparison of the acceptance rates. As indicated in [10,11,23], in the case of distinct model problems, the two-stage MCMC runs fewer fine-scale black-box simulations than the single-stage MCMC. In Fig. 4, we illustrate the reference permeability and porosity fields along with pairs of accepted fields (for permeability and porosity) for each well configuration considered here. After selecting a set of accepted realizations for permeability and porosity, we run the forward problem until 1.3 PVI. We then take the weighted average of the results of the forward problem. The weight for a particular state ψ in the MCMC chain is determined as p, if the chain stays at that state for (p − 1) proposals. This average curve is referred to as the prediction of the production curve. The prediction curves are plotted using 500, 1000, 1500 and 2000 accepted Please cite this article in press as: V. Ginting, et al. Rapid quantification of uncertainty in permeability and porosity of oil reservoirs for enabling predictive simulation, Math. Comput. Simul. (2013), http://dx.doi.org/10.1016/j.matcom.2013.04.015

+Model MATCOM-3939; No. of Pages 14 10

ARTICLE IN PRESS

V. Ginting et al. / Mathematics and Computers in Simulation xxx (2013) xxx–xxx 0.9

i

0.6

0.3

0.0 0

10

20

30

40

50

index i Fig. 3. Eigenvalues of the KLE for the Gaussian covariance with Lx = Ly = 0.2 and σY2 = 4.

Fig. 4. Characterization of underlying fields. Left: Permeability (in logarithmic scale). Right: Porosity. Top: The reference fields. Middle: A set of accepted fields for configuration A. Bottom: A set of accepted fields for configuration B.

Please cite this article in press as: V. Ginting, et al. Rapid quantification of uncertainty in permeability and porosity of oil reservoirs for enabling predictive simulation, Math. Comput. Simul. (2013), http://dx.doi.org/10.1016/j.matcom.2013.04.015

+Model MATCOM-3939; No. of Pages 14

ARTICLE IN PRESS

V. Ginting et al. / Mathematics and Computers in Simulation xxx (2013) xxx–xxx

11

Table 1 A comparison of accepted proposals for single- and two-stage MCMCs. Config.

MCMC

Proposals

σF2

σC2

Coarse acc.

Fine acc.

Acc. rate

A

Single-stage Two-stage Single-stage Two-stage

9095 16,024 8064 14,003

10−4 10−4 10−4 10−4

– 4 × 10−4 – 4 × 10−4

– 6496 – 5554

2000 2000 2000 2000

22% 31% 25% 36%

B

1.0 0.8 0.6

0.6

0.4

0.4

0.2

0.2

0.0 0.0

0.2

0.4

0.6

0.8

1.0

1.2

reference single-stage two-stage

0.8

F(t)

F(t)

1.0

reference single-stage two-stage

0.0 0.0

1.4

0.2

time (t in PVI)

1.0

0.2

0.2

0.6

0.8

1.0

1.2

0.0 0.0

1.4

0.2

0.2

0.2

0.6

0.8

1.0

1.2

0.0 0.0

1.4

0.2

0.2

0.2

0.6

0.8

1.0

time (t in PVI)

0.8

1.0

1.2

1.4

1.2

1.4

reference single-stage two-stage

0.6 0.4

0.4

0.6

0.8

0.4

0.2

0.4

1.0

F(t)

F(t)

0.6

0.0 0.0

1.4

time (t in PVI)

reference single-stage two-stage

0.8

1.2

reference single-stage two-stage

time (t in PVI)

1.0

1.0

0.6 0.4

0.4

0.8

0.8

0.4

0.2

0.6

1.0

F(t)

F(t)

0.6

0.0 0.0

0.4

time (t in PVI)

reference single-stage two-stage

0.8

1.4

reference single-stage two-stage

time (t in PVI)

1.0

1.2

0.6 0.4

0.4

1.0

0.8

0.4

0.2

0.8

1.0

F(t)

F(t)

0.6

0.0 0.0

0.6

time (t in PVI)

reference single-stage two-stage

0.8

0.4

0.0 0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

time (t in PVI)

Fig. 5. A comparison of single- and two-stage MCMCs with σF2 = 10−4 for configuration A. From top to bottom, 500, 1000, 1500 and 2000 accepted realizations used for prediction, respectively. Left: The perimeter well. Right: The corner well.

Please cite this article in press as: V. Ginting, et al. Rapid quantification of uncertainty in permeability and porosity of oil reservoirs for enabling predictive simulation, Math. Comput. Simul. (2013), http://dx.doi.org/10.1016/j.matcom.2013.04.015

+Model MATCOM-3939; No. of Pages 14 12

ARTICLE IN PRESS

V. Ginting et al. / Mathematics and Computers in Simulation xxx (2013) xxx–xxx 1.0

1.0

reference single-stage two-stage

0.6

0.6

0.4

0.4

0.2

0.2

0.0 0.0

0.2

0.4

0.6

0.8

1.0

1.2

reference single-stage two-stage

0.8

F(t)

F(t)

0.8

0.0 0.0

1.4

0.2

reference single-stage two-stage

0.6

0.4

0.2

0.2

0.2

0.4

0.6

0.8

1.0

1.2

0.0 0.0

1.4

0.2

0.6

0.8

1.0

0.4

0.2

0.2

0.6

0.8

1.0

1.2

0.0 0.0

1.4

0.2

0.4

0.6

0.8

1.0

1.0

reference single-stage two-stage

0.8

0.6

0.4

0.4

0.2

0.2

0.4

0.6

0.8

1.0

time (t in PVI)

1.4

1.2

1.4

reference single-stage two-stage

0.8

F(t)

0.6

0.2

1.2

time (t in PVI)

time (t in PVI)

1.0

1.4

0.6

0.4

0.4

1.2

reference single-stage two-stage

0.8

F(t)

0.6

0.2

1.4

time (t in PVI)

reference single-stage two-stage

0.8

F(t)

0.4

1.0

1.0

1.2

reference single-stage two-stage

time (t in PVI)

F(t)

1.0

0.6

0.4

0.0 0.0

0.8

0.8

F(t)

F(t)

0.8

0.0 0.0

0.6

1.0

1.0

0.0 0.0

0.4

time (t in PVI)

time (t in PVI)

0.0 0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

time (t in PVI)

Fig. 6. A comparison of single- and two-stage MCMCs with σF2 = 10−4 for configuration B. From top to bottom, 500, 1000, 1500 and 2000 accepted realizations used for prediction, respectively. Left: The center well. Right: The corner well.

fine-scale realizations and given in Figs. 5 and 6. The vertical line separates the measured production curves and the predicted production curves. Both figures show that the two-stage MCMC performs better than the single-stage MCMC for a limited number of MCMC runs. This may be due to the fact that in the single-stage MCMC the chain moves with the smaller steps which is attributed to the auto-tuning mechanism in the MCMC. In the ideal setting, however, the two-stage MCMC can be thought of as a special case of the single-stage MCMC and therefore both should converge to the same steady state distribution. If we study the water saturation plots in Fig. 2, we can observe that both the partial production curves in the configuration B come from the same flow channel. However, the partial production curves in the configuration A come Please cite this article in press as: V. Ginting, et al. Rapid quantification of uncertainty in permeability and porosity of oil reservoirs for enabling predictive simulation, Math. Comput. Simul. (2013), http://dx.doi.org/10.1016/j.matcom.2013.04.015

+Model MATCOM-3939; No. of Pages 14

ARTICLE IN PRESS

V. Ginting et al. / Mathematics and Computers in Simulation xxx (2013) xxx–xxx

13

from two different flow channels. That is, the partial production curves in the configuration A capture more information about the underlying field than the curves in the configuration B. As a result, the prediction curves for the configuration A are more reliable than the prediction curves for the configuration B. It highlights the fact that the location of production wells may play an important role in predicting the flow pattern. We remark that although the predictions reported in Fig. 5 are very accurate, the results of Fig. 6 could be improved by adding sparse measurements of permeability and porosity (at well locations) or more wells to the MCMC framework described here. A systematic study of the relative importance (for prediction) of the number and position of wells and the number of sparse measurements is currently under way and will appear elsewhere. 7. Concluding remarks In this paper we presented a Bayesian framework for the characterization of permeability and porosity fields of a reservoir model and for predicting production based on the characterization. We considered a two-phase problem in an oil reservoir. The partial production curves were used to characterize the reservoir, and then, using the permeability and porosity fields from the characterization, we predicted the production curves of the production wells. Numerical results indicate that this Bayesian approach with the two-stage MCMC can effectively be used to characterize the underlying fields of the reservoir model and to predict the future production curves based on the available production data. By considering two different spatial locations for a production well in an oil field, we showed that the location of the production well may significantly affect the characterization of the underlying field of the reservoir and, consequently, the prediction of the production. This suggests that we must pay attention in selecting the location of the production wells if we wish to characterize the underlying geology of a reservoir. This Bayesian approach with the two-stage MCMC can be applied to multiphase flows for prediction of migration and trapping of CO2 plumes in the subsurface. In the near future we shall focus our attention on the characterization of aquifer heterogeneity and prediction of the motion of CO2 plumes in a simplified CO2 /brine two-phase flow model in deep saline aquifers. We will make use of a recently developed high-performance simulator [7] in this investigation. On a related matter, having an informed decision making strategy for placing monitoring wells is even more crucial in CO2 geological sequestration. In particular, the monitoring of the migration of injected CO2 plumes in aquifers requires that monitoring wells be drilled in such a way that the measured data captures the relevant features of the underlying multi-scale heterogeneity. Finally, we wish to mention a potential byproduct of the characterization framework proposed here: given a set of accepted permeability and porosity realizations we can estimate directly the correlation structure of the permeability and porosity fields. We expect to speed-up considerably our MCMC scheme by taking such correlations into account. We shall pursue this line of investigation in more detail in the future. References [1] E. Abreu, J. Douglas Jr., F. Furtado, F. Pereira, Operator splitting based on physics for flow in porous media, International Journal of Computational Science 2 (3) (2008) 315–335. [2] E. Abreu, J. Douglas Jr., F. Furtado, F. Pereira, Operator splitting for three-phase flow in heterogeneous porous media, Communications in Computational Physics 6 (1) (2009) 72–84. [3] A. Boschan, B. Noetinger, Scale dependence of effective hydraulic conductivity distributions in 3D heterogeneous media: a numerical study, Transport in Porous Media (2012) 1–21. [4] Z. Chen, G. Huan, Y. Ma, Computational Methods for Multiphase Flows in Porous Media, SIAM, Philadelphia, PA, 2006. [5] G. Dagan, Flow and Transport in Porous Formations, Springer-Verlag, 1989. [6] C. Douglas, Y. Efendiev, R. Ewing, V. Ginting, R. Lazarov, Dynamic data driven simulations in stochastic environments, Computing 77 (4) (2006) 321–333. [7] C. Douglas, F. Furtado, V. Ginting, M. Mendes, F. Pereira, M. Piri, On the development of a high-performance tool for the simulation of CO2 injection into deep saline aquifers, Rocky Mountain Geology 45 (2) (2010) 151–161. [8] J. Douglas Jr., F. Furtado, F. Pereira, On the numerical simulation of waterflooding of heterogeneous petroleum reservoirs, Computational Geosciences 1 (1997) 155–190. [9] L. Durlofsky, Numerical calculation of equivalent grid block permeability tensors for heterogeneous porous media, Water Resources Research 27 (5) (1991) 699–708. [10] Y. Efendiev, A. Datta-Gupta, V. Ginting, X. Ma, B. Mallick, An efficient two-stage Markov chain Monte Carlo method for dynamic data integration, Water Resources Research 41 (2011) W12423.

Please cite this article in press as: V. Ginting, et al. Rapid quantification of uncertainty in permeability and porosity of oil reservoirs for enabling predictive simulation, Math. Comput. Simul. (2013), http://dx.doi.org/10.1016/j.matcom.2013.04.015

+Model MATCOM-3939; No. of Pages 14 14

ARTICLE IN PRESS

V. Ginting et al. / Mathematics and Computers in Simulation xxx (2013) xxx–xxx

[11] Y. Efendiev, T. Hou, W. Luo, Preconditioning Markov chain Monte Carlo simulations using coarse-scale models, SIAM Journal on Scientific Computing 28 (2) (2006) 776–803. [12] P. Frauenfelder, C. Schwab, R. Todor, Finite elements for elliptic problems with stochastic coefficients, Computer Methods in Applied Mechanics and Engineering 194 (2005) 205–228. [13] F. Furtado, V. Ginting, F. Pereira, M. Presho, Operator splitting multiscale finite volume element method for two-phase flow with capillary pressure, Transport in Porous Media 90 (3) (2011) 927–947. [14] D. Gamerman, H. Lopes, Markov Chain Monte Carlo – Stochastic Simulation for Bayesian Inference, Vol. 68 of Texts in Statistical Science., 2nd ed., Chapman & Hall, CRC, Boca Raton, 2006. [15] 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, Computational Geosciences 15 (4) (2011) 691–707. [16] W.K. Hastings, Monte Carlo sampling methods using Markov chains and their applications, Biometrika 57 (1970) 97–109. [17] D. Higdon, H. Lee, Z. Bi, A Bayesian approach to characterizing uncertainty in inverse problems using coarse and fine-scale information, IEEE Transactions on Signal Processing 50 (2002) 389–399. [18] A. Kurganov, E. Tadmor, New high-resolution central schemes for nonlinear conservation laws and convection–diffusion equations, Journal of Computational Physics 160 (2000) 241–282. [19] H. Lee, D. Higdon, Z. Bi, M. Ferreira, M. West, Markov random field models for high-dimensional parameters in simulations of fluid flow in porous media, Technical Reports, Technometrics (2002). [20] M. Liebmann, Efficient PDE solvers on modern hardware with applications in medical and technical sciences, University of Graz, 2009 (Ph.D. Thesis). [21] M. Loève, Probability Theory, Springer, Berlin, 1977. [22] S. Ma, N. Morrow, Relationships between porosity and permeability for porous rocks, in: SCA Conference Paper 9610, 1996. [23] X. Ma, M. Al-Harbi, A. Datta-Gupta, Y. Efendiev, An efficient two-stage sampling method for uncertainty quantification in history matching geological models, SPE Journal (2008 March) 77–87. [24] F. Pereira, A. Rahunanthan, Numerical simulation of two-phase flows on a GPU, in: 9th International Meeting on High Perfomance Computing for Computational Science (VECPAR ‘10), Berkeley, CA, 2010. [25] 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 (10) (2011) 2296–2306. [26] P. Raviart, J. Thomas, A mixed finite element method for 2-nd order elliptic problems, in: I. Galligani, E. Magenes (Eds.), Mathematical Aspects of Finite Element Methods, Vol. 606 of Lecture Notes in Mathematics, Springer, Berlin/Heidelberg, 1977, pp. 292–315. [27] R. Romeu, B. Noetinger, Calculation of internodal transmissivities in finite difference models of flow in heterogeneous porous media, Water Resources Research 31 (4) (1995) 943–959. [28] C.-W. Shu, S. Osher, Efficient implementation of essentially nonoscillatory shock-capturing schemes II, Journal of Computational Physics 83 (1) (1989) 32–78. [29] J. Sperl, J. Trckova, Permeability and porosity of rocks and their relationship based on laboratory testing, Acta Geodynamica et Geomaterilia 5 (1) (2008) 41–47. [30] E. Wong, Stochastic Processes in Information and Dynamical Systems, McGraw-Hill, New York, 1971.

Please cite this article in press as: V. Ginting, et al. Rapid quantification of uncertainty in permeability and porosity of oil reservoirs for enabling predictive simulation, Math. Comput. Simul. (2013), http://dx.doi.org/10.1016/j.matcom.2013.04.015