Ant colony optimization for history matching and uncertainty quantification of reservoir models

Journal of Petroleum Science and Engineering 77 (2011) 78–92

Contents lists available at ScienceDirect

Journal of Petroleum Science and Engineering j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / p e t r o l

Ant colony optimization for history matching and uncertainty quantification of reservoir models Yasin Hajizadeh ⁎, Mike Christie, Vasily Demyanov Institute of Petroleum Engineering, Heriot Watt University, Edinburgh, EH14 4AS, UK

a r t i c l e

i n f o

Article history: Received 18 May 2010 Accepted 15 February 2011 Available online 23 February 2011 Keywords: History matching Uncertainty quantification Ant colony optimization Bayesian inference Reservoir simulation

a b s t r a c t This paper introduces a new stochastic approach for assisted history matching based on a continuous ant colony optimization algorithm. Ant Colony Optimization (ACO) is a multi-agent optimization algorithm inspired by the behavior of real ants. ACO is able to solve difficult optimization problems in both discrete and continuous variables. In the ACO algorithm, each artificial ant in the colony searches for good models in different regions of parameter space and shares information about the quality of the models with other agents. This gradually guides the colony towards models that match the desired behavior — in our case the production history of the reservoir. The use of ACO for history matching has been illustrated on two reservoir simulation cases. The first case is Teal South model which is a real reservoir with a simple structure and a single producing well. History matching of this model is a low dimensional problem with eight parameters. The second case study is PUNQ-S3 reservoir which has a more complex geological structure than Teal South model. This problem entails solving a high dimensional optimization problem. © 2011 Elsevier B.V. All rights reserved.

1. Introduction The development of efficient history matching and uncertainty quantification techniques for use in prediction of reservoir performance is an important and challenging topic in petroleum engineering. History matching is a very complex non-linear and ill-posed inverse problem in which we aim to calibrate a reservoir model to reproduce historical observation data such as production rates or pressure measurements in a reservoir. Like most inverse problems, history matching has non-unique solutions, which means that different combinations of parameters can produce good matches of the observed data. It is necessary to obtain multiple good history-matched models in order to realistically quantify uncertainty in reservoir simulation models (Erbas and Christie, 2007). Obtaining multiple history matched models requires assisted history matching methods which can handle realistic full field models containing lots of production data and potentially a large number of unknown parameters. There are numerous examples of assisted history matching techniques in the literature. One of the first attempts to present an assisted framework for history matching was by Chen et al. (1974) in which the history matching was formulated as an optimal control problem. Today, various broad categories of history matching algorithms exist, namely gradient methods (Anterion and Eymard,

⁎ Corresponding author. Tel.: +44 131 4513563; fax: +44 1314513127. E-mail address: [email protected] (Y. Hajizadeh). 0920-4105/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.petrol.2011.02.005

1989), stochastic samplers, and particle filter methods such as Ensemble Kalman Filter (EnKF) (Liu and Oliver, 2005; Jafarpour and McLaughlin, 2007; Lodoen and Omre, 2008). Stochastic methods are easy to implement and they have been very popular in recent years with a number of stochastic algorithms successfully applied to history matching. Examples of recent applications of stochastic samplers include: Neighbourhood Algorithm (NA) (Subbey et al., 2003), Genetic Algorithms (GA) (Castellini et al., 2005), Simulated Annealing (SA) (Sultan et al., 1994), Scatter Search (SS) (Sousa et al., 2006), Tabu Search (TS) (Yang et al., 2007), Hamiltonian Monte Carlo (HMC) (Mohamed et al., 2009), Particle Swarm Optimization (PSO) (Mohamed et al., 2009; Kathrada 2009), Markov Chain Monte Carlo (MCMC) (Maucec et al., 2007), Simultaneous Perturbation Stochastic Approximation (SPSA) (Gao et al., 2007) and Chaotic Optimization (CO) (Mantica et al., 2002). Any optimization algorithm used for assisted history matching must be fast in navigating high dimensional search spaces and efficient in finding multiple good models with a limited number of simulations. The importance of the optimization algorithm choice and tuning becomes essential when we try to estimate large number of parameters in the presence of multiple local minima. Finding multiple local minima are essential for realistic quantification of prediction uncertainty, which is often not the case in many traditional optimization problems. The sampling performance of algorithm has a direct impact on the uncertainty estimate of future production (Erbas, 2006). Recent developments in swarm-based algorithms provide a more flexible optimization tool by interaction of simple agents. These algorithms offer advance tuning options to maintain a

Y. Hajizadeh et al. / Journal of Petroleum Science and Engineering 77 (2011) 78–92

79

Fig. 1. Solution archive in ACOR and its components. Adopted from (Socha and Dorigo, 2008).

good balance between exploration, exploitation and convergence speed while searching for optimal solutions. In this paper a new stochastic search method will be introduced for generating multiple history matching models. This method is called Ant Colony Optimization (ACO) and is one of the most successful proposed metaheuristics for optimization problems based on swarm intelligence (Dorigo et al., 1999; Bonabeau et al., 2000). The swarm intelligence (SI) paradigm originated from studying the colonies or social organisms which consist of single and simple interacting agents which can accomplish difficult tasks when forming a group. One of the major advantages of these optimization methods is their simplicity and efficiency in navigating high dimensional search spaces. In the following sections we will describe the fundamental ideas behind this method and show the application of ant colony optimization for history matching of two reservoir simulation models.

2. Ant Colony Optimization (ACO) The Ant Colony Optimization (ACO) metaheuristic was introduced by Dorigo (1992). It is inspired by studying the behavior of real ants searching for their food source. In ACO, a set of artificial ants cooperate in finding good solutions by indirect exchange of information via artificial pheromone. The main underlying idea for ACO is the parallel

search over several computational threads based on local problem data. This search is based on a dynamic memory structure of the characteristics of the problem and the quality of obtained results which are built up over time. In ACO, pheromone is typically associated with the solution components used by artificial ants to construct new solutions, guiding their decisions. This artificial pheromone is accumulated at run-time through a learning mechanism that rewards good solutions. More pheromone on a path increases the probability of that path being followed. An evaporation rule will be tied with the pheromones, which will reduce the chance for poor quality solutions. Through indirect communication with other ants via foraging behavior, a colony of searching agents can establish the optimal solution to the problem, over time with a positive feedback loop known as stygmergy. Any ACO approach must define the following parts for a new problem: (1) A heuristic function which guides the search with problem specific information. (2) A trail definition which states what information is to be stored and allows ants to share information about good solutions. (3) An update rule that defines the way in which good solutions are reinforced. (4) A fitness function which determines the quality of a particular solution.

Fig. 2. Gaussian distributions in ACOR which are sampled to obtain new models.

80

Y. Hajizadeh et al. / Journal of Petroleum Science and Engineering 77 (2011) 78–92 Table 2 Summary of parameters used in ACOR for the best history matching result.

Fig. 3. Teal South reservoir in the Gulf of Mexico.

The first ant colony optimization (Dorigo, 1992) was designed for discrete optimization problems. One of the most successful attempts to extend original ACO concept to handle continuous variables was by Socha and Dorigo (2008) within the framework of ACOR algorithm. It has been shown by the authors that ACOR algorithm outperforms other variants of ACO developed for continuous optimization.

Number of ants (m)

Archive size (k)

Search localization (q)

Pheromone evaporation (ξ)

Simulations

50

50

0.001

0.5

1000

where k is the size of the archive and q is a parameter of the algorithm that controls the balance between exploration/exploitation in ACOR. Small values of q give preference to the best-ranked solutions, while for larger q values the probability of selecting models in the archive becomes more uniform. In ACOR the selection probability is not a direct function of the fitness, but is computed indirectly from the fitness considering the effect of parameter q. In order to generate new solutions, one of the solutions in the archive is chosen according to the following probability: pl =

ð2Þ

∑ ωr

2.1. ACOR algorithm

r=1

2

ðl−1Þ 1 − pffiffiffiffiffiffi e 2q2 k2 qk 2π

ð1Þ

To construct a single model, at each step (i = 1, 2 … n) we sample a value for each unknown parameter in the problem. The mixture of Gaussian kernels form the probability density function (PDF) which will be sampled to obtain new members of the archive. For example in Fig. 2, we have 4 models in the archive to compute this mixture. The components of the solutions are used to dynamically generate probability density functions and modify their shape. For each dimension of the problem, there are k individual Gaussian functions. Each Gaussian distribution is characterized by its mean and standard deviation values. Mean value of each distribution is equal to the corresponding value in the solution archive. Considering a single column in the archive, the first number in this column will be the mean value for the first individual Gaussian function and so on. For the standard deviations, we compute the average distance between the selected member (sl) and the other members of that column in the solution archive. This computed value is multiplied by ξ which is called pheromone evaporation rate. The higher ξ value, worst solutions are forgotten faster. With lower ξ values, the search is less concentrated on the previously visited areas of the search space; hence the convergence speed of the algorithm will be higher. For

History match for best model 2500 2000

Oil rate (STB/D)

In ACOR we have a solution archive with k models in which we keep track of solutions (Fig. 1). The solutions in this archive are ranked based on their quality, which means models with lower misfit values get a higher rank and will appear at the top of the list. There are k rows and n columns in the archive, where k is the number of models that are kept in the archive and n is the number of dimensions of the problem or parameters in history matching. Each single row in this archive consists of a vector of reservoir model parameters (s) and corresponding objective function f(s) value that we obtain after running the flow simulation. The ith unknown parameter value of the lth model is denoted by sil. The misfit values show how well the proposed model can fit historical data. The last column in this archive includes the weights (ωl) of the solutions. These weights are a function of solution quality (Eq. (1)) and will be used to probabilistically build new solutions. The archive is initially filled with random solutions and the misfits (f(s)) of models are evaluated. If the number of ants evaluated at each iteration of the ACOR algorithm is m, then at each iteration, m new solutions are added to the population and from the archive which now contains k + m models, the m worst solutions are removed to keep the archive size fixed. This action simulates the pheromone update part in discrete ACO. The remaining models in the archive are sorted according to their misfit score. The aim is to bias the search process towards the good regions of the search space with low misfit models by probabilistically constructing new solutions. Next, we compute the weights for each model in the archive. The weight of each member in the archive is computed based on the following equation and will be used to probabilistically select the members of archive:

ωl =

ωl k

1500 1000 500

Table 1 Parameterization and prior ranges for Teal South model.

0 0

Parameters

Units

Prior range

kh (for each layer) kv/kh Rock compressibility Aquifer strength

mD – psi− 1 MMSTB

10{1, 3} 10{− 4, − 1} 5 × 10− 6–1 × 10− 4 10{7, 9}

200

400

600

800

1000

1200

Time (days) Observed

Simulated

Fig. 4. Simulated and observed values for the best history matched model in Teal South reservoir.

Y. Hajizadeh et al. / Journal of Petroleum Science and Engineering 77 (2011) 78–92 Table 3 Parameters used in ACOR for 10 random runs.

Table 4 Parameters used in ACOR for testing the effect of q with ξ = 0.5.

Number of ants (m)

Archive size (k)

Search localization (q)

Pheromone evaporation (ξ)

Simulations

10

50

0.01

0.5

500

obtaining new solutions, ants sample the mixture of Gaussian functions. The mixture of Gaussian kernels (Fig. 2, bold line) is the weighted sum of the individual Gaussian kernels and is computed with the following equation: ðs−μli Þ2 − 1 i G ðsÞ = ∑ ωl i pffiffiffiffiffiffi e 2σl 2 σl 2π l=1 i

81

k

ð3Þ

where μ li and σ li are the mean and standard deviations of the Gaussian distributions. In the next step the Gaussian kernel is sampled to obtain a new model and this process is continued for all dimensions of the problem by each of the m ants, and at each iteration of the algorithm. The new members replace the models with least performance and this process continues until the search terminates with a stopping condition. 3. Application of ACOR in history matching To test the application of ant colony optimization (ACOR) algorithm for history matching, we consider two reservoir models. Teal South case study is a real reservoir with a simple structure and one well. This model is used to understand the behavior of ACOR in a low dimensional history matching problem. The second case that is used in this study – PUNQ-S3 reservoir – is a synthetic problem which is widely used to benchmark the optimization algorithms for history matching. PUNQ-S3 model has a more complex geological structure than Teal South and history matching is based on multivariate data and multiple wells. 3.1. Teal South reservoir description Teal South (Fig. 3) is located in Eugene Island Block 354 in the Gulf of Mexico. It was originally operated by Texaco, later by Apache and these days by Devon Energy. Gas oil ratio of the field is exceeding 800. There are multiple sands in the Teal South reservoir and we are looking at the 4500 ft sand in this study. A single well penetrates the 4500 ft sand which was initially overpressured at 3096 psi.

Number of ants (m)

Archive size (k)

Search localization (q)

Simulations

50

50

0.01 and 0.5

1000

The reservoir simulation model is set up on an 11 × 11 × 5 corner point grid. There are five geological layers in this model with uniform properties. Key unknown parameters in history matching are horizontal permeability multipliers for each of the five layers (P1–P5), a single value for vertical to horizontal permeability ratio (P6), rock compressibility (P7) and aquifer strength (P8). Parameterization for the Teal South model and their prior range are shown in Table 1. At this stage, we chose to use all available production history for matching. History matching is only done on field oil production rate, which is included in a univariate objective function (Eq. (4)) ðqobs ðti Þ−qsim ðti ÞÞ2 2σ 2 n=1 N

M= ∑

ð4Þ

where M is the misfit value, N is the number of observations, q is the flow rate for observed and simulated data, and σ2 is the variance of the observed data. Teal South reservoir has been used in previous studies of stochastic sampling algorithms: Christie et al. (2002) used Neighbourhood Algorithm (NA) for finding multiple history matched models for the Teal South reservoir. Erbas and Christie (2007) used a genetic algorithm for this purpose and compared the results with NA and recently Mohamed et al. (2009) compared the performance of NA, Particle Swam Optimization (PSO) and Hamiltonian Monte Carlo (HMC) for finding history matched models in Teal South reservoir. In following sections we report the results of the ant colony optimization application to Teal South reservoir. ACOR algorithm is coupled with Eclipse reservoir simulator. Additional routines are used to prepare the input data for simulations based on the results of optimization and process the simulation output files for computing the misfit values. These misfit values are given to the ACOR algorithm to propose the next set of solutions. We use a random starting population in all of our simulations. However in real field applications there might be some information already available from the previous studies about the possible solutions for history matching. This information then can be used to narrow down the initial range for the parameters that will be tuned using the optimization engine. 3.1.1. Best history matching in Teal South We have tried several combinations of tuning parameters in ACOR algorithm to obtain a good match for oil production rate in the Teal South reservoir Hajizadeh et al., (2009). The tuning parameters resulted in the best match are presented in Table 2. Fig. 4 shows the best history matched model in the ensemble which captures the observed oil production data points very well. We obtain the misfit value of 14.52 for the best member in the ensemble.

Table 5 Parameters used in ACOR for testing the effect of ξ with q = 0.5.

Fig. 5. Best Misfit obtained for 10 random runs in Teal South reservoir.

Number of ants (m)

Archive size (k)

Pheromone evaporation (ξ)

Simulations

50

50

0.1 and 0.9

1000

82

Y. Hajizadeh et al. / Journal of Petroleum Science and Engineering 77 (2011) 78–92

0.01

q

0.5

Minimum misfit at the end of optimization

15

16

17

18

19

20

Misfit Fig. 6. Minimum misfits obtained for two search localization parameters at a fixed pheromone evaporation rate.

3.1.2. ACOR sensitivity to initial run conditions Any stochastic sampling algorithm is sensitive to the initial random starting samples, and conclusions about efficiency of any algorithm should be based on an overage performance over a number of initial samples. To test the performance of the ACOR algorithm, we performed 10 runs of 500 simulations. These 10 runs all had the same tuning parameters, but had a different random set of initial samples. The parameters used in ACOR algorithm for this experiment are shown in Table 3. For Teal South reservoir models with misfit values under 20 match history well, with an average error of 1 standard deviation or less. As it can be seen from Fig. 5 which shows the best misfit values obtained for each run, there are no failed runs among these 10 trials, e.g. the models that are trapped in local minima with a misfit value greater than 20. All the runs converge to regions of the search space with misfits smaller than 20. Although the number of simulations are considerably less than the ones used for finding the best model (500 comparing to 1000), three of the models have the minimum misfits between 15 and 16 which are close to the best model found in the previous section with a misfit of 14.52. Four models have converged to regions with misfits between 16 and 18 and only three models have misfits between 19 and 20.

3.1.3. Effect of search localization (q) and pheromone evaporation rate (ξ) on the performance of ACOR In ACOR there are two parameters that control the behavior of the algorithm and its sampling performance. Search localization parameter (q) controls the preference of good models in ACOR and pheromone evaporation rate indicated how quickly ACOR forgets the bad models. To understand the effect of these parameters, we have performed two separate tests. In the first test we have two search localization values (q = 0.01 and q = 0.5) with a fixed pheromone evaporation rate (ξ = 0.5). The second test keeps the search localization parameter fixed (q = 0.5) and repeats the simulations for two different pheromone evaporation rates (ξ = 0.1 and ξ = 0.9). Tables 4 and 5 summarize the tuning parameters used in these test. Each test has been repeated for 5 times. Figs. 6–9 show the results of these two tests using boxplots. Boxplots are used to show the minimum, quartiles, median and maximum of results. Circles outside the box show outlier data which are 1.5× interquartile range (IQR). Fig. 6 shows the boxplots for the minimum misfit values obtained for two different search localization values. Examining this figure, we see that using higher q values, on average, we get a lower final misfit value. Fig. 7 shows the boxplots for the first simulation which we can get a misfit value less than 20. From Fig. 7

0.01

q

0.5

Simulations to get misfit values less than 20

0

100

200

300

400

500

600

Simulation Fig. 7. Number of simulation required to a get a misfit value less than 20 for two search localization parameters at a fixed pheromone evaporation rate.

Y. Hajizadeh et al. / Journal of Petroleum Science and Engineering 77 (2011) 78–92

83

0.1

xi

0.9

Minimum misfit at the end of optimization

15

16

17

18

19

20

Misfit Fig. 8. Minimum misfits obtained for two pheromone evaporation rates at a fixed search localization value.

we understand that for higher q values, it takes longer for the algorithm to find a model with a misfit value less than 20. The median of boxplot for q = 0.01 is 59 simulations, while this value for q = 0.5 is 136 simulations. Fig. 8 compares the boxplots for the minimum misfit values in the second test, in which we keep the search localization parameter fixed (q = 0.5) and we study the effect of pheromone evaporation rate on the performance of ACOR algorithm. We can see that a lower pheromone evaporation rate (ξ = 0.1) results in a better performance of algorithm and, on average, misfit values are lower than the case with high pheromone evaporation rate (ξ = 0.9). Fig. 9 shows that with lower ξ values, the algorithm finds low misfit values (M b 20) faster. When ξ = 0.1, the median of the boxplot is 55 simulations, while for ξ = 0.5 the median value is 224 simulations. With lower pheromone evaporation rates, bad models are kept for a longer time in the solution archive and thus the algorithm concentrates on the regions with better misfit values. 3.1.4. Comparison of ACOR and NA Neighbourhood Algorithm (NA) is a recent stochastic optimization method that is proposed by Malcolm Sambridge (Sambridge, 1999a). It was originally developed for solving inverse problems in seismology

and soon it was applied in other fields including petroleum engineering (Christie et al., 2002; Subbey et al., 2003; Okano et al., 2006). We compare the results of ACOR with NA in section because 1) NA, like ACOR, is a stochastic optimization method which uses a structure of population members to find the good fitting regions in the search space and 2) NA has been proven to be an efficient tool in finding multiple history matched models even in real life problems (Valjak, 2008). The Neighbourhood Algorithm (NA) starts the optimization with an initial random generation and after calculating the fitness of those models using the objective function, the algorithm determines the nearest neighbor region of each model in the parameters space by constructing the Voronoi diagram. Then it chooses nr best models according to their objective function value and generates ns new models by a uniform random walk in the Voronoi cells of the models. In order to compare NA and ACOR algorithms in terms of their misfit reduction efficiency, we have attempted to set the NA parameters as close as possible to ACOR run (these are presented in Table 6). In this test, the minimum misfit obtained for the ACOR run is 14.71. For NA case we have tried 4 different setups with different algorithm parameters. Each of these four cases has been repeated for ten times in order to consider the effect of random initialization of NA. In each of these

0.1

xi

0.9

Simulations to get misfit values less than 20

0

200

400

600

800

Simulation Fig. 9. Number of simulation required to a get a misfit value less than 20 for two pheromone evaporation rates at a fixed search localization value.

84

Y. Hajizadeh et al. / Journal of Petroleum Science and Engineering 77 (2011) 78–92

Table 6 Parameters used for two different ACOR runs. Case

m

k

q

ξ

Total simulations

Misfit

ACO-TS-1 ACO-TS-2

25 25

50 50

0.001 0.1

0.5 0.5

1000 1000

14.71 14.75

tests, there are 50 initial models and 25 cells will be resampled at each iteration. The optimization continues for 38 iterations, which makes the total number of models equal to 1000. The only difference between different NA setups is the nr value. Increasing nr will result in a wider exploration of parameter space. The effect of the parameter nr in NA algorithm can be considered same as the parameter q in the ACOR algorithm in which higher q values makes the algorithm more explorationary. Summary of the parameters used for NA runs in Teal South reservoir (NA-TS) and the best misfit values in each setup and among 10 trials are presented in Table 7. The minimum misfit for NA algorithm in these tests has been obtained for NA-TS-2 case (M = 14.83). We have chosen the NA-TS-2 setup to compare the results with ACOR run. Comparison of the best NA misfit (14.83) with the minimum misfit obtained by ACOR with a similar setup (14.71) indicates that for this particular problem ACOR obtains marginally better models than NA in terms of misfit value. Figs. 10 and 11 show the boxplot diagrams for ACOR and NA algorithms. Here we use boxplots to show the minimum, quartiles, median and maximum of the members in each generation. These figures show that ACOR reduces the misfit value faster than NA, for the specific algorithm parameters used here. Figs. 12–14 show the sampling history of two ACOR and one NA cases. Each figure has 8 tiles, showing 8 parameters of history matching. Horizontal axis is the number of simulations and vertical axis shows the scaled values of each parameter between 0 and 1. From these figures we can see that ACO-TS-2 is able to maintain larger population diversity due to its higher value for search localization (q). On the other hand, ACO-TS-1 has a faster convergence to good fitting regions than ACO-TS-2 due to its smaller q value; however this faster convergence comes at the expense of losing population diversity. Also we can see that NA-TS-2 can simultaneously perform search in two separate regions of the parameters as it is progressing toward the end of optimization cycle. Such performance is especially noticeable for parameters P2 and P4 in Fig. 14. 3.1.5. Observation data in history matching and its effect on the performance of ACOR In this part we address the importance of observation data used in history matching of reservoir models. We illustrate how additional observation points affect the performance of the ant colony optimization. Later in this paper we will also investigate the effect of additional information of the uncertainty of the predictions made by ACOR. We consider three scenarios where the number of observation data points used in the objective function is different. In each scenario, we use sampling history figures of algorithms to show their different

Table 7 Parameters used for 4 different NA runs in Teal South reservoir. Case

nsi

ns

nr

Iterations

Total simulations

Misfit

NA-TS-1 NA-TS-2 NA-TS-3 NA-TS-4

50 50 50 50

25 25 25 25

2 5 12 25

38 38 38 38

1000 1000 1000 1000

18.71 14.83 15.03 15.16

performance during the optimization process. In the first scenario for Teal South reservoir we have used 6 oil rate measurements during 181 days of production, the second scenario involves 23 data points observed in 699 days and the third scenario has 41 data measurement points for 1247 days. First we present the results of history matching in three different scenarios. The tuning parameters used in this section are reported in Table 8. These parameters remain the same in all runs for different scenarios. Total number of simulations for each case is 1380. We start with 30 ants and continue the optimization for 45 generations. Fig. 15 shows the sampling history of ant colony optimization during assisted history matching. The left side figure shows the first case, where we use only 4 oil rate measurements. The middle figure presents the sampling behavior of the algorithm when we consider 23 measurements and the right figure shows the performance of ACOR when we have full oil production history of the field with 41 data points during 1247 days. We can see in Fig. 15 that as we use more data points in the objective function, ACOR algorithm converges faster to low misfit areas of search space. In the first case, where we use only 4 measurements corresponding to 181 days, ACOR searches for the regions of interest until the very end stages of the optimization process and it is only after around 1000 simulations that it finds low misfit models. Also we can see that final history matched models home in on different regions of the search space. 3.2. ACO in high dimensional search spaces In the previous section we applied ACOR to a simple history matching example. In this part we extend the ACOR application to a more challenging problem. We use ACOR for history matching of PUNQ-S3 reservoir model which has more unknown parameters and multivariate production data. PUNQ (Production forecasting with Uncertainty Quantification) was a joint industrial–academic project with the aim of developing efficient history matching and uncertainty quantification methods. The PUNQ-S3 reservoir model is a popular benchmark model to test and compare novel methods developed for history matching and uncertainty quantification. Many authors have published the results of their research for PUNQ-S3 reservoir model. Soleng (1999) used a steady state genetic algorithm for history matching of this model. Manceau et al. (2001) presented an integrated method for history matching and uncertainty analysis based on Fast Fourier TransformMoving Average (FFT-MA) and gradual deformation techniques. Mantica et al. (2002) coupled chaotic optimization with gradual deformation and Demyanov (2004) used Neighbourhood Algorithm (NA) coupled with a geostatistical framework. Gao et al. (2007) applied two types of Simultaneous Perturbation Stochastic Approximation (SPSA) algorithm to history matching of PUNQ-S3 reservoir. 3.2.1. PUNQ-S3 reservoir description PUNQ-S3 reservoir simulation model, as described by Floris et al. (2001) is a 5-layer model. The top depth of PUNQ-S3 reservoir is 2430 m. It has a dip angle of about 1.5° and is bounded by a fault to the east and south and a relatively strong aquifer on the north and west provides a pressure support. Because of this pressure support, no injection wells have been drilled in this reservoir. There is also a small gas cap in the PUNQ-S3 reservoir model in layer 1. There is no well completed in this layer because of effect of free gas production on recovery of the reservoir. Six production wells are marked with black dots as it can be seen in Fig. 16. These wells are located near the initial gas–oil contact. Producers 1 (PRO-1), 4 (PRO-4) and 12 (PRO-12) are perforated in layers 4 and 5. Producers 5 (PRO-5) and 11 (PRO-11) have been completed in layers 3 and 4 and producer 15 (PRO-15) has been perforated only in layer 4. There are positions for five infill wells (X1–X5) which are shown as white dots in the figure. PRO-4 has been completed near aquifer and water breakthrough has been observed in

Y. Hajizadeh et al. / Journal of Petroleum Science and Engineering 77 (2011) 78–92

85

400 0

200

Misfit

600

800

Misfit of history matching

1 2 3 4 5 6 7 8 9

11 13 15 17 19 21 23 25 27 29 31 33 35 37 39

Generation Fig. 10. Boxplot for misfit values obtained by ACOR in Teal South reservoir.

7th year. Free gas production starts in 4th and 5th year in PRO-1 and PRO-4. PUNQ-S3 model contains 19 × 28 × 5 grid blocks, which about two third of the grid blocks (1761) are active. The grid blocks have equal 180 m sides in x and y directions. The reservoir simulation case has been modeled with corner point geometry with Carter–Tracy aquifer type. Complete data set for this reservoir is available online (PUNQ, 2010). The production history of the six wells is summarized as the following periods: (1) An extended well testing period during the first year. (2) Shut-in period for the following three years (3) 12 year production period with fixed production rate at 150 sm3/day (4) Shut in for each well for the period of two weeks in each production period for the purpose of testing the well. Simulated production history for the first 8 years from six wells has been generated by the Netherlands Organization for Scientific Research (TNO). The production history includes pressure, water cuts and gas–oil ratios for each well. In order to reflect the real world measurement errors, Gaussian noise has been added to the well data. After 8 years of production, there are two different recovery scenarios. The first strategy is to continue the current

production for another 8.5 years with the same 6 original wells and the second alternative is adding 5 new wells (X1–X5) with their position being marked with white dots in Fig. 16 and continuing the production for the total period of 16.6 years. In this paper we are only considering the first scenario; continuing the production with the original wells. We use ACOR to minimize the following objective function in PUNQ-S3 case (Barker et al., 2001):    12 0 obs k sim k   oij t −oij t ; p 1 1 1 obs A ∑ ∑ ∑ @wijk SoS o ; p = nw i np j nt k σijk

where nw is number of wells with subscript i running over the wells, and np is the number of production data types with subscript j running over them. Subscript k runs over production data report times and nt is the respective number of samples. Observed data (oobs) and simulated ones (osim) for each of the parameters (p) are being reported at time steps tk with measurement error of σ. At each time step for the parameters, there are extra weighting factors denoted with w. These weights reflect the importance of some of data types at some specific

400 200 0

Misfit

600

800

Misfit of history matching

1 2 3 4 5 6 7 8 9

ð5Þ

11 13 15 17 19 21 23 25 27 29 31 33 35 37 39

Generation Fig. 11. Boxplot for misfit values obtained by NA in Teal South reservoir.

86

Y. Hajizadeh et al. / Journal of Petroleum Science and Engineering 77 (2011) 78–92

0

200

400

600

800

1000

P2

P1

1.0 0.8 0.6 0.4 0.2 0.0

P3

P4

P5

logkvkh

1.0 0.8

Parameter value

0.6 0.4 0.2 0.0 1.0 0.8 0.6 0.4 0.2 0.0

Rock compressibility

Aquifer strength

1.0 0.8 0.6 0.4 0.2 0.0 0

200

400

600

800

1000

Simulation Fig. 12. Sampling history for ACO-TS-1.

0

200

400

P1

600

800

1000

P2 1.0 0.8 0.6 0.4 0.2 0.0

P3

P4

P5

logkvkh

1.0 0.8

Parameter value

0.6 0.4 0.2 0.0 1.0 0.8 0.6 0.4 0.2 0.0

Rock compressibility

Aquifer strength

1.0 0.8 0.6 0.4 0.2 0.0 0

200

400

600

800

1000

Simulation Fig. 13. Sampling history for ACO-TS-2.

Y. Hajizadeh et al. / Journal of Petroleum Science and Engineering 77 (2011) 78–92

0

200

400

P1

600

800

87

1000

P2 1.0 0.8 0.6 0.4 0.2 0.0

P3

P4

P5

logkvkh

1.0 0.8

Parameter value

0.6 0.4 0.2 0.0 1.0 0.8 0.6 0.4 0.2 0.0

Rock compressibility

Aquifer strength

1.0 0.8 0.6 0.4 0.2 0.0 0

200

400

600

800

1000

Simulation Fig. 14. Sampling history for NA-TS-2.

time steps and are specifically indicated in the provided online dataset. 3.2.2. PUNQ-S3 parameterization As mentioned earlier, the PUNQ-S3 case has five layers. We parameterize the PUNQ-S3 model with nine homogenous regions per layer. We use ACOR to estimate the porosities in each homogenous region and layer of the reservoir. Five layers times nine regions per layer makes 45 porosity values that should be estimated in the assisted history matching framework. From least square fitting of the well data crossplots, the following deterministic relationships have been obtained (Boss, 1999). These correlations assume a linear relationship between porosity and logarithmic horizontal permeability and between horizontal and vertical permeability. lnðkh Þ = 0:77 + 9:03ϕ

ð6Þ

kv = 3:124 + 0:306kh

ð7Þ

After estimating 45 porosity values, we use the above correlations to determine the remaining horizontal and vertical permeability

Table 8 Summary of parameters used in ACOR. Number of ants (m)

Archive size (k)

Search locality (q)

Pheromone evaporation (ξ)

30

30

0.3

0.3

values in the regions and layers of the reservoir. The initial ranges for unknown parameters are given in Table 9. These values are based on the selected based on geological description of the reservoir. Layers 1, 3 and 5 are high quality and layers 2 and 4 have lower porosity and permeability values. 3.2.3. History matching results of PUNQ-S3 For history matching of PUNQ-S3 model, we consider 5 different ACOR cases. The tuning parameters for each algorithm are reported in Table 10. The best misfit value obtained using ACOR algorithm is 1.83. Fig. 17 presents the match result for selected BHP, GOR and WWC data in PUNQ-S3 reservoir. To compare the results of ACOR with NA, we also consider 5 different Neighbourhood Algorithm cases. For each case, we have tried different nr values. We report the best results for NA cases and their tuning parameters in Table 11. Further to these cases, we also compare the results of ACOR with the reported results obtained from coupling Neighbourhood Algorithm with a geostatistical framework (Demyanov, 2004). Demyanov used two different setups, one with 15 unknown parameters and the second one with 60 unknown parameters. Neighbourhood Algorithm was used to tune the geostatistical parameters with the aim of getting different realizations of the reservoir model that match the observed production data. We name these cases NA-GS-1 and NA-GS-2 which have 20,200 simulations. Fig. 18 shows the results of comparing ACOR with NA and NA-GS. Comparing ACOR with NA and NA-GS cases, one can see that there is a significant improvement in terms of the minimum misfit obtained by ACOR. In the Teal South example, the difference between best misfits obtained by ACOR and NA was marginal (ACOR = 14.71,

Y. Hajizadeh et al. / Journal of Petroleum Science and Engineering 77 (2011) 78–92

Fig. 15. Sampling history for ACOR in different scenarios.

88

Y. Hajizadeh et al. / Journal of Petroleum Science and Engineering 77 (2011) 78–92

89

Fig. 16. PUNQ-S3 reservoir simulation model and top surface map with well positions.

NA = 14.83). As we go to higher dimensional search spaces, this difference is becoming larger and this shows the success of ACOR even when the number of unknown parameters is high. In all of the 5 setups, ACOR obtained lower misfit values in direct comparison with NA from same number of simulations (3000). Also looking at the results of coupling NA with geostatistics (NA-GS), we can see that the best misfit value of ACOR with 3000 simulations (1.83) is lower than the best result obtained by NA-GS after 20,200 simulations (2.92). 4. Uncertainty analysis using NA-Bayes algorithm After exploring the search space and generating the ensemble of models by ACOR and NA, the next step is the calculation of the probabilities of the models in ensemble. In this work, the NA-Bayes (NAB) algorithm (Sambridge, 1999b) is used for posterior inference. NAB is a Markov Chain Monte Carlo (MCMC) method which builds an approximation for the posterior probability distribution (PPD) using a Gibbs sampler. It uses Voronoi cells to represent the model space and to interpolate the PPD of unknown points in the search space. This interpolation is based on the assumption that misfit value is constant over each Voronoi cell. The main benefits of using NAB for this purpose are that: (1) It infers the information from all the models in the ensemble to evaluate the posterior probability of the ones contributing to the

Table 9 Initial ranges for parameters in PUNQ-S3 reservoir. Layer 1 2 3 4 5

Porosity 0.15–0.3 0.05–0.15 0.15–0.3 0.1–0.2 0.15–0.3

Horizontal permeability (md)

Vertical permeability (md)

133–3013 16–133 133–3013 47–376 133–3013

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

confidence prediction. Hence, models with different goodness of fit contribute to the differently to the uncertainty prediction. (2) There is no need to run forward reservoir simulations for all the models generated by the sampling algorithm, but only for the ones resampled by NAB. This helps to save computational resources.

4.1. Uncertainty analysis in Teal South In this part we compare three cases of history matching using ACOR algorithm with different number of oil rate measurements. We submitted the ensemble of 1380 models generated by ACOR to NAB algorithm in order to determine the Bayesian credible intervals (p10– p50–p90) for the oil production after 1247 days. For the first and second scenarios where history matching is done for 181 and 699 days, the cumulative oil recovery is predicted after 1247 days, running forecast simulations. Fig. 19 shows the total recovery prediction from Teal South reservoir for the three scenarios studied in this paper. The blue dashed line is the true total production after 1247 days which is computed using the observed production data in the field. As we see in Fig. 19, p10–p90 range in all scenarios cover the true oil production obtained after 1247 days. In the case of the first scenario we have wide uncertainty bands. This is because the history matching is done with having access to only 6 data points during

Table 10 Algorithm setup and final misfit values for different ACOR cases in PUNQ-S3 reservoir. Case

m

k

q

ξ

Total simulations

Misfit

ACO-PU-1 ACO-PU-2 ACO-PU-3 ACO-PU-4 ACO-PU-5

50 50 50 100 50

50 50 50 100 500

0.4 0.7 0.5 0.9 0.05

0.7 0.7 0.9 0.5 0.2

3000 3000 3000 3000 3000

1.83 1.90 2.26 3.05 2.75

90

Y. Hajizadeh et al. / Journal of Petroleum Science and Engineering 77 (2011) 78–92

160

GOR (Well 4)

120 100 80

Gas/Oil Ratio

140

ACOR Observations

60

ACOR Observations

0

1000

2000

3000

4000

5000

0

6000

3000

4000

GOR (Well 11)

Water cut (Well 11)

ACOR Observations

6000

5000

6000

0.4 0.1

0.2

0.3

Water cut

100 80

5000

ACOR Observations

0.5

120

2000

Time (days)

0.0

60

Gas/Oil Ratio

1000

Time (days)

0.6

BHP (Bar)

100 120 140 160 180 200 220 240

BHP (Well 1)

0

1000

2000

3000

4000

5000

6000

0

1000

Time (days)

2000

3000

4000

Time (days)

Fig. 17. BHP, GOR and WWC match results for selected wells in PUNQ-S3 reservoir.

181 days. As we move to the second and third scenarios and use more measurement points in history matching, the uncertainty intervals become narrower and the p50 value of the predictions becomes closer to the true total production. 4.2. Uncertainty analysis in PUNQ-S3 History matching of PUNQ-S3 is not the only challenge of this example. Prediction of the total recovery from this reservoir after 16.5 years is another important objective for the PUNQ-S3 model. Many history matching and uncertainty quantification methods fail to correctly predict the total recovery value (Floris et al., 2001). Methods used for history matching and uncertainty quantification of PUNQ-S3 are different in the following 3 aspects. The first difference is the parameterization of the model, the second is that they use different methods to find good matching models (gradient, stochastic) and the third is the approach they take to quantify the

uncertainty. These three aspects are equally important in reservoir engineering studies. For quantification of uncertainty in prediction of total oil recovery in PUNQ-S3 reservoir, we follow the same procedure for Teal South reservoir. We submit the ensemble of models generated by ACOR to the NAB algorithm and obtain the Bayesian credible intervals. Fig. 20 shows the comparison between our work and the approaches that use a gradient-based optimization technique for finding good fitting models (Floris et al., 2001) and Fig. 21 compares the stochastic-based methods used in previous publications (Floris et al., 2001; Demyanov, 2004; Holmes et al., 2007) and this work. In Figs. 20 and 21 we see that the results of five ACOR cases cover the truth production rate and their p-50 values are very close to the exact answer (3.87 × 106). We should note that some methods cannot predict the correct production value, although they have a relatively wide uncertainty range.

5. Conclusions Table 11 Tuning parameters and final misfit values obtained in each setup of NA algorithm for PUNQ-S3. Case

nsi

ns

nr

Iterations

Total simulations

Misfit

NA-PU-1 NA-PU-2 NA-PU-3 NA-PU-4 NA-PU-5

50 50 50 100 500

50 50 50 100 50

10 25 50 100 10

60 60 60 30 51

3000 3000 3000 3000 3000

4.26 5.26 4.07 4.32 4.08

In this paper ant colony optimization (ACO) has been successfully applied to the problem of reservoir history matching and uncertainty quantification. We used ACOR algorithm to generate multiple history matched models for two reservoir modes: a simple model of a real reservoir (Teal South), and a more complex synthetic model (PUNQ-S3) with more degrees of freedom and variables. We studied the value of information in ant colony optimization and how additional data points affect the history

Y. Hajizadeh et al. / Journal of Petroleum Science and Engineering 77 (2011) 78–92

91

Fig. 18. Comparison of best misfit values in ACOR, NA and NA-GS for PUNQ-S3 case.

matching results and uncertainty of the prediction. The specific conclusions from this study are: (1) ACOR is capable of generating multiple history matched models which makes it a promising tool for using in reservoir engineering. (2) ACOR converges to good fitting regions of the search space even when it is started from different initial points. (3) For both reservoir case studies, ACOR is able to find better models in terms of misfit value comparing with NA. Also ACOR can do this from fewer simulations. In the Teal South reservoir the difference between the performance of the ACOR and NA was marginal due to a relatively simple model with just eight free parameters and the univariate objective function based on the data from a single-well. In high dimensional problems (PUNQ-S3), ACOR obtains much better models in comparison with NA. (4) Using too few data points in evolutionary optimization may cause convergence problems and mislead the search of these algorithms (5) Assessing the effect of additional data clearly showed the impact of this information on narrowing down the range of likely model parameters. The additional data also resulted in

Fig. 20. Uncertainty intervals in PUNQ-S3 model using ACOR and its comparison with gradient methods.

×1000 STB

750

Total recovery (STB)

725 700 675 Truth

650 625 600 575 550 ACO-181

ACO-699 P-10

p-50

ACO-1247 P-90

Fig. 19. Total recovery prediction for ACOR and NA runs in Teal South reservoir.

Fig. 21. Uncertainty intervals in PUNQ-S3 model using ACOR and its comparison with stochastic methods.

92

Y. Hajizadeh et al. / Journal of Petroleum Science and Engineering 77 (2011) 78–92

reduced uncertainty bands in forecasting cumulative oil produced. (6) ACOR provides the confidence bounds which comfortably include the true solution both for Teal South and PUNQ-S3 problems. Acknowledgments The support for Uncertainty Quantification project at Heriot Watt University is provided by BG, BP, ConocoPhillips and JOGMEC which are gratefully acknowledged here. We also thank Schlumberger GeoQuest for providing ECLIPSE reservoir simulation software. References Anterion, F., Eymard, F., 1989. Use of parameter gradients for reservoir history matching, SPE 18433. SPE Symposium on Reservoir Simulation, Houston, Texas, U.S.A. 6–8 February. Barker, J.W., Cuypers, M., Holden, L., 2001. Quantifying uncertainty in production forecasts: another look at the PUNQ-S3 problem, SPE 74707. SPE J. 6 (4), 433–441. Bonabeau, B., Dorigo, M., Theraulaz, G., 2000. Inspiration for optimization from social insect behavior. Nature 406 (6791), 39–42. Boss, C., 1999. Production forecasting with UNcertainty Quantification. Netherlands Institute of Applied Geoscience TNO. 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. Chen, W.H., Gavalas, G.R., Seinfeld, J.H., Wasserman, M.L., 1974. A new algorithm for automatic history matching, SPE 4545. SPE J. 14 (6), 593–608. Christie, M., MacBeth, C., Subbey, S., 2002. Multiple history-matched models for Teal South. Lead. Edge 21 (3), 286–289. Demyanov, V., 2004. Uncertainty assessment using geostatistics in reservoir forecasting: use of NA on PUNQ-S3 problem, progress report. Institute of Petroleum Engineering, Heriot Watt University, Edinburgh, UK. Dorigo, M., 1992. Learning and natural algorithms, PhD thesis, Politecnico di Milano, Milan, Italy. Dorigo, M., DiCaro, G., Gambardella, L., 1999. Ant algorithms for discrete optimization. Artif. Life 5 (2), 137–172. Erbas, D., 2006. Sampling strategies for uncertainty quantification in oil recovery prediction, PhD Thesis, Institute of Petroleum Engineering, Heriot Watt University, Edinburgh, UK. 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. 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. Gao, G., Li, G., Reynolds, A.C., 2007. A stochastic optimization algorithm for automatic history matching, SPE 90065. SPE J. 12 (2), 196–208. Hajizadeh, Y., Christie, M., Demyanov, V., 2009. Ant colony optimization algorithm for history matching. SPE 121193, EUROPEC/EAGE Conference and Exhibition, 8–11 June, Amsterdam, The Netherlands. 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.

Jafarpour, B., McLaughlin, D.B., 2007. History matching with an ensemble kalman filter and discrete cosine parameterization, SPE 108761. SPE Annual Technical Conference and Exhibition, Anaheim, California, U.S.A. 11–14 November. Kathrada, M., 2009. Uncertainty evaluation of reservoir simulation models using particle swarms and hierarchical clustering, PhD thesis, Institute of Petroleum Engineering, Heriot Watt University, Edinburgh, United Kingdom. Liu, N., Oliver, D.S., 2005. Critical Evaluation of the Ensemble Kalman Filter on History Matching of Geologic Facies, SPE 92867. SPE Reservoir Simulation Symposium. The Woodlands, Texas, U.S.A. 31 January-2 February. Lodoen, O.P., Omre, H., 2008. Scale-corrected ensemble kalman filtering applied to production-history conditioning in reservoir evaluation, SPE 111374. SPE J. 13 (2), 177–194. Manceau, E., Mezghani, M., Zabala, I., Roggero, F., 2001. Combination of experimental design and joint modeling methods for quantifying the risk associated with deterministic and stochastic uncertainties — an integrated test study, SPE 71620. SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, U.S.A. 30 September-3 October. Mantica, S., Cominelli, A., Mantica, G., 2002. Combining global and local optimization techniques for automatic history matching production and seismic data, SPE 78353. SPE J. 7 (2), 123–130. 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. 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. SPE Reservoir Simulation Symposium, Woodlands, Texas, U.S.A. 2–4 February. Okano, H., Pickup, G., Christie, M., Subbey, S., Sambridge, M., 2006. Quantification of uncertainty in coarse-scale relative permeabilities due to sub-grid heterogeneity. 10th European Conference on the Mathematics of the Oil Recovery, paper A034, Amsterdam, The Netherlands. 4–7 September. PUNQ-S3, 2010. Department of earth science and engineering. Imperial College Londonhttp://www3.imperial.ac.uk/earthscienceandengineering/research/perm/ punq-s3model. Last Accessed October 2010. 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. Socha, K., Dorigo, M., 2008. Ant Colony optimization for continuous domains. Eur. J. Oper. Res. 185 (3), 1155–1173. Soleng, H.H., 1999. Oil reservoir production forecasting with uncertainty estimation using genetic algorithms. Evolutionary Computation, CEC 99, Washington D.C, U.S.A. 6–9 July. Sousa, S.H.G., Maschio, C., Schiozer, D.J., 2006. Scatter search metaheuristic applied to the history matching problem, SPE 102975. SPE Annual Technical Conference and Exhibition. San Antonio, Texas, U.S.A., 24–27 September. Subbey, S., Christie, M., Sambridge, M., 2003. A strategy for rapid quantification of uncertainty in reservoir performance prediction, SPE 79678. SPE Reservoir Simulation Symposium, Houston, Texas, U.S.A. 3–5 February. Sultan, A.J., Ouenes, A., Weiss, W.W., 1994. Automatic history matching for an integrated reservoir description and improving oil recovery, SPE 27712. Permian Basin Oil and Gas Recovery Conference, Midland, Texas, U.S.A. 16–18 March. Valjak, M., 2008. History matching and forecasting with uncertainty, challenges and proposed solutions for real life field application, PhD Thesis, Institute of Petroleum Engineering, Heriot Watt University, United Kingdom. Yang, C., Nghiem, L., Card, C., 2007. Reservoir model uncertainty quantification through computer-assisted history matching, SPE 109825. SPE Annual Technical Conference and Exhibition, Anaheim, California, U.S.A. 11–14 November.