key: cord-0045988-sja0xbc0
authors: Ali, Alsadig; Al-Mamun, Abdullah; Pereira, Felipe; Rahunanthan, Arunasalam
title: Markov Chain Monte Carlo Methods for Fluid Flow Forecasting in the Subsurface
date: 2020-05-25
journal: Computational Science - ICCS 2020
DOI: 10.1007/978-3-030-50436-6_56
sha: b428b67cbb25184a8b9f7d38c3f234024b3a9218
doc_id: 45988
cord_uid: sja0xbc0

Accurate predictions in subsurface flows require the forecasting of quantities of interest by applying models of subsurface fluid flow with very little available data. In general a Bayesian statistical approach along with a Markov Chain Monte Carlo (MCMC) algorithm can be used for quantifying the uncertainties associated with subsurface parameters. However, the complex nature of flow simulators presents considerable challenges to accessing inherent uncertainty in all flow simulator parameters of interest. In this work we focus on the transport of contaminants in a heterogeneous permeability field of a aquifer. In our problem the limited data comes in the form of contaminant fractional flow curves at monitoring wells of the aquifer. We then employ a Karhunen-Loève expansion to truncate the stochastic dimension of the permeability field and thus the expansion helps reducing the computational burden. Aiming to reduce the computational burden further, we code our numerical simulator using parallel programming procedures on Graphics Processing Units (GPUs). In this paper we mainly present a comparative study of two well-known MCMC methods, namely, two-stage and DiffeRential Evolution Adaptive Metropolis (DREAM), for the characterization of the two-dimensional aquifer. A thorough statistical analysis of ensembles of the contaminant fractional flow curves from both MCMC methods is presented. The analysis indicates that although the average fractional flow curves are quite similar, both time-dependent ensemble variances and posterior analysis are considerably distinct for both methods.

In subsurface characterization and flow forecasting, one could characterize the subsurface using a Bayesian framework. The Bayesian framework consists essentially of a Markov Chain Monte Carlo (MCMC) algorithm, in which we repeatedly solve a flow numerical simulator that models the porous media problem of interest [11] . In this paper we use a single-phase flow numerical simulator that simulates the contaminant transport in an aquifer. In the MCMC algorithm we characterize the heterogeneous permeability field of the aquifer using the simulator. Our C/C++ flow simulator takes advantage of Graphics Processing Units (GPUs), which have the computational capacity to speedup the single-phase flow simulation [12] . The GPU flow simulator runs on a computational grid of the permeability field that is divided into several thousand elements. However, populating each computational element by a random permeability value and changing those values in each MCMC iteration is impracticable. In this case the dimension of the stochastic space is that of the computational domain and the dimension reduction is achieved by a Karhunen-Loève expansion (KLE) [15] . Due to the serial nature of the MCMC algorithm, the computational burden in solving the problem recurrently using the GPU numerical simulator is still huge. This can make the Bayesian framework less attractive for our problem. Parallelizations of the MCMC algorithm to speedup the characterization were considered in [10, 12] . The simulation of several parallel MCMCs reduces the computational cost drastically. However, the convergence of such parallel MCMCs should be carefully analyzed.

In this paper we consider two-stage and DiffeRential Evolution Adaptive Metropolis (DREAM) MCMC methods for the subsurface characterization. The two-stage MCMC was introduced in [6, 8] and it has been investigated more recently (see [7, 14] and references therein). The two-stage procedure is of particular interest to subsurface flow problems because samples can be rejected with inexpensive simulations on coarse grids. The DREAM is an extension of Differential Evolution (DE) MCMC, which integrates the essential ideas of DE genetic algorithm and MCMC algorithm [3] . Inspired by DE MCMC, the DREAM was proposed in [22] and applied to several interesting problems (see [19] [20] [21] and references therein). The DREAM is well known to converge relatively faster when compared to earlier procedures. It employs subspace sampling and outlier chain correction to accelerate the convergence towards the stationary distribution.

In the current work we run several parallel simulations for each MCMC method. In this approach, we need to determine when it is safe to stop the MCMC simulations for a reliable characterization of the permeability field. There are several convergence diagnostics available for this purpose and those diagnostics fall into two categories: the first category of diagnostics entirely depends on the output values of the MCMC simulation and those in the second category use not only the output values but also the information on the target distribution. In the first category, Brooks and Gelman [4] proposed a convergence diagnostic that uses the Multivariate Potential Scale Reduction Factor (MPSRF) to decide when to terminate MCMC simulations. Very recently in [2] we proposed a stopping criterion using a statistical analysis for single-phase flow prediction. In the analysis we considered ensembles of fractional flow curves to decide when it is safe to stop MCMC simulations for a reliable characterization and prediction. In this work we use the criterion in [2] for both MCMC methods and compare them for the characterization of the permeability field of the aquifer. Our results show that the two-stage MCMC provides a good estimate for the average of the quantity of interest (fractional flow curves), but the DREAM MCMC method reveals that the posterior distribution is not well characterized by the two-stage method.

We organize the present study as follows. We discuss the computational modeling of the single-phase flow of the aquifer in Sect. 2. The reduction of the parameter space dimension by the KLE is discussed in Sect. 3. Section 4 contains the statistical framework using two-stage and DREAM MCMC methods for the characterization of the permeability field in the aquifer. The MPSRF, a frequently used criterion to estimate convergence of the MCMC methods, is presented in Sect. 5. Results obtained from our numerical experiments are discussed in Sect. 6. Concluding remarks appear in Sect. 7.

We consider a unit square-shaped subsurface aquifer Ω with a heterogeneous permeability field. The aquifer contains two monitoring wells: one of which is located at the top right corner of the domain (corner well) and the other well (center well) is placed at the middle of the right boundary. An accidental contamination occurs at the spill site and the contaminated water flows naturally through the aquifer. The spill well is positioned at the bottom left corner through which the contaminated (or tracer) water is discharged (Fig. 1) . In this singlephase flow model we assume a relatively low concentrations of the contaminant, thus it does not affect the velocity field.

Let us denote the Darcy velocity and the pressure of the fluid by v(x, t), and p(x, t), respectively, where x ∈ Ω is the location at time t. We also denote the absolute permeability of the rock, porosity of the rock, fluid viscosity, and contaminant concentration in the fluid by k(x), φ(x), μ, and ρ(x, t), respectively. We consider that the pore space of the aquifer is filled by water. Applying Darcy's law and mass conservation, the governing equations describing the single-phase flow can be written as the following [5] :

In this study the porosity of the rock is considered a constant throughout the domain with φ(x) = 0.2. We assume no-flow boundary conditions and we take ρ(x, t = 0) = 0. The system of partial differential equations in (1) does not contain any source or sink because all three wells are modeled through appropriate boundary conditions. The coupled system consists of an elliptic problem and a hyperbolic problem. After applying an operator splitting technique, we solve each problem separately by an appropriate numerical scheme [13, 18] . The permeability field is characterized by using measured data in the form of fractional flow curves, F (t) which are defined as

where ∂Ω out represents the well outflow boundary, and v n (x, t) is the component of the velocity field normal to the well boundary. The dimensionless time is denoted by t, which is measured in Pore Volume Injected (PVI) and is computed using the following integral:

where V p denotes the total pore volume of the reservoir and T represents the total time the contaminated water entered through the spill well.

If we consider the aquifer in a 128 × 128 computational domain, we have 16, 384 elements. In our Bayesian framework the numerical simulator uses the permeability value in each element. Therefore, we need to start with a random permeability value in each element and change one or more of those values in each MCMC iteration. Thus, the dimension of the parameter space is 16, 384 and it presents a far-fetched framework for the characterization of the permeability field. In order to reduce the number of uncertainty parameters, we therefore use the KLE. Below we reduce the parameter space from 16, 384 to d, which is 20 for our study, using the KLE. The KLE has been explained in [10] [11] [12] . A very short discussion on the KLE is being presented below for the sake of completeness.

is a Gaussian field and the covariance function Cov(Y k (x 1 ), Y k (x 2 )) is given by the following formula:

where L x and L y are the correlation lengths in x− and y−direction, respectively, and

can be expressed as the following:

where Y k i are random coefficients. On the other hand, since L 2 is a complete space, thus ϕ i (x) is an eigenfunction satisfying

and the corresponding eigenvalue (5) can be expressed as the following:

If the eigenvalues decrease, a truncated KLE can be written as (8) decay quickly, the Y k d will be a good approximation of Y k .

We set L x = L y = 0.2 and σ 2 Y = 4 in (4). Figure 2 shows that the eigenvalues decay fast for these values. In this paper we thus consider the first twenty eigenvalues in the KLE to model the permeability field.

In this subsection we discuss how to characterize the permeability field using the fractional flow curves at the monitoring wells of the aquifer. Let us denote the fractional flow data by F m and the corresponding permeability field by ψ ψ ψ. Using the Bayes' theorem we can write the posterior probability

where P (ψ ψ ψ) denotes the prior distribution and the normalizing constant is ignored due to the iterative search in the MCMC algorithm. The ψ ψ ψ is generated through the KLE, for which the vector θ θ θ is used as input in the expansion.

In the remainder of the discussion we use ψ ψ ψ = KLE[θ θ θ]. Moreover, we consider a Gaussian distribution as in [9] for the likelihood function, i.e.,

where the simulated fractional flow data F ψ ψ ψ is obtained by the numerical solution from the GPU simulator for each permeability distribution ψ ψ ψ in the MCMC algorithm. We denote the covariance matrix by Σ, which is defined as Σ = I I I/2σ 2

where I I I and σ 2 F are the identity matrix and the precision parameter, respectively. We sample data from the posterior by using the Metropolis-Hasting (MH) MCMC and create a Markov chain, which has the posterior distribution as target distribution. We consider an instrumental distribution q(ψ ψ ψ p |ψ ψ ψ), where ψ ψ ψ denotes the previously accepted proposal, in order to propose ψ ψ ψ p = KLE[θ θ θ p ] at every iteration. We use the following acceptance probability in MH MCMC and the probability value is computed by solving the forward problem for a given permeability distribution on the numerical simulator:

We now describe the two-stage and DREAM MCMCs that use the MH algorithm.

Here we present the two-stage MCMC method. The method has been widely used for porous media applications [8, 16] . The two-stage MCMC consists of a screening procedure, which relies on a coarse-scale model approximating the governing equations (1). The coarse-scale discretization is done in a similar way as in the fine-scale discretization. The main idea lies on a rigorous projection of k on the coarse-scale that is obtained from the fine-scale resolution. For this reason, an upscaling method is used so that the effective permeability values on the coarse-scale yield a similar average response as that of the underlying fine-scale problem locally [11] . We then run the numerical simulator on the coarse-scale model and get the numerical solution F c . Now we can compute the coarse-scale and fine-scale acceptance probabilities α c (ψ ψ ψ, ψ ψ ψ p ) = min 1, q(ψ ψ ψ|ψ ψ ψ p )P c (ψ ψ ψ p |F m ) q(ψ ψ ψ p |ψ ψ ψ)P c (ψ ψ ψ|F m ) and α f (ψ ψ ψ, ψ ψ ψ p ) = min 1,

In the two-stage MCMC the following Random Walk Sampler (RWS) is used:

where θ and θ p represent the previously accepted proposal and the current proposal, respectively. The symbol stands for a N (0, 1)-random variable and β (= 0.75) is a tuning parameter [1] . The two-stage MCMC algorithm is presented in Algorithm 1. The convergence diagnostic to break the "for loop" in the algorithm is described in Sect. 5.

The DREAM, which is an extension of DE MCMC [3] , runs multiple MCMCs simultaneously for a thorough exploration of the posterior, and has an in-built mechanism to adapt the scale and orientation of the proposal distribution during the evolution to the posterior distribution [22] . We describe the DREAM for our Compute the upscaled permeability on the coarse-scale using ψ ψ ψp.

Solve the forward problem on the coarse-scale to get Fc.

Compute the coarse-scale acceptance probability αc(ψ ψ ψ, ψ ψ ψp).

if ψ ψ ψp is accepted then 8:

Use ψ ψ ψp in the fine-scale simulation to get F f . 9:

Compute the fine-scale acceptance probability α f (ψ ψ ψ, ψ ψ ψp). 10:

if ψ ψ ψp is accepted then ψ ψ ψ = ψ ψ ψp. 

where δ and d denote the number of pairs that are used to generate the proposal and the number of parameters that are updated jointly in each iteration, respectively. Two randomly chosen chains are denoted by r 1 

Use ψ ψ ψp to get F f (on the fine-scale). 6:

Compute the acceptance probability α(ψ ψ ψ, ψ ψ ψp) using (11). 7:

if ψ ψ ψp is accepted then ψ ψ ψ = ψ ψ ψp. 8: end if 9: end for 10: end for We simultaneously run the DREAM MCMC with m = 11 parallel chains. In (14) we set δ = 5, b = 0.1 and b * = 10 −6 . The jump rate is given by

, where the constant β 0 = 1 16 . The user should select the value of β 0 in such a way the MCMC method has an acceptance rate between 15-35%. In the present study we set d = 5, i.e., we update five parameters at every iteration. In addition to setting those parameters, we also set γ = 1.0 at every tenth iteration to encourage a jump between two disconnected posterior modes.

The DREAM MCMC requires a parallel simulation of multiple chains simultaneously. However, the two-stage MCMC does not require the same. Due to the computational burden in repeatedly solving the numerical simulator, we still run the parallel simulation of several two-stage MCMCs. Now we need to investigate the convergence of multiple chains in each method for a reliable characterization of the permeability field of the aquifer. In this section we discuss convergence diagnostics for that purpose. The Potential Scale Reduction Factor (PSRF) and its multivariate extension MPSRF are used to measure the convergence of multiple MCMCs [4] .

We set the number of parameters as d = 20 in θ as discussed in Sect. 3. Let us consider m chains and n posterior draws of θ in each chain. θ t i refers to the vector θ at iteration t in the ith chain of multiple MCMCs. Then the posterior variance-covariance matrix in higher dimensions is computed by

The within-covariance-matrix W is given by

and the between-chain-covariance-matrix B is computed by

whereθ i. denotes the mean of θs within the chain andθ .. represents the mean of θs in all the chains. The PSRF is defined as follows:

The PSRF values close to one indicate that the samples in multiple chains are being generated from the same limiting distribution and thus confirm the convergence. The MPSRF is computed as follows [17] :

where λ 1 is the greatest eigenvalue of W −1 B/n. If the MPSRF approaches one for a reasonably large n, the convergence of the multiple chains is ensured. Moreover, a relationship between the PSRFs and MPSRF is given by [4] max of PSRFs ≤ MPSRF. (20) 6 Numerical Results

In this subsection we present a simulation study of the characterization of the heterogeneous permeability field of the aquifer. See Fig. 1 . The aquifer is not contaminated initially. The contaminated water flows through a spill well into the aquifer at a rate of one pore-volume every five years. The synthetic reference permeability field is constructed on a fine-grid of size 128 × 128. The fractional flow curves for this reference field at the monitoring wells are shown in Fig. 5 . These curves are generated by running the numerical simulator on the reference field until t = 1.0 PVI. The time evolution of the contaminant flow on the reference permeability field is shown in Fig. 1 . We use the two-stage and DREAM MCMC methods to generate samples from the posterior. A coarse-grid of size 32 × 32 is chosen for the two-stage algorithm, which runs four times faster than the fine-grid simulation, but still manages to capture the general trend of the flow. See Fig. 3 . 

This subsection contains the MPSRF and PSRF analysis for both two-stage and DREAM MCMCs. We run twelve and eleven chains for two-stage and DREAM MCMCs, respectively. Using an equal number of total proposals in all chains, we then compute the maximum of PSRFs and the MPSRF against the number of iterations. Figure 4 shows that the MPSRF for each method is the upper bound of the maximum of the PSRFs, which is consistent with the inequality shown in (20) . Moreover, it is also observed that both MPSRF and the maximum of PSRFs for the DREAM MCMC have a faster downward trend than the two-stage MCMC method. Thus, we can conclude that the DREAM MCMC samples from the posterior faster than the two-stage MCMC and converges faster towards the stationary distribution. However, Fig. 4 demonstrates that both MCMC methods need a large number of iterations to achieve a complete convergence. To achieve a complete convergence, the curves should approach the numerical value of one. Next we focus on a statistical analysis of two ensembles of accepted fractional flow curves. In the analysis we consider the variances of the ensembles of fractional flow curves as well as the posterior distributions.

We compute the average production curves using 24000 accepted proposals and compare those curves with the reference fractional flow curves. See Fig. 5 . From the comparison we can say that both two-stage and DREAM MCMCs produce very similar fractional flow curves. Table 1 shows the acceptance rates for both MCMCs, where σ 2 F and σ 2 C are precision parameters in (10) for coarse-and finescale simulations, respectively. Note that the acceptance rates for both MCMC methods are the same, however, the convergence rate in the DREAM MCMC is considerably higher than that in the two-stage MCMC (see Fig. 4) .

We now construct two ensembles by taking the same number of fractional flow curves in each MCMC method: The first ensemble contains 24000 samples and the second one has 36000 samples. Figure 6 displays the variances of the fractional flow curves for the center and corner wells of the aquifer. The variance curves differ between not only both ensembles but also both MCMC methods. Furthermore, we sketch the posterior densities in Fig. 7 for the same ensembles that we considered for the variance analysis. The normalized frequencies reveal that the posterior densities are also different for both MCMC methods. 

Using a GPU-based single-phase flow simulator we have compared two frequently used MCMCs, the two-stage algorithm based on a random walk sampler and the DREAM. We have confirmed that the DREAM converges much faster than the two-stage MCMC by comparing the corresponding PRSF and MPRSF curves. Moreover, a careful statistical analysis of ensembles of accepted fractional flow curves, produced by the two MCMCs, reveals that such ensembles share essentially the same average behavior. However, significant differences have been observed in the time-dependent variance curves as well as in the posterior distributions for those ensembles. This provides an indication that, for the purpose of making Monte Carlo predictive simulations, one might observe considerable differences in the results. One could combine both MCMC methods in a two-stage version of DREAM to achieve good convergence along with reduced computational cost in line with the work of [14] . As a future work we intend to combine both methods for our problem and study the convergence.

A Bayesian framework for the validation of models for subsurface flows: synthetic experiments

Contaminant transport forecasting in the subsurface using a Bayesian framework

A Markov chain Monte Carlo version of the genetic algorithm differential evolution: easy Bayesian computing for real parameter spaces

General methods for monitoring convergence of iterative simulations

Computational Methods for Multiphase Flows in Porous Media

Markov chain Monte Carlo using an approximation

Bayesian calibration of a large-scale geothermal reservoir model by a new adaptive delayed acceptance metropolis hastings algorithm

An efficient twostage Markov chain Monte Carlo method for dynamic data integration

Preconditioning Markov chain Monte Carlo simulations using coarse-scale models

Multiple Markov chains Monte Carlo approach for flow forecasting in porous media

A multi-stage Bayesian prediction framework for subsurface flows

A prefetching technique for prediction of porous media flows

Multi-physics Markov chain Monte Carlo methods for subsurface flows

Efficient posterior exploration of a high-dimensional groundwater model from two-stage Markov chain Monte Carlo simulation and polynomial chaos expansion

Probability Theory

An efficient two-stage sampling method for uncertainty quantification in history matching geological models

Markov switching models: an application to roadway safety

A semi-discrete central scheme for the approximation of two-phase flows in three space dimensions

PyDREAM: high-dimensional parameter inference for biological models in python

Dream (D) : an adaptive Markov Chain Monte Carlo simulation algorithm to solve discrete, noncontinuous, and combinatorial posterior parameter estimation problems

Markov chain Monte Carlo simulation using the DREAM software package: theory, concepts, and MATLAB implementation

Treatment of input uncertainty in hydrologic modeling: doing hydrology backward with Markov chain Monte Carlo simulation

A. Rahunanthan's research was supported by National Science Foundation under Grant No. HRD-1600818. The computer simulations were performed on the NSF-funded GPU computing cluster housed in the Department of Mathematics and Computer Science at Central State University.