A multiobjective Markov chain Monte Carlo approach for history matching and uncertainty quantification

Journal of Petroleum Science and Engineering 166 (2018) 759–777

Contents lists available at ScienceDirect

Journal of Petroleum Science and Engineering journal homepage: www.elsevier.com/locate/petrol

A multiobjective Markov chain Monte Carlo approach for history matching and uncertainty quantification Feyi Olalotiti-Lawal *, Akhil Datta-Gupta Texas A&M University, United States

A R T I C L E I N F O

A B S T R A C T

Keywords: Subsurface model calibration Uncertainty quantification Markov chain Monte Carlo Low-rank parameter representation Multiobjective optimization

We present a probabilistic approach for integrating multiple data types into subsurface flow models. Our approach is based on a Bayesian framework whereby we exhaustively sample the multi-dimensional posterior distribution to define a Pareto front which represents the trade-off between multiple objectives during history matching. These objectives can be matching of water-cut, GOR, BHP and time-lapse seismic data. For field applications, these objectives do not necessarily move in tandem because of measurement errors and also interpretative nature of the seismic data. Our proposed method is built on a Differential Evolution Markov Chain Monte Carlo (DEMC) algorithm in which multiple Markov Chains are run in parallel. First, a dominance relationship is established amongst multiple models. This is followed by construction of the posterior distribution based on a hypervolume measure. A unique aspect of our method is in the parameter proposal generation which is based on a random walk on two arbitrarily selected chains. This promotes effective mixing of the chains resulting in improved convergence. We illustrate the algorithm using a nine-spot waterflood model whereby we use water-cut and bottomhole flowing pressure data to calibrate the permeability field. The permeability field is re-parameterized using a previously proposed Grid Connectivity Transform (GCT) which is a model order reduction technique defined based only on the decomposition of the grid Laplacian. The compression power of GCT allows us to reconstruct the permeability field with few parameters, thus significantly improving the computational efficiency of the McMC approach. Next, we applied the method to the Brugge benchmark case involving 10 water injectors and 20 producers. For both cases, the algorithm provides an ensemble of models all constrained to the history data and defines a probabilistic Pareto front in the objective space. Several experimental runs were conducted to compare the effectiveness of the algorithm with Non-Dominated Sorting Genetic Algorithms (NSGA-II). Higher hypervolume was constantly measured using our algorithm which indicates that more optimal solutions were sampled. Our method provides a novel approach for subsurface model calibration and uncertainty quantification using McMC in which the communication between parallel Markov chains enhances adequate mixing. This significantly improves the convergence without loss in sampling quality.

1. Introduction In all fields of modeling, it is essential to calibrate (and re-calibrate) models to reproduce available data for reliable predictions and prudent decision making. Hydrocarbon reservoir models are typically calibrated by constraining them to different sources of data including conceptual depositional models, time lapse seismic data, well logs and production history which are acquired at varying scales and resolution. Although there are multiple sources of data, they are often sparse and provide little information of the subsurface, resulting in a wide range of uncertainties

especially in the flow and storage capacity. Integrating multiple data sets into reservoir models helps to reduce these uncertainties as well as improve the predictive capability of the updated models (Haldorsen, 1986; Omre and Tjelmeland, 1997). The process of integrating these data sets into reservoir models is conventionally referred to as history matching. It has evolved from the less reliable tedious and time-consuming manual method (Saleri and Toronyi, 1988) to assisted approaches which often involve the application of re-parameterization and sometimes, model order reduction techniques, together with optimization algorithms/workflows to obtain

* Corresponding author. E-mail address: [email protected] (F. Olalotiti-Lawal). https://doi.org/10.1016/j.petrol.2018.03.062 Received 15 August 2017; Received in revised form 10 January 2018; Accepted 14 March 2018 Available online 16 March 2018 0920-4105/© 2018 Elsevier B.V. All rights reserved.

F. Olalotiti-Lawal, A. Datta-Gupta

Journal of Petroleum Science and Engineering 166 (2018) 759–777

the Weighted Sum Method as the target Pareto front gets more concave, an intrinsic shortcoming of gradient-based methods are their relatively poor parameter space exploration capability. In other words, global optimum solutions are rarely obtained when starting (guess) solutions are far from the target. Schulze-Riegert et al. (2007) applied the Strength Pareto Evolutionary Algorithm (SPEA) in resolving the trade-off between a priori established conflicting objectives in a synthetic reservoir history matching problem. Hajizadeh et al. (2011), with a synthetic reservoir model calibration problem, demonstrated the superior robustness of a Differential-Evolution-based multi-objective optimization algorithm over the weighted sum method. Likewise, Park et al. (2015) applied the Non-Dominated Sorting Genetic Algorithm (NSGA-II) (Deb et al., 2002) in integrating multiple sources of data into reservoir models using an illustrative synthetic case and an application to the Brugge benchmark reservoir model. A multiobjective Particle Swarm Optimization algorithm, another heuristic-based algorithm, was applied to perform a history match of a field case (Christie et al., 2013). Most of these algorithms are designed to obtain the Pareto set of solutions in a deterministic sense, and the set of solutions only provide trade-off among objectives deterministically. Multiobjective optimization, together with uncertainty quantification algorithms are, therefore, required in providing Pareto-optimal solutions while quantifying uncertainties. In this approach, solutions in the vicinity of the optimal trade-off solutions are sampled as well so that each tradeoff solution comes with a range of uncertainty (Zhang et al., 2008). This becomes particularly useful for certain reservoir model calibration problems in which the data sets exhibit varying levels of uncertainties. Some applications of probabilistic multiobjective algorithms have been published in the water resources and hydrology literature (Singh and Minsker, 2008). proposed and applied a Probabilistic Multiobjective Genetic Algorithm (PMOGA) in the study of simultaneous effects of uncertainty on multiple objectives for groundwater remediation design. On another account, W€ ohling and Vrugt (2008) proposed and applied a technique referred to as Bayesian Model Averaging (BMA) for groundwater model calibration in which the sampled probability distribution function (PDF) is obtained by a weighted sum of individual PDFs of each objective function. We present a multiobjective Markov chain Monte Carlo approach based on a Differential Evolution algorithm for probabilistic history matching of subsurface flow models. Similar to Li (2012), our approach applies the Differential Evolution Markov Chain Monte Carlo (DEMC) algorithm (ter Braak, 2006; ter Braak and Vrugt, 2008) in sampling a probability distribution function defined by a hypervolume, which represents the dominance relationships among solutions in the objective space (Zitzler and Thiele, 1999). The hypervolume is an indicator known to be fully sensitive to Pareto dominance (Bader and Zitzler, 2011) and also avoids a weighted sum approach in probability distribution calculations as in the BMA algorithm discussed earlier. The DEMC involves running multiple Markov chains in parallel such that information is exchanged between the chains during the exploration phase. This improves mixing and enhances the overall performance compared to the traditional McMC algorithm (ter Braak, 2006). For our case studies, we consider bi-objective optimization problems in which we seek to obtain multiple reservoir models constrained to two sets of data at the trade-off region of the probabilistic Pareto front. We must unequivocally state here that not all reservoir model calibration problems require a multiobjective optimization routine, since it is not always the case that all data sets are conflicting. For all of our applications, we re-parameterize the permeability field using a model order reduction and image compression technique called the Grid Connectivity-based Transform (GCT) (Bhark et al., 2011). By this, only few coefficients are required for generating parameter proposals during the sampling process, which significantly reduces the dimensionality of the problem and grossly improve the sampling efficiency and overall convergence of the algorithm. The paper is organized as follows. First we describe the multiobjective optimization problem, highlighting the contrast between the

updated subsurface models which preserve geologic realism. Assisted history matching workflows have been developed and applied in different forms including gradient-based (Dong and Oliver, 2008), sensitivity-based (Cheng et al., 2007; Tanaka et al., 2015), evolutionary (Yin et al., 2011), and ensemble-based methods (Aanonsen et al., 2009), with different levels of computational cost and ease/difficulty of implementation. Gradient/sensitivity based methods possess certain desirable Newton/quasi-Newton characteristics which generally promote faster convergence to optimal solutions. However they typically only generate local solutions which means that, depending on the degree of non-linearity of the objective function, they are most effective when initial guesses are close to the optimal solutions in the parameter space. Population-based or stochastic methods being gradient-free on the other hand, are more expensive but perform better at locating global solutions. Combining the strengths of both methods in a hierarchical fashion can lead to improved assisted history matching techniques (Yin et al., 2011). Advancements in computing technology in recent years has spurred a growing interest in the application of population-based and stochastic algorithms to reservoir history matching and uncertainty quantification problems. Romero et al. (2000), by building structured chromosomes to handle multiple data sets, applied a modified Genetic Algorithm (GA) to integrating multiple data into a synthetic reservoir model. In another work, Castellini et al. (2005) applied the GA using an iterative generation of proxy models to obtain multiple reservoir models constrained to history data in a West Africa field case study. The Particle Swarms Optimization (PSO) algorithm has also been applied to history matching and uncertainty quantification (Mohamed et al., 2010) and geophysical inversion problems (Fern
andez-Martínez et al., 2008). Hajizadeh et al. (2010) demonstrated the applicability of the Differential Evolution (DE) algorithm to reservoir history matching problems. In their work, they applied variants of the algorithm to both synthetic and field cases, studying the trade-offs between parameter space exploration and rate of convergence. Beyond population-based heuristics, the Markov chain Monte Carlo (McMC) algorithm samples directly from the posterior distribution constructed from a combination of the prior and likelihood (hard data) distributions. McMC provides the true representation of model uncertainty in a Bayesian probabilistic framework and has been applied to subsurface model calibration and uncertainty quantification problems (Maucec et al., 2007; Oliver et al., 1997). Two-stage McMC methods have been introduced and applied (Ma et al., 2008; Olalotiti-Lawal, 2013) to reduce the number of rejected simulation runs, thereby improving the overall efficiency of the sampling process. Typically, multiple data sets such as production, well test and seismic data have been integrated by minimizing an objective function defined by a simple weighted sum of respective misfits in each data set. While intuitive, this sometimes causes poor exploration of the parameter space, resulting in premature convergence (Deb et al., 2002). Furthermore, results obtained can be misleading from the decision making stand point when the objective functions tend to trade-off. In such cases, multiobjective optimization algorithms are applied to obtain multiple trade-off solutions. Each of these solutions is optimal in at least one objective function and cannot be optimized further without degrading the optimality of other objective functions. This set of solutions is often referred to as the Pareto-optimal set (Zitzler and Thiele, 1999). Multiobjective optimization problems have likewise been approached with gradient-based as well as population-based methods. Izui et al. (2015) proposed an adaptive weighting of solutions in a linear programming approach to establish the relative dominance in the objective space. In an oil reservoir waterflood optimization problem, Liu and Reynolds (2016) investigated the application of the Normal Boundary Intersection (NBI) method, previously proposed by Das and Dennis (1997). In this approach, the solution is initialized with a ‘utopia’ line which represents a locus of linear combinations of the objective extremes (initial guess of trade-off solutions), and at every iterations, these points are shifted following normal directions to this line. While these gradient-based methods elegantly improves on the known weakness of 760

F. Olalotiti-Lawal, A. Datta-Gupta

Journal of Petroleum Science and Engineering 166 (2018) 759–777

Fig. 1. Illustration of the concept of dominance.

deterministic and probabilistic perspectives. Next we discuss the solution approach from the construction of probability distribution function construction to details of the sampling algorithm. Then we illustrate the concept and the utility of the algorithm using a simple two dimensional synthetic reservoir model calibration problem. Application to the Brugge benchmark case (Peters et al., 2010) follows, by which we demonstrate the practical feasibility of the approach. Finally we present the discussion of results and conclusion. 2. Problem description In this work, we attempt to solve a continuous multiobjective optimization problem defined as follows: minfðxÞ ¼ ½f1 ðxÞ; f2 ðxÞ; :::; fm ðxÞT

Fig. 2. Illustration of the Pareto optimality concept in deterministic and probabilistic perspectives.

(1) sum method typically suffers certain deficiencies. Das and Dennis (1997) clearly showed the limitation of the weighted sum method with non-convex Pareto fronts, for which case the method typically fails to cover the needed Pareto solution spectrum entirely. Furthermore, issues arise with scaling and weights assignment to objectives especially with history matching problems where we seek to integrate different data sets at different resolution and scales (Hajizadeh et al., 2011). The multiobjective optimization problem formulation will be necessary in such cases to circumvent the apparent limitations of the weighted sum method. To explain some basic concepts, a bi-objective optimization problem is depicted in Fig. 2 showing a Pareto-optimal set in the parameter space (Fig. 2(a) and (c)) and in the objective space (Fig. 2(b) and (d)). In these pictures, F is a vector-valued function that maps some parameter vector X to the objective space. The majority of algorithms have been designed to target the Pareto set deterministically as shown in Fig. 2(a) and (b), while a few approach the problem in a probabilistic sense. The basic idea in a typical probabilistic approach is to employ statistical methodologies for characterizing sampled solutions with an end goal of obtaining a population of decision vectors which are uniformly scattered around the Pareto-optimal set (Zhang et al., 2008). This is illustrated in Fig. 2(c) and (d). Our approach employs a Bayesian framework using a parallel Markov chain Monte Carlo (McMC) algorithm for sampling the posterior distribution defined by the dominance relationship among solutions in the objective space.

subject to x 2 X⊂ℝ

n

where x denotes the decision vector in the n-dimensional decision variable space X while fðxÞ denotes the m-dimensional vector of objective functions. For subsurface model calibration purposes, the decision vector may be represented by the dynamic model parameters such as permeability, whereas the misfit in model responses such as bottomhole pressure and production water cut make up the objective function vector. A common strategy for approaching multiobjective optimization problems involves obtaining a set of trade-off solutions, otherwise known as the Pareto-optimal set, which represents decision vectors whose corresponding objective functions cannot be improved in any direction without degrading the other (Zitzler and Thiele, 1999). Usually Pareto based optimization algorithms require the evaluation of the dominance relationship among a set of solutions fxI 2 XjI ¼ 1; 2; :::; Ng, where N is the number of solutions in X. Typically, according to Deb et al. (2002), a solution xI is said to dominate another solution xJ , that is xI  xJ if: 8i 2 f1; 2; :::; mg : fi ðxI Þ  fi ðxJ Þ ^ 9j 2 f1; 2; :::; mg : fj ðxI Þ < fj ðxJ Þ

(2)

An illustration of the dominance concept is shown in Fig. 1 where solutions have been characterized relative to a specific solution O in the objective space. Solutions belonging to Quadrant C have higher objective function values along both axes, and therefore are worse than solution O. Solutions in Quadrant A however, have lower objective function values along both axes, making them better solutions compared to solution O. Therefore, any solution located in the Quadrant A is said to dominate solution O and conversely, any solution in Quadrant C is dominated by solution O. Quadrants B and D on the other hand, contain solutions that are comparable to solution O since there are relative improvements in a single objective function. A common critique is to question the necessity of a multiobjective optimization problem formulation when simple weighted sum approach might just suffice. A multiobjective optimization problem formulation will absolutely not be required for cases where the objectives are known to sufficiently correlate. However when there are indications of conflicts among different objectives, a scenario which is often observed with integrating multiple sources of data in reservoir models, the weighted

3. Solution approach In solving the multiobjective optimization problem, we seek a set of solutions in the vicinity of the Pareto-optimal set such that uncertainties (mainly due to errors or noise in the data to be integrated) are captured for every trade-off solution along the Pareto-optimal set. For this, we construct a probability distribution function (PDF) based on the dominance relationship among solutions in the objective space and then sample using a McMC algorithm. In this paper we adopt the Differential Evolution Markov chain Monte Carlo (DEMC) algorithm (ter Braak, 2006) for sampling the complex PDF. Our approach follows a previously proposed multiobjective Markov chain Monte Carlo algorithm (Li, 2012), however our proposal generation protocol satisfies detailed balance 761

F. Olalotiti-Lawal, A. Datta-Gupta

Journal of Petroleum Science and Engineering 166 (2018) 759–777

3.2. Proposal generation and distribution sampling

(Robert and Casella, 2013), providing an assurance of converging to the target stationary distribution. We also improved the solution diversity by incorporating a crowding distance function in the PDF construction. The details of the PDF construction and the sampling algorithm are as follows.

Due to the complex nature of the PDF described in the foregoing subsection and documented difficulties of traditional McMC algorithm with sampling such distributions, we adopt a relatively robust and efficient sampler known as the Differential Markov chain Monte Carlo algorithm. As the name suggests, it is a combination of two independent algorithms: Differential Evolution (DE) (Storn and Price, 1997) and the Markov chain Monte Carlo (McMC) sampling algorithm (Hastings, 1970). Differential Evolution (DE), a fairly recent population-based optimization algorithm which has rapidly gained popularity in the technical computing world, has found applications in petroleum engineering problems including assisted history matching (Mirzabozorg et al., 2013) and seismic waveform inversion (Wang et al., 2011). Hajizadeh et al. (2011) applied the DE algorithm in approaching reservoir model calibration tasks posed as multiobjective optimization problems, but this was without regard to possible existence of noise in the data sets at varying scales. For true probabilistic data integration and uncertainty quantification, the Markov chain Monte Carlo algorithm (Robert and Casella, 2013), as a versatile probability distribution sampling algorithm, fits the purpose for rigorous sampling from the posterior distribution of optimal solutions (Xie et al., 2011). It generally requires large amounts of functional evaluations, and is thus computationally expensive for subsurface model calibration problems. However, recent exponential improvements in computing power, coupled with development of hybrid and more efficient McMC algorithms (Ma et al., 2008; Olalotiti-Lawal, 2013), makes the method attractive for the required high fidelity representation of all possible scenarios of subsurface characterization. The Differential Evolution Markov Chain Monte Carlo (DEMC), initially proposed by (ter Braak, 2006) is a population-based McMC algorithm in which multiple chains are run in parallel. As the chains explore the specified parameter domain, information is exchanged among the chains using the DE-type proposal generation protocol. In contrast to running independent Markov chains in parallel for convergence diagnostic purposes (Maucec et al., 2007), DEMC significantly enhances mixing and thus, improves the performance of the sampling. A simple study was conducted to compare the relative sampling efficiencies and robustness between the DEMC and a traditional McMC algorithm. Details and results of the experiment are provided in Appendix A. The results show clearly that DEMC is superior at consistently sampling the target distribution compared to the traditional McMC algorithm. This is a direct implication of the improved mixing achievable with DEMC which helps to avoid some of the deficiencies of a traditional McMC algorithm when sampling for multimodal distributions. In the DEMC algorithm, a proposal xp at the t th state of a Markov chain of the I th Markov chain is given by the sum of its current state with the scaled difference between the current states of two randomly selected chains (R1 and R2), mathematically expressed as:

3.1. Probability distribution function We are interested in sampling the likelihood distribution of possible models defined with a proportionality relation in the parameter space as:

π ðxÞ∝expðγðxÞ=TÞ; T > 0

(3)

where T is the simulated temperature, a parameter that controls the acceptance rate of the proposals. This form of distribution is frequently encountered in the field of statistical physics (Marinari and Parisi, 1992) where it is of interest to sample potential energy distribution of thermodynamic systems, otherwise known as the Gibbs distribution. Such distributions, especially when they are multimodal, are typically tempered to improve mixing and overall performance of sampling (Geyer and Thompson, 1995). This will be discussed in more detail in subsequent sections. In our implementation, the exponent γðxÞ contains two kinds of information about the solution x: its dominance relationship with other solutions measured in the objective space ηðxÞ, and its diversity level Δ½fðxÞ. The dominance function value for the I th solution ηðxI Þ depends on whether it is dominated or non-dominated in the generation of solutions. If xI is non-dominated, that is fxJ jxJ  xI g ¼ ∅, then ηðxI Þ is the fraction of non-dominated solutions in the generation. If otherwise xI is dominated, ηðxI Þ is expressed as the sum of the pairwise hypervolumes (measured in the objective space) between xI and the set of solutions that dominate it, concisely expressed as:

ηðxÞ ¼ 1 þ

xX J xI

HYPðf I ; f J Þ

(4)

J6¼I

where f I ¼ fðxI Þ is the objective function vector (described for history matching problems in more detail later) corresponding to the I th solution in the parameter space for a certain generation. The pairwise hypervolume HYPðf I ; f J Þ is the size of objective space covered between objective vectors f I and f J (Zitzler and Thiele, 1999). To improve the diversity in the Pareto optimal set, we augmented the likelihood exponent with a measure of diversity Δ½fðxÞ, so that the exponent becomes: γðxÞ ¼ maxð0; ηðxÞ  Δ½fðxÞÞ

(5)

This measure of diversity is simply the minimum normalized distance between a solution i and the rest of the members of a generation, measured in the objective space. That is,   Δ½f I  ¼ minðf I  f J Þ: J6¼I

(6)

  xpI ¼ xtI þ κ xtR1  xtR2 þ e

Since we are interested in obtaining a set of all possible optimal solutions that span the Pareto front of the optimization problem, an effective sampling algorithm must avoid clustered solutions. The measure of diversity assigns higher weights to optimal solutions that are spaced apart in the objective space, thus enhancing the solution diversity. Using Eq. (6), we increase the fitness of solution I by the minimum Euclidean distance to any solution in the objective space. By this reasoning, a solution with a higher diversity measure will have higher chance of survival for the next generation over a solution of comparable dominance characteristics, but with lover diversity measure. Generally, favorable solutions must have γðxÞ close to zero. Note that small values of γðxÞ does not necessarily mean that all objective values to be minimized will all tend to zero, as this will automatically defeat the entire purpose of a multiobjective optimization problem formulation. Rather, Pareto solutions are identified with small values of γðxÞ.

(7)

where xtI and xpI respectively denote the current and proposed state for the I th Markov chain. Current states of randomly selected chains are denoted by xtR1 and xtR2 . Here, the vector e is essentially some Gaussian noise consistent with the problem dimensionality d, i.e. e is drawn from Nð0; bÞd , where b is orders of magnitude smaller than the parameter values. The scale factor κ takes a suggested optimal value of pffiffiffiffiffiffiffiffiffi 2:38= ð2dÞ (ter Braak, 2006). DEMC forms an ergodic Markov chain with unique stationary distribution (ter Braak and Vrugt, 2008). A criterion required for convergence of a Markov chain to a target stationary distribution is the detailed balance. This simply ensures equal forward and reverse transition probabilities at any given state in the Markov chain. Therefore, following Aanonsen et al. (2009), detailed balance requires that:

762

F. Olalotiti-Lawal, A. Datta-Gupta

Journal of Petroleum Science and Engineering 166 (2018) 759–777

  π ðxp Þ ρðxp jxÞ ¼ min 1; π ðxÞ

(10)

3.3. Probabilistic multiobjective optimization algorithm A simplified graphical flow chart is provided in Fig. 3 and the details of the steps involved in the algorithm are described as follows: 1. Initialize N chain at the 0th states in the parameter space. For our purpose, the Latin Hypercube sampler is applied for this purpose. 2 For all chains: a. At the t th state of a chain I, make a proposal using Eq. (7). b Evaluate posterior distribution function using Eq. (3). These calculations are carried out independently in each chain using the current temperatures of each chain. c. Calculate acceptance probability and evaluate acceptance of the proposal using Eq. (10). 3. If the desired number of generations is not reached, go back to step 2, otherwise stop. To keep track of the Pareto front, we define a global hypervolume which is the volume of the objective space enclosed between the prior state and the non-dominated solutions. Referring to our earlier description of dominance relationships in Fig. 1, the reader would recall that the non-dominated solutions only refer to the category of solutions that do not have any solution in their respective A-Quadrants. The global hypervolume is graphically shown in Fig. 4(a), in which the nondominated solutions are highlighted in red rings and the prior state is located at the top right vertex of the objective space. This is not to be confused with the pairwise hypervolumes discussed in Eq. (4) which are little volumes of the objective space between two solutions. The global hypervolume, scaled between 0 and 1 is efficiently computed using the Monte Carlo approach (Bader and Zitzler, 2011) to avoid the computational complexities associated with the exact computations using the hypervolume by slicing objectives approach (While et al., 2006). The higher the value of the global hypervolume, the more optimal the non-dominated set of sampled solutions are. The typical trend in this hypervolume measure, as shown in Fig. 4(b), is an initial rise (which translates to the burn-in period) followed by a stabilization (during which samples are collected from the posterior distribution).

Fig. 3. Flow Chart for Multiobjective Markov chain Monte Carlo Algorithm.

∫ π ðxÞqðxp jxÞρðxp ; xÞdx ¼ ∫ π ðxp Þqðxjxp Þρðx; xp Þdxp

(8)

where π ðxÞ represents the density of the parameter vector x, qðxp jxÞ is the probability of proposing a jump from state x to xp , and finally ρðxp ; xÞ is the probability of accepting that proposal. Assuming both forward and inverse transformations x→xp are differentiable, then the jump acceptance probability is calculated as: 

ρðxp jxÞ ¼ min 1;





p π ðxp Þqðxjxp Þ  ∂x   p π ðxÞqðx jxÞ ∂x 

(9)

It is clear from (7) that the Jacobian of transformation in (9) is identically unity since both x and xp exist in the same space. Note that the instrumental distributions qðxp ; xÞ and qðx; xp Þ are equivalent, since the two chains are selected with equal probabilities in both forward and reverse jumps, the acceptance probability function identically reduces to the Metropolis-type acceptance criterion (Metropolis and Ulam, 1949):

4. Applications We apply the Multiobjective Markov chain Monte Carlo algorithm to subsurface model calibration problems. Often times, different sets of data

Fig. 4. Global Hypervolume measure (a) Calculation and (b) Typical trend. 763

F. Olalotiti-Lawal, A. Datta-Gupta

Journal of Petroleum Science and Engineering 166 (2018) 759–777

Fig. 5. (a) Prior and (b) Reference Permeability felds; (c) Basis Functions.

appear to be conflicting and they cannot all be matched to the same level of accuracy. Common examples are time-lapse seismic and production data because of the different length scales of heterogeneity scanned by the data and the subjectivity involved in the interpretation of the seismic data. In such cases, a set of calibrated models examining the trade-offs in data misfits are needed for prudent decision making. We applied the algorithm first to a synthetic case and then the Brugge model. For each case the observed data are the producing water cut (WWCT) and the producing bottomhole pressure (BHP).To start with, we define our history matching objective function and our model reparameterization scheme.

where u denotes the permeability field of equal size as the model property field whereas v is the M-length vector of basis coefficients while Φ is a NxM matrix, each column representing a member of the set of orthonormal basis. In our application, instead of parameterizing the permeability field, we parameterize permeability multiplier field for the prior permeability. Thus the permeability multiplier field is initially set to unity for all cells and adjusted during history matching. The desired flow simulation property m field can be constructed by performing an element-wise multiplication of the prior model m0 with the multiplier field: m ¼ m0 ∘Φv

4.1. Objective function

For history matching purposes, we perturb the basis coefficients v in the multiobjective McMC algorithm to generate updated multiplier fields which are imposed on the prior model for a global update of the absolute permeability field. A unique advantage of the GCT parameterization scheme is that all M basis coefficients are uncorrelated. Furthermore, compared to other classical parameterization schemes such as the Principal Component Analysis, relatively fewer basis functions (and coefficients) are needed for sufficient resolution in the property field due to the high compression ratio (Bhark et al., 2011).

The objective functions for each parameter set are follows a standard least square form: f ¼

N wells X

ðdi;obs  gi ðmÞÞT C1 D ðdi;obs  gi ðmÞÞ

(11)

i¼1

where CD represents the data covariance, di;obs denotes the observed data for well i. The simulated model response, denoted by gi ðmÞ, is a function of the model parameter set which is restricted to the permeability field for all cases considered for our demonstration.

4.3. Synthetic model Here we present an application of the Multiobjective McMC algorithm to a two-dimensional, 21  21-cell synthetic reservoir. The three phase model has a nine-spot pattern with a water injector at the center and oil producers at the edges as shown in Fig. 5(a). Also, Fig. 5(a) and (b) plots the initial and reference logarithm of permeability fields respectively. The initial permeability field is homogeneous with a value of 10mD while the reference permeability field has a logarithm of permeability ranging from 0.0 to 2.0. The model porosity is however kept constant at 0.1 throughout to simplify the problem. For this problem, we utilize 10 leading order basis functions in constructing multiplier fields used for updating the absolute permeability field. Thus, we are able to achieve a significant model order reduction

4.2. Model reparameterization We parameterized the permeability field using the Grid Connectivity Transform (GCT) (Bhark et al., 2011) which is based on a general transform basis suitable for structured and unstructured grid geometry. An eigen-decomposition of the grid connectivity Laplacian results in a set of orthonormal basis Φ. By projecting the permeability field onto the orthonormal basis system, we obtain a linear transformation of spatial parameters between physical and transform domains: v ¼ ΦT u ⇔ u ¼ Φv

(13)

(12)

Fig. 6. (a) Overall acceptance ratio (b) A plot comparing hypervolumes obtained using Multiobjective McMC and NSGA-II. 764

F. Olalotiti-Lawal, A. Datta-Gupta

Journal of Petroleum Science and Engineering 166 (2018) 759–777

0:05. In this problem, we seek a set of trade-off solutions that are constrained to well bottomhole pressure (BHP) and water cut data (WWCT). The probabilistic Pareto-optimal set of solutions will then be characterized depending on their relative bias towards a certain data set. Assuming constant data covariance for well BHP and WWCT data respectively, we sampled the parameter space using 45 parallel Markov chains through 100 generations using the proposed algorithm. The algorithm was compared with the Non Dominated Sorting Genetic Algorithm (NSGA-II) (Deb et al., 2002), which was designed for multiobjective optimization problems but in the deterministic sense, and has been applied to model calibration and optimization problems both within and outside the petroleum engineering community. For consistency, the NSGA-II algorithm was run on the problem with a population size of 45 through 100 generations, with genetic operation parameters simlar to (Park et al., 2015). Acceptance ratios, computed as the fraction of accepted parameter proposals through all the sampling chains at a particular generation, are shown in Fig. 6(a). Acceptance rations generally fall within the range of normal values for an McMC algorithm (Geyer and Thompson, 1995). Hypervolumes computed for both algorithms are compared in Fig. 6(b). While an agreement in the trends (monotonic increase) can be noticed, the hypervolume obtained using multiobjective McMC algorithm rises faster than that from the NSGA-II algorithm. The implication is that optimal solutions are sampled much faster using the multiobjective McMC algorithm. This improvement results from the adoption of a single measure of dominance (hypervolume) on one hand, and effective mixing between the Markov chains that improves exploration of the parameter

Fig. 7. Pareto-optimal set of solutions from multiobjective McMC and NSGAII algorithms.

from 441 correlated parameters (cell permeability values) to 10 uncorrelated parameters (basis coefficients). Reducing problem dimensionality and eliminating redundancies is essential in any model calibration or optimization problem as it helps improve adequate exploration of the parameter space. For this problem, the temperature ladder (Earl and Deem, 2005; Vousden et al., 2016) was designed logarithmically with values ranging from 0 to 4 and each Markov chain assigned a unique temperature. Both well bottomhole pressure and water cut data were assumed uncorrelated in time, so that their respective covariance functions took the form of σ 2p I and σ pw I respectively, where σ p ¼ 100 and σ w ¼

Fig. 8. Optimal Trade-off sample responses for (a) Well Water Cut and (b) bottomhole pressure. 765

F. Olalotiti-Lawal, A. Datta-Gupta

Journal of Petroleum Science and Engineering 166 (2018) 759–777

Simulation responses for optimal trade-off solutions are shown in Fig. 8. A reasonable match to the observed data can be seen from these plots. Also, the optimal trade-off samples exhibit significant diversity for uncertainty quantification purposes. Simulation responses were also obtained for the set of samples obtained from the other two extremes in the trade-off: the best well BHP match and the best WWCT match. The simulation responses for these are shown in Figs. B1 and B2 of the Appendix. Comparing the matches for all category of ensembles, it is apparent that while the quality of the matches depend on the type of data bias (BHP or WWCT), the optimal trade-off solutions provides a compromise between the two extremes. This is graphically illustrated in Fig. 9 which displays 3 random samples from each category of the solution ensemble and Fig. 10 which compares the ensemble averages of logarithm of permeability distribution with the reference permeability. Looking at the solutions obtained from the best WWCT matches, it turns out that the low permeability region in the lower right corner of the model is not captured. Likewise the best BHP match only weakly captures the high permeability region of the top left corner of the model. Both of these subtle features are however rectified in the set of optimal trade-off solutions. Furthermore, looking at the histograms in Fig. 10, the updated permeability fields obtained from the best BHP match tend to be more conservative than that of the best WWCT match solutions. In other words on the average, more cells in the best BHP match solutions are likely to have lower permeability values compared to the best WWCT match. The optimal trade-off solutions, as shown in the histogram, once again tend to compromise between the two extremes. Aside parameter uncertainty quantification, we also attempted to validate the quality of our parameter estimation using the proposed algorithm. Fig. 11 shows the normalized distribution of the basis coefficients obtained from all the three categories of solution ensembles. In each of these plots we overlay the true basis coefficients that are obtained by projecting the reference permeability into the spectral domain defined by a set of orthonormal basis functions, as described in Eq. (12). Clearly, we obtain estimates with all three groups of solutions. However, the best BHP match solutions provide the best parameter estimation while the

Fig. 9. Comparing log permeability fields of samples from Best BHP match, Best WWCT match and Optimal Trade-off ensembles.

space on the other hand. Comparing the Pareto-optimal sets of solutions obtained from both algorithms likewise shows a good agreement as shown in Fig. 7. The probabilistic Pareto solutions were obtained by collecting samples from the Markov chains, assuming a burn-in period of 50 states along each chain. Both algorithms clearly show a significant trade-off between the well BHP and WWCT data sets. We applied the local least square regression method in selecting the optimal trade-off point in the Pareto solution sets. In this method, first the entire solution set is partitioned into bins along Pareto front. Then from the centroids of each partition, we fit the solutions towards the left and right independently and compute the sum of individual regression errors. This is done for every discrete partition along the Pareto front and the bin that gives the minimum sum of errors is taken as the best compromise solution set. This location typically agrees with (non-dominated) solution with the highest hypervolume (W€ ohling and Vrugt, 2008).

Fig. 10. Comparing reference log of permeability with the mean of log permeability of samples in Best BHP match, Best WWCT match and Optimal Tradeoff ensembles. 766

F. Olalotiti-Lawal, A. Datta-Gupta

Journal of Petroleum Science and Engineering 166 (2018) 759–777

Fig. 11. Basis coefficient distribution for (a) Best BHP match samples, (b) Best WWCT match samples and (c) Optimal Trade-off samples; compared with true values (green triangles). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 12. A comparison between the Base Case, the Large CDW and Small CDP cases, in terms of mean logPerm field and distribution as well as quality of parameter matches and estimates of parameter uncertainty.

best WWCT match solutions provide the worst. This underscores the importance of matching pressures in reservoir model calibration problems.

algorithm to a history matching problem with constant data variance. In this section, we investigate the performance of the algorithm in cases of varying values of noise through production history. This is considered closer to reality as, for instance, production water cut and Gas Oil Ratio (GOR) are typically known with greater level of confidence at breakthrough time than at later times. In our experiment, we ran the base case with σ 2w increasing exponentially from 5  104 to 2:5  103 over the 38

4.4. Problems with varying data variance The last section illustrates the application of the multiobjective McMC 767

F. Olalotiti-Lawal, A. Datta-Gupta

Journal of Petroleum Science and Engineering 166 (2018) 759–777

Fig. 13. A comparison between the Base Case, the Large CDW case in terms of water cut responses.

Fig. 14. A comparison between the Base Case, the Small CDP case in terms of bottomhole pressure responses.

forward simulation time steps, while σ 2P ¼ 10000 as in the previous example. Then two cases were created: First a case with σ 2P ¼ 2000 and default σ 2w ; and second, a case with default σ 2P while σ 2w increases exponentially from 5  104 to 5  103 . Note that the WBHP data noise was reduced in the first case, which practically translate to the fact that we have higher confidence in the WBHP data compared to the default case. Whereas, in the second case, we have increase the WWCT data noise, which can also be interpreted to be a less confidence in the data compared to the default case. For all runs, the temperature ladder design was kept the same as in the previous illustration, and likewise the number of Markov chains and number of generations kept the same at 45 and 100 respectively. The result of the parameter estimation for all three cases is shown in Fig. 12. Note that only optimal trade-off (compromise) solutions are shown for our comparison, and displayed permeability fields represent the mean of the ensemble of those solutions for all three cases. We see clearly that parameters (basis coefficients) are very well captured, although slight variations in the update permeability fields can be observed. It gets more interesting to see that, by comparing respective histograms, the case with less confidence in the WWCT data results in an almost perfect regeneration of the reference permeability field distribution shown in Fig. 10. This further supports the observation made in the previous example on how important matching the reservoir pressure is

for this model. Comparing the WWCT responses between the base and the Large WWCT data variance cases (Fig. 13) shows increased spread in the ensemble response, consistent with the large data noise specification for the algorithm. Similarly, as shown in Fig. 14, reducing the data noise results in smaller spreads of the well bottomhole pressure responses of the generated ensemble of compromise models. Based on our observation from these numerical experiments, the multiobjective McMC algorithm handles variable data noises are in a consistent manner while accurately capturing the parameters and quantifying the uncertainties in them. 4.5. The Brugge benchmark model Here we further demonstrate the utility and robustness of the multiobjective McMC algorithm in a more realistic subsurface model calibration. Our choice model is the Brugge reservoir (Peters et al., 2010) which is a synthetic benchmark case developed to appraise closed-loop reservoir management workflows and algorithms. The model, shown in Fig. 15, has a half-dome structure with a single fault cutting through the reservoir at one edge. The model was built with 60,048 corner point grid cells with nine layers which follow the underlying stratigraphy of the reservoir. The model consists of four geologic zones: Shelde, Waal, Maas and Schie; each with different flow and storage properties. The reservoir model is made up of 7 sand facies, each with its unique rock-fluid physics. 768

F. Olalotiti-Lawal, A. Datta-Gupta

Journal of Petroleum Science and Engineering 166 (2018) 759–777

Fig. 15. Permeability field of the 103rd realization of the Brugge model used as the prior in our application.

freedom will be necessary to effectively explore all possible solutions. Olalotiti-Lawal and Datta-Gupta (2015) as well as Park et al. (2015) performed the GCT re-parameterization based on each geologic zone. This will however be limiting for solution space exploration purposes since, as previously noted (Bhark et al., 2011), significant loss of vertical heterogeneity might result during permeability updates. It gets more severe with the Brugge model in which the median ratio of vertical to horizontal permeability values is less than 0.01. For this reason, we carry out layer-wise re-parameterization of the model retaining only 9 leading order bases per layer. Our selection of the basis functions for each layer depends on the underlying heterogeneity in the prior model. We followed the procedure

All of these features brings the reservoir model, and thus the history matching problem closer to reality. The reservoir prediction began with early primary depletion, followed by a peripheral waterflood program over the 10 years of history. Here, we apply the multiobjecitve McMC algorithm in probabilistically characterizing all possible updated models that are all constrained to the well bottomhole pressure (WBHP) and water cut (WWCT) data. The 103rd realization of 104 Brugge model realizations was selected as our prior model. Note that the initial model response of this realization happens to be closer to the observed data, unlike the fluvial realization that was used in the earlier version of this paper (Olalotiti-Lawal and Datta-Gupta, 2015). A clear implication of this is that a higher degree of

Fig. 16. Pareto set of solutions obtained for the Brugge case. 769

F. Olalotiti-Lawal, A. Datta-Gupta

Journal of Petroleum Science and Engineering 166 (2018) 759–777

Fig. 17. Samples of multiplier fields from the optimal trade-off solutions.

covariance for well bottomhole pressures and water cut were assumed uncorrelated and therefore took the form of σ 2p I and σ pw I respectively,

outlined by Bhark et al. (2011). For each layer, we projected the permeability field onto a 100 dimensional spectral domain to obtain the corresponding basis coefficients which are sorted in a descending order of magnitude. Then the first 9 corresponding basis functions that follow the same order are selected for that layer. This way, the initial model heterogeneity features are already captured in the basis selection, thereby improving the sampling robustness and eventual geologic realism of updated models. This basis selection routine was not necessary in the synthetic model since we started with a homogeneous permeability as our prior model. For our application, we estimate 81 parameters using 300 Markov chains. This is in compliance with the general guidance provided by ter Braak (2006), regarding the number of chains required for high dimensional problems. The temperature ladder limits for this case were specified equal to that employed for the synthetic model case. However, note here that the numerical differences in temperature between chains are reduced, as there are many more chains compared to the synthetic case. This is desirable and ultimately favors better communication between chains thereby promoting improved mixing. Similarly here, data

where σ 2p ¼ 5000 and σ 2w ¼ 0:01. The algorithm was likewise ran through 100 generations. The sampled Pareto solutions are shown in Fig. 16 in which also color-coded the different classifications of the Pareto solution spectrum. Optimal tradeoff (or compromise) solution region was selected applying the local least square regression method discussed earlier in the synthetic model example. We provide the layer images of the permeability multiplier field samples at the three regions of interest along the Pareto front in Fig. 17 for comparison purposes. While permeability updates in certain layers share some resemblance in pattern, a holistic comparison of models in each partition of the Pareto front indicate the desired diversity of plausible models, all conditioned to production data. When comparing these models, it is important to keep in mind that these maps show layer-wise multiplier fields on the logarithmic scale, which means that slight differences in these values for a layer can result in marked variations in the model response. We present the parameter estimate and quantified uncertainties there-in in

Fig. 18. Normalized distributions of basis coefficients for optimal trade-off models from layer 1 (left end) to layer 9 (right end). 770

F. Olalotiti-Lawal, A. Datta-Gupta

Journal of Petroleum Science and Engineering 166 (2018) 759–777

Fig. 19. Optimal Trade-off model responses.

multiobjective probability distribution function to the description of the Differential Evolution Markov Chain Monte Carlo (DEMC) algorithm which plays a pivotal role in the multiobjective McMC algorithm. We illustrated our algorithm with a simple 2D reservoir example which demonstrated the relative strength of the algorithm compared to a standard NSGA-II algorithm. Although our applications have purely dwelled on biobjective optimization and uncertainty quantification problems, we expect comparable performance with problems with three dimensions in the objective space. The following conclusions can be drawn from this work:

the normalized distribution plot shown in Fig. 18. Note that in all these plots, the prior parameter is distributed between 0 and 1. Also in the updated parameter distribution plots, the red line in each box plot represents the median value, while the left and right edges of each box plot represent the 25th and 75th percentiles respectively for the corresponding parameter. The parameter distribution provided here is only for the optimal trade-off solution ensemble. Similar plots are provided for the other two extremes of the Pareto front in Appendix C. It is easy to notice relatively larger range of parameter values in the compromise solution set, compared with the other two extremes, especially for layers 1 and 7 in which the first (high energy) bases appear to be poorly resolved. This is also reflected in Fig. 17 where a marked difference exists between the first layer multiplier fields of the two samples shown. Integration of other available data such as seismic might also help to resolve these parameters. Model responses from the optimal trade-off models are shown in Fig. 19. Similar plots are provided in Appendix C for the other extreme ensembles. Although the matches are not as perfect as those obtained in the synthetic model discussed earlier, comparing these plots, it is clear to see the compromise in the model matches shown in Fig. 19. Quality of a history match hinges on the re-parameterization approach, which in our case is the GCT which significantly reduced the problem dimensionality. For the sake of simplicity, we only parameterized the permeability field in all of our applications here, however a more sophisticated approach will be to also re-parametrize the reservoir structure using GCT-based pore volume multiplier field, together with the GCT-based permeability multiplier field (Kam et al., 2016). Since introducing the reservoir pore volume calibration smoothly modifies the model storativity which the reservoir pressure is highly impacted by; we surmise that this will result in a significant improvement of the pressure match in the Brugge model example.

1 The Multiobjective McMC algorithm can serve as a viable option for sampling multiple trade-off solutions along the Pareto front, just like most classical multiobjective optimization algorithm. However, the approach allows for quantifying uncertainties in the process as well. In other words, each of the trade-off solutions obtained with our algorithm carries with it the range of underlying uncertainty, which often result from uncertainties in the individual objective functions (data sets in a history matching problem). 2 The hypervolume-based dominance criterion in a multiobjective optimization approach helps to improve convergence compared to the rank-based criterion in a classical NSGA-II algorithm. The strength of the proposed algorithm comes from combining the hypervolumebased dominance relationship with an efficient search algorithm like the DEMC. In our implementation, we adopted a hybrid parallel tempering methodology which further improves the sampling efficiency of the multi-modal posterior distribution. 3 Adopting an McMC approach in sampling intricate multidimensional probability distribution functions, like the one constructed for the probabilistic multiobjective optimization problem requires improved mixing among Markov chains. The DEMC provides the needed level of mixing for accurate sampling compared to independent parallel Markov chains.

5. Concluding remarks We presented a Multiobjective Markov chain Monte Carlo Algorithm for subsurface model calibration and uncertainty quantification. We began by introducing the concept of probabilistic multiobjective optimization and how it contrasts with the deterministic perspective. Next we described the algorithm, starting from the hypervolume-based

Acknowledgement The authors acknowledge the sponsors of the Model Calibration and Efficient Reservoir Imaging (MCERI) JIP at Texas A&M University.

771

F. Olalotiti-Lawal, A. Datta-Gupta

Journal of Petroleum Science and Engineering 166 (2018) 759–777

Nomenclature d Δ½:

η

f f HYPð:Þ m Φ

π ρ

T u v x

Observation data vector Crowding distance function Quantified dominance relationship Single objective function Objective function vector Hypervolume function Reservoir parameter vector GCT Basis functions Sampled probability distribution function Acceptance probability Simulated Temperature Grid cell property vector Basis coefficient vector Solution vector

Operators/Symbols 8 For all 9 There exists ^ And  Dominates ∅ Null set Subscripts I; J i R1; R2 obs

Index of solutions points in parameter space Index of objective function Randomly selected Markov chains Observed

Superscripts I; J Index of solutions points in objective space t State of Markov chain Appendix A. Comparing DEMC with Parallel Markov Chains Here we illustrate the superior sampling effectiveness of the Differential Evolution Markov Chain Monte Carlo (DEMC) compared with independent Parallel Markov Chain Monte Carlo (PMCMC) algorithms. Both DEMC and PMCMC involve running multiple Markov chains in parallel through the parameter space. The chains in DEMC however exchange information as they sample the space but chains are completely independent in the PMcMC algorithm. Running a traditional Markov chains in parallel was mainly for convergence diagnostics (Maucec et al., 2007), which means conclusions from this experiment also hold for single Markov chains. To be fair, several cases of traditional McMC sampling of varying step size using the Random Walk Metropolis proposal generation were conducted. The step size which gave the best performance was compared with the DEMC algorithm. We compare the two algorithms by evaluating their performance in sampling a 6-D bimodal probability distribution function (PDF). The distribution was generated by summing up two separate 6-D Gaussian distributions defined as Nðμ1 ; σ1 Þ and Nðμ2 ; σ2 Þ; where μ1 ¼  μ2 ¼ ð7; 7; 7; 7; 7; 7Þ, σ1 ¼ 1:5I6 and σ2 ¼ 2:5I6 . Both algorithms were executed with the same set of parameters declared in Table A1. Table A1 Algorithm Parameters used for DEMC and PMCMC samplers Algorithm Parameter

DEMC

PMCMC

Number of Chains Maximum Iterations Per Chain Burn in time Total Number of Function Evaluations

12 2000 500 24,000

12 2000 500 24,000

The samples collected from both algorithms are analyzed and presented as follows. The comparisons between the samples with the true PDFs using both algorithms is shown in Fig. A1, while Fig. A2 shows the same comparisons, but with the Cumulative Distribution Function (CDF). The DEMC algorithm clearly obtained samples more representative of the true PDF than PMcMC. The PMcMC, like traditional McMC algorithms, tends to perform poorly sampling disconnected PDFs like this since adequate mixing is not guaranteed. The experiment was repeated 500 times while computing the misfits in the CDFs obtained from both algorithms. Fig. A3 presents the box plot of the CDF misfits, clearly showing a superior sampling performance using the DEMC algorithm.

772

F. Olalotiti-Lawal, A. Datta-Gupta

Journal of Petroleum Science and Engineering 166 (2018) 759–777

Fig. A1. Comparing true (red line) and estimated (green bars) PDFs using (a) DEMC and (b) PMcMC samplers.

Fig. A2. Comparing true (in blue) and estimated (in red) CDFs using (a) DEMC and (b) PMCMC samplers.

773

F. Olalotiti-Lawal, A. Datta-Gupta

Journal of Petroleum Science and Engineering 166 (2018) 759–777

Fig. A3. Comparing over CDF misfits between (a) DEMC and (b) PMcMC samplers over 500 trials.

Appendix B. BHP and WWCT Match Model Reponses for Synthetic Case

Fig. B1. Best Well Water Cut (WWCT) match sample responses for (a) Well Water Cut and (b) bottomhole pressure.

774

F. Olalotiti-Lawal, A. Datta-Gupta

Journal of Petroleum Science and Engineering 166 (2018) 759–777

Fig. B2. Best bottomhole pressure (BHP) match sample responses for (a) Well Water Cut and (b) bottomhole pressure.

Appendix C. BHP and WWCT Match Model Reponses for the Brugge Model

Fig. C1. Normalized distributions of basis coefficients for WBHP match models from layer 1 (left end) to layer 9 (right end).

775

F. Olalotiti-Lawal, A. Datta-Gupta

Journal of Petroleum Science and Engineering 166 (2018) 759–777

Fig. C2. WBHP match model responses.

Fig. C3. Normalized distributions of basis coefficients for WWCT match models from layer 1 (left end) to layer 9 (right end).

Fig. C4. WWCT match model responses.

776

F. Olalotiti-Lawal, A. Datta-Gupta

Journal of Petroleum Science and Engineering 166 (2018) 759–777

References

Mirzabozorg, A., Nghiem, L., Chen, Z., et al., 2013. Differential evolution for assisted history matching process: sagd case study. In: Proc., SPE Heavy Oil Conference Canada. Canada: Society of Petroleum Engineers, Calgary, Alberta. Mohamed, L., Christie, M.A., Demyanov, V., et al., 2010. Application of Particle Swarms for history matching in the Brugge reservoir. In: Proc., SPE Annual Technical Conference and Exhibition Society of Petroleum Engineers. Olalotiti-Lawal, F., 2013. Application of Fast Marching Methods for Rapid Reservoir Forecast and Uncertainty Quantification. Texas A&M University. Olalotiti-Lawal, F., Datta-Gupta, A., 2015. A multi-objective Markov chain Monte Carlo approach for history matching and uncertainty quantification. In: Proc., SPE Annual Technical Conference and Exhibition Society of Petroleum Engineers. Oliver, D.S., Cunha, L.B., Reynolds, A.C., 1997. Markov chain Monte Carlo methods for conditioning a permeability field to pressure data. Math. Geol. 29 (1), 61–91. Omre, H., Tjelmeland, H., 1997. Petroleum geostatistics. Geostat. Wollongong 96, 41–52. Park, H.-Y., Datta-Gupta, A., King, M.J., 2015. Handling conflicting multiple objectives using pareto-based evolutionary algorithm during history matching of reservoir performance. J. Petrol. Sci. Eng. 125, 48–66. Peters, L., Arts, R., Brouwer, G., et al., 2010. Results of the Brugge benchmark study for flooding optimization and history matching. SPE Reserv. Eval. Eng. 13 (03), 391–405. Robert, C., Casella, G., 2013. Monte Carlo Statistical Methods. Springer Science & Business Media. Romero, C., Carter, J., Gringarten, A., et al., 2000. A modified genetic algorithm for reservoir characterisation. In: Proc., International Oil and Gas Conference and Exhibition in China Society of Petroleum Engineers. Saleri, N., Toronyi, R., 1988. Engineering control in reservoir simulation: Part I. In: Proc., SPE Annual Technical Conference and Exhibition Society of Petroleum Engineers. Schulze-Riegert, R.W., Krosche, M., Fahimuddin, A., et al., 2007. Multi-objective optimization with application to model validation and uncertainty quantification. In: Proc., SPE Middle East Oil and Gas Show and Conference Society of Petroleum Engineers. Singh, A., Minsker, B.S., 2008. Uncertainty-based multiobjective optimization of groundwater remediation design. Water Resour. Res. 44 (2). Storn, R., Price, K., 1997. Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces. J. Glob. Optim. 11 (4), 341–359. Tanaka, S., Kam, D., Datta-Gupta, A., et al., 2015. Streamline-based history matching of arrival times and bottomhole pressure data for multicomponent compositional systems. In: Proc., SPE Annual Technical Conference and Exhibition Society of Petroleum Engineers. ter Braak, C.J., 2006. A Markov chain Monte Carlo version of the genetic algorithm differential evolution: easy bayesian computing for real parameter spaces. Stat. Comput. 16 (3), 239–249. ter Braak, C.J., Vrugt, J.A., 2008. Differential evolution Markov chain with snooker updater and fewer chains. Stat. Comput. 18 (4), 435–446. Vousden, W., Farr, W.M., Mandel, I., 2016. Dynamic temperature selection for parallel tempering in Markov chain Monte Carlo simulations. Mon. Not. R. Astron. Soc. 455 (2), 1919–1937. Wang, C., Gao, J., Yang, H., et al., 2011. Waveform inversion of cross-well data with cooperative coevolutionary differential evolution algorithm. In: Proc., 2011 SEG Annual Meeting Society of Exploration Geophysicists. While, L., Hingston, P., Barone, L., et al., 2006. A faster algorithm for calculating hypervolume. Evol. Comput. IEEE Trans. 10 (1), 29–38. W€ ohling, T., Vrugt, J.A., 2008. Combining multiobjective optimization and bayesian model averaging to calibrate forecast ensembles of soil hydraulic models. Water Resour. Res. 44 (12). Xie, J., Efendiev, Y., Datta-Gupta, A., 2011. Uncertainty quantification in history matching of channelized reservoirs using Markov chain level set approaches. In: Proc., SPE Reservoir Simulation Symposium Society of Petroleum Engineers. Yin, J., Park, H.-Y., Datta-Gupta, A., et al., 2011. A hierarchical streamline-assisted history matching approach with global and local parameter updates. J. Petrol. Sci. Eng. 80 (1), 116–130. Zhang, Q., Zhou, A., Jin, Y., 2008. Rm-meda: a regularity model-based multiobjective estimation of distribution algorithm. Evol. Comput. IEEE Trans. 12 (1), 41–63. Zitzler, E., Thiele, L., 1999. Multiobjective evolutionary algorithms: a comparative case study and the strength Pareto approach. Evol. Comput. IEEE Trans. 3 (4), 257–271.

Aanonsen, S.I., Nævdal, G., Oliver, D.S., et al., 2009. The ensemble kalman filter in reservoir engineering–a review. Spe J. 14 (03), 393–412. Bader, J., Zitzler, E., 2011. Hype: an algorithm for fast hypervolume-based manyobjective optimization. Evol. Comput. 19 (1), 45–76. Bhark, E.W., Jafarpour, B., Datta-Gupta, A., 2011. A generalized grid connectivity–based parameterization for subsurface flow model calibration. Water Resour. Res. 47 (6). Castellini, A., Gullapalli, I., Hoang, V.T., et al., 2005. Quantifying uncertainty in production forecast for fields with significant history: a West African case study. In: Proc., International Petroleum Technology Conference International Petroleum Technology Conference. Cheng, H., Oyerinde, A.S., Datta-Gupta, A., et al., 2007. Compressible streamlines and three-phase history matching. SPE J. 12 (04), 475–485. Christie, M., Eydinov, D., Demyanov, V., et al., 2013. Use of multi-objective algorithms in history matching of a real field. In: Proc., SPE Reservoir Simulation Symposium Society of Petroleum Engineers. Das, I., Dennis, J.E., 1997. A closer look at drawbacks of minimizing weighted sums of objectives for Pareto set generation in multicriteria optimization problems. Struct. Optim. 14 (1), 63–69. Deb, K., Pratap, A., Agarwal, S., et al., 2002. A fast and elitist multiobjective genetic algorithm: nsga-Ii. Evol. Comput. IEEE Trans. 6 (2), 182–197. Dong, Y., Oliver, D.S., 2008. Reservoir simulation model updates via automatic history matching with integration of seismic impedance change and production data. In: Proc., International Petroleum Technology Conference International Petroleum Technology Conference. Earl, D.J., Deem, M.W., 2005. Parallel tempering: theory, applications, and new perspectives. Phys. Chem. Chem. Phys. 7 (23), 3910–3916. Fern
andez-Martínez, J., García-Gonzalo, M., Fern
andez-Alvarez, J., et al., 2008. Particle Swarm optimization (pso): a simple and powerful algorithm family for geophysical inversion. In: Proc., 2008 SEG Annual Meeting Society of Exploration Geophysicists. Geyer, C.J., Thompson, E.A., 1995. Annealing Markov chain Monte Carlo with applications to ancestral inference. J. Am. Stat. Assoc. 90 (431), 909–920. Hajizadeh, Y., Christie, M.A., Demyanov, V., 2010. History matching with differential evolution approach; a look at new search strategies. In: Proc., SPE EUROPEC/EAGE Annual Conference and Exhibition Society of Petroleum Engineers. Hajizadeh, Y., Christie, M.A., Demyanov, V., 2011. Towards multiobjective history matching: faster convergence and uncertainty quantification. In: Proc., SPE Reservoir Simulation Symposium Society of Petroleum Engineers. Haldorsen, H.H., 1986. Simulator parameter assignment and the problem of scale in reservoir engineering. Reserv. Charact. 6. Hastings, W.K., 1970. Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57 (1), 97–109. Izui, K., Yamada, T., Nishiwaki, S., et al., 2015. Multiobjective optimization using an aggregative gradient-based method. Struct. Multidiscip. Optim. 51 (1), 173–182. Kam, D., Han, J., Datta-Gupta, A., 2016. Streamline-based rapid history matching of bottomhole pressure and three-phase production data. In: Proc., SPE Improved Oil Recovery Conference Society of Petroleum Engineers. Li, Y., 2012. Momcmc: an efficient Monte Carlo method for multi-objective sampling over real parameter space. Comput. Math. Appl. 64 (11), 3542–3556. Liu, X., Reynolds, A.C., 2016. Gradient-based multi-objective optimization with applications to waterflooding optimization. Comput. Geosci. 20 (3), 677–693. Ma, X., Al-Harbi, M., Datta-Gupta, A., et al., 2008. An efficient two-stage sampling method for uncertainty quantification in history matching geological models. SPE J. 13 (01), 77–87. Marinari, E., Parisi, G., 1992. Simulated tempering: a new Monte Carlo scheme. EPL Europhys. Lett. 19 (6), 451. Maucec, M., Douma, S.G., Hohl, D., et al., 2007. Streamline-based history matching and uncertainty–markov-chain Monte Carlo study of an offshore turbidite oil field. In: Proc., SPE Annual Technical Conference and Exhibition Society of Petroleum Engineers. Metropolis, N., Ulam, S., 1949. The Monte Carlo method. J. Am. Stat. Assoc. 44 (247), 335–341.

777