A hybrid differential evolution algorithm approach towards assisted history matching and uncertainty quantification for reservoir models

Journal of Petroleum Science and Engineering 142 (2016) 21–35

Contents lists available at ScienceDirect

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

A hybrid differential evolution algorithm approach towards assisted history matching and uncertainty quantification for reservoir models Emil C. Santhosh, Jitendra S. Sangwai n Gas Hydrate and Flow Assurance Laboratory, Petroleum Engineering Program, Department of Ocean Engineering, Indian Institute of Technology Madras, Chennai 600036, India

art ic l e i nf o

a b s t r a c t

Article history: Received 19 May 2015 Received in revised form 19 January 2016 Accepted 29 January 2016 Available online 30 January 2016

History matching is an important process in the reservoir model development. In the process of history matching, the most significant uncertain model parameters are identified and adjusted to get an acceptable match between the simulated production with the historical field production data. In the past decade, many population based algorithms have been applied for history matching. In this paper, a novel population based stochastic algorithm called hybrid differential evolution (HDE) is applied for the assisted history matching process. An adaptive mechanism for the control parameters is incorporated in the algorithm which automatically adjusts the control parameters according to the problem. The performance of the algorithm is tested on a 3-D reservoir model called PUNQ-S3 which is a benchmark model for the comparison of different history matching and uncertainty quantification techniques. Since history matching is an inverse problem, multiple models can give good match. So, prediction using a single history matched model involves more risk because of the parameter uncertainty. One of the methods to solve this problem is to quantify the uncertainty in the predictions. In this paper, the neighbourhood approximation Bayes (NAB) algorithm is applied to quantify the uncertainty in reservoir forecast which is a Bayesian extension of neighbourhood algorithm. The NAB algorithm quantifies the uncertainty in the predictions using multiple models generated during history matching phase and this does not require additional simulations. The main focus of this paper is to study about how HDE algorithm can be used when coupling with the NAB algorithm in predicting the true forecast with minimum uncertainty range under limited number of simulations. The influence of population size on the performance of the algorithm in history matching and forecast is analyzed. The HDE provides wide sampling of the search space and the truth case was comfortably included within the predicted confidence bounds. The results show that HDE can be used as a promising tool for assisted history matching of the reservoir models. & 2016 Elsevier B.V. All rights reserved.

Keywords: Control parameters History matching Hybrid differential evolution Reservoir simulation Uncertainty quantification

1. Introduction In modern oilfield practices, reservoir simulation plays an important role because of its ability to predict the performance of the reservoir and helps to prepare optimal oilfield development plans. A good knowledge about the reservoir is required to build an effective reservoir model which can mimic the real reservoir. The petroleum reservoirs are usually geologically complex and highly heterogeneous in nature. Our knowledge about the underground reservoir is limited due to the sparse and low resolution data obtained from seismic and other well-log and well-testing methods. This leads to uncertainties in the value of parameters used in the reservoir model. In modern reservoir engineering practices, n

Corresponding author. E-mail address: [email protected] (J.S. Sangwai).

http://dx.doi.org/10.1016/j.petrol.2016.01.038 0920-4105/& 2016 Elsevier B.V. All rights reserved.

once the reservoir model is developed with the available data, it has to get validated before using it for the prediction of reservoir fluid flows and other reservoir properties. This can be done by comparing the simulation results with the oilfield production history. In history matching, the reservoir model is conditioned to the historical observations from the field such as oil production rate, water cut, bottom-hole pressure, etc. The process of history matching in reservoir modeling started almost six decades ago. In early days, the history matching started as a manual method in which the user manually adjusts the parameters until a good match between the reservoir predictions and the actual production profiles is found (Sheldon et al., 1960; Jacquard, 1965). However, this is very time consuming and intuition-based process which may not always provide a global solution. In 1972, there came one of the first attempts of automatic history matching (Thomas et al., 1972). The authors considered the history matching problem as an optimal control problem.

22

E.C. Santhosh, J.S. Sangwai / Journal of Petroleum Science and Engineering 142 (2016) 21–35

Following this, several authors have tried different optimization techniques for history matching. The main methods include gradient based methods, Monte Carlo methods, ensemble Kalman filter, etc. Gradients based methods showed less performance in history matching because of their tendency to get trapped in local minima. In stochastic methods, many authors have tried different algorithms like simulated annealing (SA) (Sultan et al., 1994), neighbourhood algorithm (NA) (Subbey et al., 2003), genetic algorithms (GA) (Castellini et al., 2005; Sangwai et al., 2007), scatter search (SS) (Sousa et al., 2006; Erbas and Christie, 2007), Markov chain Monte Carlo (McMC) (Maucec et al., 2007), particle swarm optimization (PSO) (Mohamed et al., 2009), ant colony optimization (Hajizadeh et al., 2011), differential evolution (DE) (Hajizadeh et al., 2010). It is observed that the differential evolution algorithm has shown good results for the history matching but the performance of the algorithm was very much sensitive to the value of control parameters such as crossover rate. In this work, a hybrid differential evolution (HDE) (ReynosoMeza et al., 2011) algorithm is being introduced for the history matching problem in which the algorithm itself adapts the best value for the control parameter according to the type of the problem. Since the history matching is an inverse problem, multiple models may show good match with the oilfield data. So, planning reservoir development decisions based on a single reservoir model carries more risk. This can be solved by taking information from multiple history matched models. The uncertainty in the predictions can be quantified from the posterior probability distribution of the models. In this paper, hybrid differential evolution algorithm is applied to the history matching of a 3-D petroleum reservoir called PUNQ-S3 which is a benchmark model for the comparison of different history matching and uncertainty quantification techniques. The uncertainty in the oil production forecast is quantified using neighbourhood approximation Bayes (NAB) algorithm (Sambridge, 1999a, 1999b) which is Bayesian extension of neighbourhood algorithm. The approach taken in this work for assisted history matching and uncertainty quantification using HDE algorithm and NAB routine is shown in Fig. 1. The data

collected through several field tests and measurements gives us some idea about the reservoir. It helps us to get information about the initial ranges of uncertain parameters. With this, the history matching loop is initiated in which the misfit between the simulated production data and field production data is minimized using the hybrid differential evolution (HDE). The history matching loop will continue until the termination criteria are met which is either an acceptable misfit or maximum number of simulations. This step will generate an ensemble of history matched models. This ensemble of models is submitted to the inference step. The uncertainty in reservoir predictions is quantified from the posterior probability distributions of the models calculated using NAB algorithm.

2. Theory and methodology 2.1. Differential evolution Recently, differential evolution has gained much importance because of its efficacy in solving real parameter optimization problems (Das and Suganthan, 2010). It is a powerful population based stochastic optimization algorithm with three operators which are: 2.1.1. Mutation In differential mutation, the difference between two vectors is multiplied by a constant which is called scaling factor and this is added to a third vector in the population. This will produce a new mutant vector which makes some perturbation in the models in each generation. Based on the selection of vectors for mutation makes different variants for DE.

(

miG = p1G + F × p2G −p3G

)

miG

(1)

where, is the mutant vector generated, p is the parent vector in the generation G. F is called the scaling factor and is always positive and ranges between 0 and 2. The magnitude of the

Fig. 1. Approach taken for history matching and uncertainty quantification.

E.C. Santhosh, J.S. Sangwai / Journal of Petroleum Science and Engineering 142 (2016) 21–35

difference vector will depend upon the scaling factor. The selection of parents used in this paper is random. 2.1.2. Crossover There are different types of crossover operators. The main idea behind this operator is to make diversity in the population. In this operator, the mutant vector (m) is mixed with the parent vector. One of the main crossover operators is the binomial crossover.

⎧ ⎪ miG uiG = ⎨ G ⎪ ⎩ pi

if

( rand (j) ≤ Cr )

if (rand (j ) > Cr )

(2)

where, uiG is the trial vector produced after the crossover operator, Cr is the crossover probability in the range of [0,1] which controls the fraction of parameter values that are copied from the mutant vector. All parent vectors are uniquely indexed from 1 to Np (population size) and the generated mutant vectors are also indexed from 1 to Np such that the mutant vector 1 ( m1G ) will be generated using the parent vector 1 ( p1G ) as the base vector and so on. The mutant vector having the same index number as that of the parent vector will be compared at this step. For the jth parameter value of the vector uiG , a random variable rand(j) is generated in the range of [0,1] and the parameter value of the parent vector is replaced by its pair in the mutant vector if rand(j) is less than or equal to Cr. 2.1.3. Selection New individuals are selected based on the fitness value. The child will replace the parent vector if the child vector has a better fitness value than the parent vector.

⎧ u G if f u G ≤ f p G ⎪ i i i piG + 1 = ⎨ G ⎪ else ⎩ pi

( )

( ) (3)

where, piG + 1 is the vector selected to the next generation and f is the fitness function.

23

2.2. Hybrid differential evolution Some of the authors have applied differential evolution algorithm for history matching process (Wang and Buckley, 2006; Hajizadeh et al., 2010). They observed that the performance of the differential evolution algorithm is very sensitive to the control parameters such as crossover rate and selecting the proper control parameter for a particular problem is an issue. For example, lower values of crossover rate are recommended for separable problems and higher values for non-separable problems (Storn, 2008). A separable problem is the one in which all of its variables are separable and a non-separable problem is the one in which some of its variables are correlated. Here, hybrid differential evolution algorithm (Reynoso-Meza et al., 2011) is applied for the history matching problem which is incorporated with a local search technique and an adaptive mechanism for the crossover rate. The main idea is that a proper balance of the exploration abilities of the DE and exploitation abilities of a local searcher (LS) can lead to an algorithm with better performance. The adaptive mechanism combines the binary crossover and the line recombination. To avoid stagnation, a population refreshment mechanism is also incorporated in the algorithm. Fig. 2 shows the pseudo code for hybrid differential evolution algorithm used in this work. 2.2.1. Local search strategy Local search strategies have improved the convergence speed of the DE algorithm (Noman and Iba, 2008; Neri and Tirronen, 2009; LaTorre et al., 2010). The local search mechanism incorporated to improve the convergence is a sequential quadratic programming (SQP) routine. Sequential Quadratic Programming (SQP) (Reklaitis et al., 1983; Boggs and Tolle, 2000) is a gradient based non-linear optimization method. The method is identical to Newton's method for constrained optimization. The SQP method seems to be the best nonlinear programming method for constrained optimization. It shows good performance in terms of accuracy and efficiency when compared with other nonlinear programming methods (Boggs and Tolle, 2000). The SQP routine is executed in every child of the population with a probability αLS. If

Fig. 2. Pseudo code for HDE algorithm used for history matching. F: scaling factor; Cr: crossover operator.

24

E.C. Santhosh, J.S. Sangwai / Journal of Petroleum Science and Engineering 142 (2016) 21–35

a random number generated during a generation is more than that of the αLS, algorithm will perform a local search. Details of the SQP procedure are presented by Fletcher (1987). 2.2.2. Adaptive parameter values The values for scaling factor (F) and crossover rate (Cr) of a generation are calculated from the triangular distributions ΛF, ΛCr respectively from

If rand < (MedCr − MinCr)/(MaxCr − MinCr) Cr = MinCr +

(rand*(MaxCr − MinCr)*(MedCr − MinCr)

else Cr = MaxCr− (1−rand)*(MaxCr − MinCr)*(MaxCr − MedCr) end

if it shows a good match within the history matching period. One of the ways to solve this problem is to calculate the production forecast within the confidence bound using the information collected from all the models developed during the history matching phase. A model which shows a poor match will also convey some amount of information as the other models do. So, it is important to take information from all the models in the ensemble. The uncertainty in production forecast is quantified by computing the posterior mean prediction and the posterior confidence intervals around the mean from the posterior probability distribution (PPD) of the models. The posterior probability distribution of the models is calculated using neighbourhood approximation Bayes algorithm (NAB) as described below.

(4)

where rand is a random number. MinCr, MedCr, MaxCr indicates the minimum, median and maximum values of the triangular distribution ΛCr. The same procedure is adopted for the calculation of the scaling factor. The triangular distribution ΛF is fixed and it is used to make some dither in the algorithm. The triangular distribution ΛCr is incorporated with an adaptive mechanism in which the triangular distribution ΛCr is adapted after a number of recorded successes. If the child vector is having a lesser or equal fitness value than that of parent vector, then parent vector will be substituted by the child which is recorded as a success. The values for the triangular distribution are taken from the minimum, medium and maximum values of such set of recorded successes. The main purpose of the adaptive mechanism is to identify the problem type. The adaptive mechanism identifies the problem type then chooses low Cr values for separable problems and high Cr values for non-separable problem for the binomial crossover operator. Upon detecting strong-dependency among decision variables, the algorithm uses a non-rotationally invariant line recombination (Price, 1999). 2.3. Exploration and exploitation Exploitation is the ability of the optimization algorithm to search in a promising region of the search space in order to improve a solution that we already have while exploration is the ability of the optimization algorithm to search different areas of the search space in order to find promising solutions that are yet to be refined. The exploration relates to diversification which indicates how different the individuals within a population are. The algorithm should have a good diversification in order to explore the search space properly. A low level of diversification means a convergence. So, it is necessary for the algorithm to have high diversification at the beginning of the search process and a lower one at the end. A lower level of diversification may lead to premature convergence towards the local optima. The exploration capability of the algorithm helps to diversify the search thus minimizing the tendency to get trapped in local optima. Zaharie (2001) provides details on the explorative power of differential evolution algorithms. The details about the critical values for the control of differential evolution algorithms and their associated parameter tunings are discussed in (Zaharie, 2002, 2003). The details about balancing the exploration and exploitation capabilities of the differential evolution algorithm are given in Epitropakis et al. (2008). 2.4. Uncertainty quantification using NAB algorithm The history matching is an ill-posed inverse problem. Since illposed inverse problems have non-unique solutions, we cannot rely on a single history matched model for a good prediction even

2.4.1. Neighbourhood approximation Bayes algorithm (NAB) NAB is a Markov chain Monte Carlo (McMC) method using Gibbs sampler. The NAB algorithm was originally developed for the inversion of receiver functions for crustal seismic structure (Sambridge, 1999a, 1999b). In NAB algorithm, the model search space is represented using Voronoi cells. Each point in the model search space corresponds to a model in the ensemble and the Voronoi cells are the neighbor regions about each point in the model space. The size and shape of the Voronoi cells are defined by a distance norm. The volume of Voronoi cell is inversely proportional to the density of points in the model search space. Each Voronoi cell represents the area of influence of a particular model in the ensemble. The misfit value inside each Voronoi cell is taken as constant and based on this assumption the PPD of unknown points in the model space is interpolated. The NAB algorithm generates new sets of points (models) whose distribution tends towards posterior probability distribution from an irregular distributed ensemble. In this case the irregular distributed ensemble is generated using the HDE algorithm. The newly generated points are called the resampled ensemble. The resampled ensemble is generated by using a standard statistical technique called Gibbs sampler (Geman and Geman, 1984; Smith and Roberts, 1993). The Gibbs sampler importance samples the irregular distributed ensemble to generate the resampled ensemble. The Gibbs sampler will generate a random walk in the model space and generate models after many iterations which asymptotically tends towards the posterior probability distribution (Gelman et al., 1995). The steps for importance sampling the search space using Gibbs sampler in a 2-dimensional space is explained in the section 3.3 of Sambridge (1999b). The advantages of using NAB are: 1. The posterior probability is evaluated from all the models instead of a subset from the ensemble, thus models with different misfit will contribute differently to the uncertainty prediction. 2. The cost for further forward simulations for all the models generated by sampling algorithm is eliminated by approximating constant misfit values inside the Voronoi cells.

3. Reservoir description and parameterization The performance of HDE in history matching is tested on the PUNQ-S3 reservoir model. The PUNQ-S3 is a synthetic benchmark model which is used for comparing the performance of different history matching and uncertainty quantification techniques (PUNQ-S3, 2015). The same has been used by various researchers for history matching and uncertainty quantification studies (Floris et al., 2001; Demyanov et al., 2004; Holmes et al., 2007; Hajizadeh et al., 2010).

E.C. Santhosh, J.S. Sangwai / Journal of Petroleum Science and Engineering 142 (2016) 21–35

25

Fig. 3. PUNQ-S3 reservoir model showing top surface and well positions.

3.1. PUNQ-S3 reservoir description The PUNQ-S3 case has been taken from a reservoir engineering study on a real field performed by Elf Exploration Production (PUNQ-S3, 2015). PUNQ stands for Production forecasting with Uncertainty Quantification. PUNQ-S3 is a 5-layer reservoir model with a dip angle of 1.5°. The east and south of the reservoir is surrounded by a fault and the west and north are bounded by a strong aquifer. Six production wells are drilled in the reservoir (see Fig. 3). As the aquifer provides good pressure support, there are no injection wells drilled in the reservoir. There is a small gas cap in layer one and no wells are drilled in this area to prevent free gas production and reduction in the reservoir pressure. The reservoir is having a top depth of 2430 m. The 3 production wells (PRO-1, PRO-4 and PRO-12) are perforated in the layers 4 and 5. The other producer wells (PRO-5 and PRO-11) are perforated in the layers 3 and 4 and PRO-15 is only perforated in the layer 5. The PUNQ-S3 reservoir is modeled using a black oil simulator (Schlumberger, 2010) which consisting of 19
28
5 grid cells which become a total of 2660 grid cells out of which only 1761 grid cells are active. The grid blocks are having dimensions of 180 m in both x and y direction. The geological model is modeled by corner point geometry. The production history for the first 8 years is available which includes well oil production rate (WOPR), well water cut (WWCT), bottom hole pressure (BHP) and gas–oil ratio (GOR) of all the six wells. Since measurement errors are unavoidable, Gaussian noise is incorporated into the well porosity/permeability data and to the synthetic field production data (PUNQ-S3, 2015). The history includes total of eight years with 1 year of extended well testing followed by 3 years of field shut-in, and covers 4 years of actual field production at 150 sm3 /day. After eight years of production there are two production strategies. First one is to continue the same production strategy for the next 8.5 years with the same six production wells and the second strategy is to add 5 more wells and continuing the production for

another 16.6 years. In this work, only the first strategy is followed. The data set for PUNQ is available online (PUNQ-S3, 2015). The porosity and permeability distribution of the reservoir are unknown. The PUNQ-S3 reservoir model with top surface map and well positions are shown in Fig. 3. 3.2. PUNQ-S3 parameterization The PUNQ-S3 model consists of 2660 grid cells and each grid cell contains 4 unknown parameters. The unknown parameters are porosity and permeability in x, y and z directions. That will become 10,640 unknown parameters in total. With techniques like sensitivity analysis, the most sensitive parameters in a history matching problem are identified and sometimes pair-wise and triple-wise combinations often yield better impact than separate changes. Dogru and John (1981) provides information regarding the comparison of sensitivity coefficient calculation methods in automatic history matching. The advanced methods to find the most sensitive parameters (decision parameters) in the history matching study are beyond the scope of this paper. In this study, we use the parameterization of PUNQ-S3 described by Hajizadeh et al. (2010). The 5 layered reservoir is parameterized into 9 uniform regions per layer. The correlations between the porosity and horizontal permeability and between horizontal permeability and vertical permeability are obtained from the study performed by Total (Elf Exploration) (Boss, 1999).

ln (kh ) = 0.77 + 9.03φ

(5)

k v = 3.124 + 0.306kh

(6)

The horizontal permeabilities in x and y directions are taken as same. Thus, the total number of unknown parameters is downscaled to 45. The initial ranges for the unknown parameters are shown in Table 1. These values are obtained from the geological description of the reservoir (Hajizadeh et al., 2010, 2011).

26

E.C. Santhosh, J.S. Sangwai / Journal of Petroleum Science and Engineering 142 (2016) 21–35

Table 1 Initial ranges for unknown reservoir parameters of PUNQ-S3 reservoir model (Hajizadeh et al., 2010, 2011).

Table 3 Tuning parameters and best misfit for DE algorithm (Hajizadeh et al., 2010). Cases

S. No. Layer Porosity

1 2 3 4 5

1 2 3 4 5

Horizontal permeability (mD)

0.15–0.3 133–3013 0.05–0.15 16–133 0.15–0.3 133–3013 0.1–0.2 47–376 0.15–0.3 133–3013

Vertical permeability (mD)

DE-Best 50 DE-Rand 50

44–925 8–44 44–925 17–118 44–925

4. Mathematical formulation of the history matching problem In assisted history matching, the process is considered as an optimization problem where the misfit between the field production data and simulated data is minimized. The PUNQ-S3 model has 6 wells and each well is having 4 data types which are WOPR, WBHP, WWCT and WGOR. By using a misfit definition which indicates the difference between field production data and simulated data, the problem is converted to a single objective optimization problem. The misfit function is defined as

1 M= nw

∑i

1 np

∑j

1 nt

∑k

⎛ Oij t k − Sij t k ⎜ w ijk ⎜ σijk ⎝

( )

Population size, NP

2

( ) ⎞⎟ ⎟ ⎠

(7)

where, M represents the mismatch between the simulated production data and the observed production data, here M is used as the fitness function f as described in the Section 2.1.3. nw , np , nt are the number of wells, number of types of production data and number of samples, respectively. i, j, k runs over the number of wells, production types and the samples respectively. O represents the observed data and S represents the simulated data. σ is the measurement error and W is the weight for different data types at different time steps.

5. Results and discussion 5.1. History matching of PUNQ-S3 The hybrid differential evolution algorithm is integrated with the black oil reservoir simulator (Schlumberger, 2010) as shown in the Fig. 1. One of the aims of this study is to analyze the influence of population size of HDE algorithm in history matching and uncertainty quantification. 5.1.1. Influence of population size of HDE in history matching We took many cases of HDE algorithm for history matching by altering the population size. To check the influence of population size, the total number of simulations for all the cases is kept Table 2 Population vs. best misfit for HDE algorithm for this work. Cases

Population size, Np

Number of simulations

Best misfit value

HDE-1 HDE-2 HDE-3 HDE-4 HDE-5 HDE-6 HDE-7 HDE-8 HDE-9 HDE-9 HDE-10

12 14 15 20 25 30 40 45 50 50 80

3400 3400 3400 3400 3400 3400 3400 3400 3400 5950 3400

2.71 2.05 1.86 2.04 2.10 2.4 3.09 3.27 4.25 1.94 5.52

Scaling factor, F

Crossover rate, Cr

Generations Best misfit

0.5 0.5

0.7 0.5

60 60

1.45 1.95

constant (3400). That means a case with a larger population size will get less generations and vice versa for the reason that the total number of the simulations is the product of population size and number of generations. Table 2 shows the values of the final best misfit in 10 different cases. From the results, it can be seen that the value of best misfit first decreases and then increases with increase in population size. The minimum misfit obtained is 1.86 with a population size of 15. 5.1.2. Performance comparison between HDE and DE in history matching Hajizadeh et al. (2010) applied 2 variants (DE rand and DE best) of differential evolution for history matching of PUNQ-S3. He conducted studies on control parameters. Table 3 shows the best results obtained for the two differential evolution variants with optimized control parameters. The DE selection strategy used in this paper is random selection. The parameter tuning is one of the drawbacks of every population based algorithm (Yang and He, 2013). It takes several sets of simulations to find which control parameters fits for the problem. The HDE gives a solution by its adaptive mechanism. Hajizadeh et al. (2010) got a best value of 1.95 for the DE rand variant after conducting several sets of simulation runs for parameter tuning with 3000 simulations each meanwhile HDE algorithm gives a best value of 1.86 with a total number of 3400 simulations. DE got the best value by using a population size of 50 whereas HDE used a population size of 15 for getting the best misfit value. When the comparison is done on the basis of number of iterations, HDE needs more iterations than DE to reach the best misfit value. This will become a limitation for HDE when using parallel computing where number of iterations is more restrictive than number of simulations. HDE with the population size of 50 took 5950 simulations to reach comparable misfit values with DE which is having optimized control parameters (Table 2). That means, with a population size of 50, an additional 59 iterations is required by HDE to get identical performance with parameter tuned DE algorithm which is having optimized control parameters. Since every problem is different, the choice of control parameters will be confusing for the user, and hence it is better to have an adaptive mechanism for control parameters than going for inappropriate constant control parameters which will give poor performance to the algorithm. Using HDE, user doesn't have to bother about the control parameters in the expense of some additional iterations or reduced population size. 5.1.3. Influence of population size in sampling the search space For further analysis regarding the performance of the algorithm, we took 3 out of the ten cases. The 3 cases chosen are HDE3, HDE-5 and HDE-7 (see Table 2 for nomenclature). The sampling histories of the algorithm in all generations for HDE-3, HDE-5, and HDE-7 are given in Fig. 4a–c respectively. These figures show the value of 45 unknown parameters corresponding to the best misfit model in all the generations. The values shown in these figures are scaled to [0, 1]. Among the three, the HDE-7 shows the maximum exploration of the search space but have a comparatively high misfit value (Table 2). HDE-7 spent a lot of simulations in searching good models throughout the entire search space but

E.C. Santhosh, J.S. Sangwai / Journal of Petroleum Science and Engineering 142 (2016) 21–35

Fig. 4. Sampling history of (a) HDE-3; (b) HDE-5; (c) HDE-7. y-Axis represent unknown parameter values scaled to [0,1] and x-axis represent generations.

27

28

E.C. Santhosh, J.S. Sangwai / Journal of Petroleum Science and Engineering 142 (2016) 21–35

Fig. 5. Best misfit value in 136 generations for HDE 5.

Fig. 6. Total field oil production during history match period (observed vs. simulated).

Fig. 7. Total field water production during history match period (observed vs. simulated).

E.C. Santhosh, J.S. Sangwai / Journal of Petroleum Science and Engineering 142 (2016) 21–35

29

Fig. 8. Well level history matching results of well bottom hole pressure (WBHP) and well water cut (WWCT) for the three wells (PRO-5, PR0-11, and PR0-12). Black lines shows the initial models, red line shows the best history matched model and blue dots shows the history matching data. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

didn't give much importance in refining the solution which reduced its exploitation quality. On the other side, HDE-3 has the lowest misfit value among the 3 cases (Table 2), but shows a poor exploration of the search space. HDE-3 spent a lot of simulations in refining the solution, thus giving less importance to the exploration of the search space. HDE-3 and HDE-7 shows the best and least convergence among the three cases respectively, but a lot of parameters hit the constraints in the case of HDE-3, which is not in the case of HDE-7. Considering the factors such as misfit, convergence of parameters, exploration of the search space, exploitation of the search space, HDE-5 shows a good trade-off. So, HDE-5 case is taken for further study.

5.1.4. History matching results for PUNQ-S3 In the HDE-5 case, iterations were continued for 136 generations with a population size of 25. The rationale behind the choice of 136 generation (¼3400/25) is based on the population size (25 for HDE-5 case, see Table 2) and number of simulations (3400). The best misfit value obtained is 2.1. Fig. 5 shows the value of best misfit function in each generation for the 136 generations. Figs. 6 and 7 show the match between the observed field production data and simulated production data of total well oil production and total well water cut, respectively, for the history matching period. Fig. 8 shows the well level history matching results of three wells (PRO-5, PRO-11, and PRO-12). The results show a good match which indicates the efficacy of HDE in history matching. Fig. 9a–c shows the porosity distribution,

30

E.C. Santhosh, J.S. Sangwai / Journal of Petroleum Science and Engineering 142 (2016) 21–35

Fig. 9. Distribution of important variable for maximum likelihood function based on history match (a) porosity distribution; (b) horizontal permeability distribution; (c) vertical permeability distribution.

E.C. Santhosh, J.S. Sangwai / Journal of Petroleum Science and Engineering 142 (2016) 21–35

31

Fig. 9. (continued)

Fig. 10. Histogram of the misfit distributions of the models generated from HDE algorithm and resampled by NAB algorithm.

horizontal permeability distribution and vertical permeability distribution in all the 5 layers of PUNQ-S3 corresponding to the maximum likelihood model obtained from HDE-5 after history matching. This completes the history matching part. But the risk associated with the prediction using a single history matched model still remains. In order to solve this problem, we used the information from all the models to predict the reservoir forecast using NAB algorithm. 5.2. Uncertainty quantification of PUNQ-S3 In this appraisal step, the ensemble of models (all the 3400) generated from the history matching step are submitted to the

NAB algorithm. The NAB algorithm builds a high dimensional parameter space using Voronoi cells with the ensemble of models. In order to calculate the posterior probability distribution of the models, a 15,000 step McMC walk is carried out in the search space using Gibbs sampler. The algorithm was modified to monitor the frequency of visits to each Voronoi cell during the McMC walk. The PPD of the ith model is calculated from the equation.

Pi =

fi n

∑ j =m1 f j

(8)

32

E.C. Santhosh, J.S. Sangwai / Journal of Petroleum Science and Engineering 142 (2016) 21–35

Fig. 11. Cumulative probability density function of FOPT after 16.5 years for HDE-3, HDE-5, and HDE-7.

where, nm is the number of models, f shows the frequency of visit in the Voronoi cell. j runs from 1 to number of models. 5.2.1. Resampling with NAB algorithm In HDE-5 case, the NAB algorithm resampled 1220 unique models from the ensemble of models generated during the history matching loop. Fig. 10 shows the misfit value distribution of both models sampled by hybrid differential evolution and resampled by NAB algorithm. The NAB selected unique models with low misfit value. Cumulative probability distribution function (CDF) is calculated from the posterior probability distribution of the resampled models. Fig.11 shows the cumulative probability distribution of field oil production total (FOPT) after 16.5 years for HDE-3, HDE-5, HDE-7, respectively (analyzed in detail subsequently in Section 5.2.2). From this, the Bayesian credible intervals (p10–p50–p90) for FOPT prediction after 16.5 years are calculated. 5.2.2. Influence of population size on uncertainty quantification Fig. 12 shows the comparison of the FOPT prediction with uncertainty estimates after 16.5 years between the three HDE cases along with the true case. It shows the influence of the population size of HDE (see Table 2 for details on population size) on the FOPT prediction with the uncertainty estimate. Even though HDE-3 gives FOPT prediction with less uncertainty, it failed to predict the

actual field oil production meanwhile the p50 value of the HDE-7 prediction comes very close to the truth case but the range of uncertainty is very large. The reason behind this is that, HDE-3 gives a good exploitation of the search space but very poor exploration (see Fig. 4a). HDE-7 gives good exploration of the search space but very poor exploitation (see Fig. 4c). From these observations, we can infer that a good exploration of search space gives good prediction, however without a proper exploitation of the search space, the prediction will contain large uncertainty range. In real life situations, each simulation will take much time, so there is not enough luxury to run the simulations until we get proper exploration and exploitation of the search space. In order to get a good prediction with minimum uncertainty range, a good trade-off is required between the exploration and exploitation within limited number of simulations. HDE-5 shows a good tradeoff between the exploration and exploitation of the search space. From this, it can be inferred that the population size which gives the best misfit need not give a good prediction. Exploration of the search space is as much important as finding a best misfit model. HDE-5 has this quality, so this case is taken for further studies. 5.2.3. Comparison between single and multiple history matched model prediction The cumulative probability distribution functions of FOPT for all the time steps for HDE-5 are calculated and from this, the Bayesian credible intervals (p10-p50-p90) for all time steps are determined. Fig. 13 shows the truth case, FOPT prediction by maximum likelihood model and Bayesian credible intervals for all time steps. The red line shows the prediction of the most likelihood model obtained from HDE-5 case. The history matching phase is up to 8 years and the prediction is done up to 16.5 years. It can be seen that the most likelihood model shows deviation from the truth case and the maximum deviation is shown near 4000 days. Even at the time of maximum deviation of the maximum likelihood model from the truth case, the P10–P90 confidence bound comfortably holds the truth case in all the time steps. So the predicted approach with multiple history matched models shows a promising future in field management. 5.2.4. Comparison with other approaches Different authors have applied different approaches to predict the FOPT after 16.5 years with uncertainty estimates. Fig. 14 shows a comparison of these approaches with this work of the same PUNQ-S3 model (Floris et al., 2001; Demyanov et al., 2004; Holmes et al., 2007; Hajizadeh et al., 2010). With this study, it can be

Fig. 12. FOPT prediction after 16.5 years for the case HDE-3, HDE-5, HDE-7.

E.C. Santhosh, J.S. Sangwai / Journal of Petroleum Science and Engineering 142 (2016) 21–35

33

Fig. 13. FOPT prediction (P50) bounded by p10/p90 vs. truth case and maximum likelihood model.

Fig. 14. Uncertainty estimates calculated by HDE-NAB framework and comparison with other approaches. Other from various literatures (Floris et al., 2001; Demyanov et al., 2004; Holmes et al., 2007; Hajizadeh et al., 2010); HDE-NAB (this work).

inferred that the approaches which show large deviation of FOPT prediction with the truth case with small uncertainty range may lack proper exploration of the search space. The approaches which show large uncertainty ranges may have explored the search space without proper exploitation. The Fig. 14 shows that HDE gives good FOPT prediction with minimum uncertainty range among the approaches which predicted FOPT close to the true value. HDE gave better prediction results than the parameter tuned DE variants (Hajizadeh et al., 2010) without many series of simulations for parameter tuning. The test was repeated three times in which the approach displayed identical results which makes it robust. It is also observed that the approach used in this work is comfortable and reliable in the prediction of the production forecast as compared to other methodologies used in the open literature. 6. Conclusion The process of history matching plays an important role in the reservoir model development. Since the simulations are time

consuming, the approaches taken for history matching should use limited number of simulations. This can be achieved by efficient optimization algorithms which can navigate through high dimensional search space and find multiple history matched models. Choosing the optimal control parameters was one of the problems faced by population based algorithms and this difficulty is solved by using an adaptive mechanism for control parameters. In this work, the influence of population size of hybrid differential algorithm in history matching is studied where the value of best misfit first decreased and then increased with increase in the population size. When compared with parameter tuned DE algorithm, the HDE algorithm was able to achieve similar misfit values of that of parameter tuned DE algorithm without conducting several sets of simulation runs for parameter tuning. When the comparison is done on the basis of number of iterations, HDE needs more iterations than DE to reach the best misfit value. This will become a limitation for HDE when using parallel computing where number of iterations is more restrictive than number of simulations. HDE gave equally good results even when the search started from

34

E.C. Santhosh, J.S. Sangwai / Journal of Petroleum Science and Engineering 142 (2016) 21–35

different models. The history matching of PUNQ-S3 using HDE showed good results. So, HDE can be used as a promising tool in history matching. Due to the non-uniqueness of the inverse problems, the decision making for field development using a single reservoir model carries more risk. Forecasting within confidence interval, thus quantifying the uncertainty showed a promising way to solve this problem. The influence of population size on uncertainty estimates was carried out. Results showed that as the population size decreases the uncertainty range also decreases, but the algorithm fails to predict the true FOPT because of the poor exploration of the search space. On the other hand, as the population size increases the algorithm is able to make a good prediction but the uncertainty range in the prediction also increases. The reason behind this is because of the strategy adopted in the NAB algorithm is constant misfit in the Voronoi cells. The accuracy of the interpolation depends upon how well the optimization algorithm samples the high dimensional model space. The efficiency of the results will depend upon the structure of model space created by the NAB algorithm using the ensemble of models generated during history matching phase. A proper sampling of the search space is needed in the generation of models while history matching to build a good model space in NAB. So a good trade-off between the exploration and exploitation of the search space by the optimization algorithm is required to make a good prediction. Thus, the accuracy in the production forecast with the uncertainty estimate not only depends on the NAB algorithm but also on the optimization algorithm. It can also be concluded that the population size which gives the best misfit value need not give a good prediction. The HDE-NAB approach gave a good mean FOPT prediction with decent uncertainty range. Hence, HDE-NAB framework proves to be an efficient tool for history matching and quantifying the uncertainty.

Acknowledgments We thank Indian Institute of Technology, Madras for providing facilities to conduct this research and Schlumberger for providing Eclipse reservoir simulation software. We also thank Prof. Malcolm Sambridge of the Australian National University; Canberra, Australia for providing a basic NAB routine. We show appreciation for the valuable comments from the experts in this field especially to Dr. Nigel Goodwin (Essence Products and Services Ltd., UK) and Prof. Mike Christie (Heriot-Watt University, Edinburg, UK). We sincerely thank three anonymous reviewers and Editor-in-Chief Dr. Vural Sander Suicmez for rigorous review and suggestions to improve the work. Authors would also like to acknowledge the receipt of ‘Prof. C.S. Krishnamoorthy Endowment Prize’ from IIT Madras for the best project in the area of genetic algorithms and evolutionary computation for the year 2015.

References Boggs, P., Tolle, J., 2000. Sequential quadratic programming for large scale nonlinear optimization. J. Comput. Appl. Math. 124 (1–2), 123–137. http://dx.doi.org/ 10.1016/S0377-0427(00)00429-5. Boss, C., 1999. Production Forecasting with Uncertainty Quantification. Netherlands Institute of Applied Geoscience TNO, Netherlands. Castellini, A., Gullapalli, I., Hoang, V., Condon, P., 2005. Quantifying uncertainty in production forecast for fields with significant history: a West African case study, IPTC 10987. International Petroleum Technology Conference, Doha, Qatar. 21–23 November 2005. Das, D., Suganthan, P.N., 2010. Differential evolution: a survey of the state-of-theart. Evol. Comput., IEEE Trans. 15, 4–31.

Demyanov, V., Subbey S., Christie, M.A., 2004. Uncertainty Assessment in PUNQS3– neighbourhood algorithm framework for geostatistical modelling. In: Proceedings of the 66th EAGE Conference & Exhibition-Workshops. Dogru, Ali H., John, H., 1981. Comparison of sensitivity coefficient calculation methods in automatic history matching (includes associated paper 10873). Soc. Pet. Eng. J. 21 (05), 551–557. Epitropakis, Michael G., Plagianakos, Vassilis P., Vrahatis, Michael N., 2008. Balancing the exploration and exploitation capabilities of the differential evolution algorithm. Evolutionary Computation, CEC. (IEEE World Congress on Computational Intelligence). IEEE Congress on. IEEE, 2008. Erbas, D., Christie, M., 2007. Effect of sampling strategies in prediction uncertainty estimation, SPE 106229. SPE Reservoir Simulation Symposium, Houston, Texas, U.S.A. 26–28 February. Fletcher, R., 1987. Practical Methods of Optimization, 2nd ed. Wiley, New York. Floris, F.J.T., Bush, M.D., Cuypers, M., Roggero, F., Syversveen, A.-R., 2001. Methods for quantifying the uncertainty of production forecasts. Pet. Geosci. 7, 87–96. Gelman, A., Carlin, J.B., Stern, H.S., Rubin, D.B., 1995. Bayesian Data Analysis. Chapman & Hall, London. Geman, S., Geman, D., 1984. Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images. IEEE Trans. Pattern Anal. Mach. Int. 6, 721–741. Hajizadeh, Y., Michael, A.C., Demyanov, V., 2011. Ant colony optimization for history matching and uncertainty quantification of reservoir models. J. Pet. Sci. Eng. 77, 78–92. Hajizadeh, Y., Michael, A.C., Demyanov, V., 2010. Comparative study of novel population-based optimization algorithms for history matching and uncertainty quantification: PUNQ-S3 revisited. Society of Petroleum Engineers-14th Abu Dhabi International Petroleum Exhibition and Conference, ADIPEC-2010, SPE136861-MS. Holmes, J.C., McVay, D.A., Senel, O., 2007. A system for continuous reservoir simulation model updating and forecasting, SPE 107566. Digital Energy Conference and Exhibition, Houston, Texas, U.S.A. 11–12 April. Jacquard, P., 1965. Permeability distribution from field pressure data. Soc. Pet. Eng. J. 5, 281–294. LaTorre, A., Muelas, S., Peña, J.-M., 2010. A mos-based dynamic memetic differential evolution algorithm for continuous optimization: a scalability test. Soft Comput. – Fusion Found. Methodol. Appl. 15, 1–13. http://dx.doi.org/10.1007/ s00500-010-0646-3. Maucec, M., Douma, S., Hohl, D., Leguijt, J., Jimenez, E.A., Gupta, A.D., 2007. Streamlinebased history matching and uncertainty, Markov-Chain Monte Carlo study of an offshore oil field, SPE 109943. In: Proceedings of the SPE Annual Technical Conference and Exhibition. Anaheim, California, U.S.A. 11–14 November. Mohamed, L., Christie, M., Demyanov, V., 2009. Comparison of stochastic sampling algorithms for uncertainty quantification, SPE 119139. In: Proceedings of the SPE Reservoir Simulation Symposium, Woodlands, Texas, U.S.A. 2–4 February. Neri, F., Tirronen, V., 2009. Scale factor local search in differential evolution. Memet. Comput. 1, 153–171. Noman, N., Iba, H., 2008. Accelerating differential evolution using an adaptive local search. Evol. Comput., IEEE Trans. 12 (1), 107–125. Price, K.V., 1999. An Introduction to Differential Evolution. McGraw-Hill Ltd., UK, Maidenhead, UK, England, pp. 79–108. PUNQ-S3, 2015. Department of Earth Science and Engineering. Imperial College London 〈https://www.imperial.ac.uk/engineering/departments/earth-science/ research/research-groups/perm/standard-models/〉 (last accessed 15.05.15.). Reklaitis, G.V., Ravindran, A., Ragsdell, K.M., 1983. Engineering Optimization: Methods and Applications. Wiley-Interscience, United States. Reynoso-Meza. G., Sanchis, J., Blasco, X., Herrero M.J., 2011. Hybrid DE algorithm with adaptive crossover operator for solving real-world numerical optimization problems. In: Proceedings of the IEEE Congress on Evolutionary Computation (CEC-2011), New Orleans, LA, 5–8th June. Sambridge, M., 1999a. Geophysical inversion with a Neighbourhood Algorithm — I searching a parameter space. Geophys. J. Int. 138, 479–494. Sambridge, M., 1999b. Geophysical inversion with a Neighbourhood Algorithm — II appraising the ensemble. Geophys. J. Int. 138, 727–745. Sangwai, J.S., Bhat, S.A., Saraf, D.N., Gupta, S.K., 2007. An experimental study on online optimizing control of free radical bulk polymerization in a rheometer-reactor assembly under conditions of power failure. Chem. Eng. Sci. 62, 2790–2802. Schlumberger, Eclipse reservoir simulator, 2010. Sheldon, J.W., Harris, C.D., Bavly, D., 1960. A method for general reservoir behavior simulation on digital computers. Fall Meeting of the Society of Petroleum Engineers of AIME, 2–5 October, Denver, Colorado. Smith, A.F.M., Roberts, G.O., 1993. Bayesian computation via the Gibbs sampler and related Markov chain Monte Carlo methods. J. R. Stat. Soc., B 55, 3–23. Sousa, S.H.G., Maschio, C., Schiozer, D.J., 2006. Scatter search metaheuristic applied to the history matching problem, SPE 102975. In; Proceedings of SPE Annual Technical Conference and Exhibition, San Antonio, Texas, U.S.A., 24–27 September. Storn, R., 2008. Differential Evolution Research–Trends and Open Questions. Advances in Differential Evolution. Springer, Berlin, Heidelberg, pp. 1–31. Subbey, S., Christie, M., Sambridge, M., 2003. A strategy for rapid quantification of uncertainty in reservoir performance prediction, SPE 79678. In: Proceedings of SPE Reservoir Simulation Symposium, Houston, Texas, U.S.A., 3–5 February.

E.C. Santhosh, J.S. Sangwai / Journal of Petroleum Science and Engineering 142 (2016) 21–35

Sultan, A.J., Ouenes, A., Weiss, W.W., 1994. Automatic history matching for an integrated reservoir description and improving oil recovery, SPE 27712. In: Proceedings of Permian Basin Oil and Gas Recovery Conference, Midland, Texas, U. S.A., 16–18 March. Thomas, K.L., Hellums, L.J., Reheis, G.M., 1972. A nonlinear automatic history matching technique for reservoir simulation models. Soc. Pet. Eng. J. 12, 6. Wang, J., Buckley, J.S., 2006. Automatic history matching using differential evolution algorithm. In: Proceedings of International Symposium of the Society of Core Analysts, Trondheim, Norway, 12–16 September, 2006. Yang, Xin-She, He, Xingshi, 2013. Bat algorithm: literature review and applications. Int. J. Bio-Inspired Comput. 5 (3), 141–149. Zaharie, D., 2001. On the explorative power of differential evolution algorithms, In:

35

Proceedings of the 3rd International Workshop. Symbolic and Numeric Algorithms of Scientific Computing, SYNASC, Romania. Zaharie, D., 2002. Critical values for the control parameters of differential evolution algorithms. In: Proceedings of the 8th International Conference of Soft Computing, pp. 62–67. Zaharie, D., 2003. Control of population diversity and adaptation in differential evolution algorithms. In: Matouek, R., Omera, P. (Eds.), Proceedings of Mendel, Ninth International Conference on Soft Computing, pp. 41–46.