Abstract

The optimization of oil field development scheme considering the uncertainty of reservoir model is a challenging and difficult problem in reservoir engineering design. The most common method used in this regard is to generate multiple models based on statistical analysis of uncertain reservoir parameters and requires a large number of simulations to efficiently handle all uncertainties, thus requiring a huge amount of computational power. In order to reduce the computational burden, a method which combines reservoir simulation, an economic model, polynomial chaos expansion with response surface methodology, and Levy flight particle swarm optimization (LFPSO) algorithm is proposed to determine the optimal injection-production parameters with reservoir uncertainties at a reasonable computational cost. This approach is applied to a five-spot well pattern optimization design for obtaining the optimal parameters, including oil-water well distance, injection rate, and bottom hole pressure, while considering the uncertainties of porosity, permeability, and relative permeability. The results of the case study indicated that the integrated approach is practical and efficient for performing reservoir optimization with uncertain reservoir parameters.

1. Introduction

The main goal of oil field development is to obtain maximum waterflooding efficiency and ultimate reservoir recovery with minimum capital investment and operational cost, which, in other words, means to obtain the greatest economic benefits using an optimized design scheme [1, 2]. The optimal field development strategy depends on reservoir geometry and physical parameters of rocks and fluids. The subsurface is usually poorly understood, as only limited information about the reservoir can be obtained from well logging, well testing, and laboratory experiments conducted on core samples. Due to the sparsity and noise of field information, sources of uncertainties in reservoir engineering are almost infinite and can be found anywhere within the reservoir modeling workflow. Reservoir modeling workflow includes the static model, upscaling, fluid flow modeling, production data integration, production scheme development, and economic evaluation [3]. Uncertainty in reservoir parameters must be accounted for during optimization due to the huge cost involved in oil field development.

In recent years, the oil industry has given great importance to optimization of oil field schemes and reservoir uncertainty analysis [47]. With the advent of computer technology, most researchers adopt the method of combining numerical simulations with relevant mathematical optimization algorithms or orthogonal experiment design when dealing with production optimization problems of waterflooding. In order to avoid a time-consuming, tedious, and complicated problem of large number of simulations, the proxy model and more effective algorithms are introduced into the optimization of injection-production parameters. Some progress has also been made with regards to the quantification of uncertainty of reservoir parameters. Geostatistics provides a set of methods to study the spatial variability and uncertainty through parameter estimation and distribution tests. Many well-known methods have been proposed for performing the propagation of uncertainty. Some of them include nested simulations, decision tree analysis, Monte Carlo sampling, differential analysis, response surface methodology, Bayesian analysis, fuzzy logic, and interval analysis [815]. These methods are rigorous and systematic and ensure reliable results in terms of quantifying uncertainty. However, a common characteristic of all these approaches is the very large number of simulations required to efficiently handle multiple reservoir models generated according to the uncertainty statistics, which will then require huge computational power and slow convergence. Therefore, it is difficult to use these approaches in evaluating the uncertainty benefit scheme in an actual reservoir. Furthermore, it is difficult to realize the optimal design of reservoir development scheme under the condition of reservoir uncertainty.

The combination of stochastic modeling theory, fuzzy mathematics, and other methods with response surface methodology provides an effective means for studying uncertainty optimization design in oil field development. Dejean and Blanc [10] proposed integration of experimental design, response surface, and Monte Carlo methods to optimize production scheme and applied their proposed methodology to a field in assessing the effects of uncertainties on oil production. Cullick et al. [16] integrated reservoir simulation, an economic model, and a Monte Carlo algorithm with a global optimization search algorithm to identify optimal reservoir planning and management decision alternatives under uncertainties. In order to reduce the sensitivity and risk of production strategy to account for parametric uncertainties in reservoir models, robust optimization (RO) [1721] and proper orthogonal decomposition (POD) method [2224] have been established based on an ensemble of realizations reflecting the uncertain characteristics of reservoirs. Since simulation cost associated with constructing uncertain response surface model is not realistic for risk assessment in complex reservoir modeling, it is difficult to achieve by a single RO method, and a compromise has to be made between the reliability of approaches and their related costs.

At present, the methods to quantify uncertainty that are based on polynomial chaos expansion (PCE) have been successfully applied in the area of stability and control [25], solid mechanics [26, 27], electronic circuits [28], and computational fluid dynamics [29]. There are also numerous applications of PCE in petroleum engineering problems [3036]. PCE techniques can be seen as an efficient and mathematically optimal way to construct the response surface of a model with uncertain parameters and has its origin in the Wiener chaos expansion [37], which is based on Hermite polynomials and normally distributed random variables.

In this paper, polynomial chaos expansion method is used to evaluate the uncertainty of waterflooding oil field scheme design. A method coupling reservoir numerical simulation models with polynomial chaos expansion as well as response surface methodology and particle swarm optimization (PSO) algorithm with Levy flight (LFPSO) is proposed to determine the optimal value of the injection-production parameters under uncertainty. The proposed method can take a significantly reduced computational cost to quantify the uncertainty of economic efficiency with accuracy similar to that of computer numerical simulations models using the Monte Carlo method.

2. Methodology

In order to optimize the injection-production parameters of oil field development with uncertainty, the simplest and direct method is the Monte Carlo simulation method. Compared with other techniques related to the propagation of uncertainty, Monte Carlo simulation [9, 38, 39] is based on mathematical and statistical theory and is the most widely used uncertainty method. However, in order to ensure the accuracy of calculations, Monte Carlo simulation requires a large number of samples. For complex reservoir numerical models, Monte Carlo simulation requires a long time to execute optimal injection-production parameters. In order to reduce the computational complexity of injection-production parameters’ optimization with uncertainties, polynomial chaos expansion of net present value (NPV) with response surface methodology is coupled with LFPSO to find the optimal injection-production parameters for reservoir development.

The paper is structured as follows: in Section 2.1, the uncertain parameters of the reservoir as well as injection-production that need to be optimized are analyzed. Also, the performance objective function (i.e., NPV), based on uncertain injection-production parameters, is established. Section 2.2 presents a method of constructing a surrogate model using PCE. This includes the determination of orthogonal polynomial bases, the coefficients of PCE, a method to evaluate the convergence of PCE, and methods for further establishing the response surface model based on the converged PCE of NPV. In Section 2.3, the LFPSO optimization algorithm is introduced and the procedure and methods employed in optimizing the injection-production parameters based on the converged PCE of NPV with response surface methodology and LFPSO are discussed in detail.

2.1. Parameter Analysis and the Developed NPV Function

Most traditional optimization methods for the development of oil field waterflooding focus on finding the optimal design parameters to maximize some of the predefined reservoir performance metrics, such as the NPV or the cumulative oil production [6, 4042]. NPV is the difference in the present value of cash inflows and outflows associated with a project, where the future value of money is discounted to its present value. It is the most widely used objective function in reservoir optimization. Because there are many uncertain parameters in a reservoir, the optimization and evaluation of oil field development waterflooding schemes need to deal with very complex models. According to the material balances of the two phases and following Darcy’s law for multiphase fluids, the governing equation of oil-water two-phase flow in reservoir system can be given aswhich is constrained by the equation

The subscripts and o stand for water phase and oil phase, respectively; is the porosity; and are the water-phase saturation and oil saturation, respectively; K is the absolute permeability; and are the relative permeability of water and oil, respectively, which are functions of and ; and are the water density and oil density, respectively; is the gravity vector; stands for time; and and are the sources/sinks, respectively.

According to the governing equations, porosity, absolute permeability, and relative permeability as a function of intrinsic properties of reservoir and fluid have a significant impact on the recovery of oil field waterflooding. However, these parameters are usually sampled at few locations, for which reason uncertainty becomes dominant, especially considering the spatial heterogeneity of reservoirs. Therefore, the uncertain parameters considered in waterflooding design in this paper are porosity, permeability, and relative permeability. According to some previous studies, the relative permeability of each phase can be described using a function of saturation of each phase in the porous medium. The most common type of such a relationship is the Corey model [43], which takes the following form:wherewhere and denote the water and oil relative permeability, is the residual oil saturation, and is the connate water saturation. Furthermore, is the water saturation, is the normalized water saturation considering the irreducible saturations of both phases, and and are the indices of water and oil, respectively.

According to the Corey model, the relative permeability curve of oil-water can be characterized using six parameters of , , , , , and . In this paper, only and are considered as uncertain parameters, while other parameters are the values determined from experiments. One hundred realizations of oil-water relative permeability curves are shown in Figure 1. They are obtained using the Corey model considering and that satisfies the uniform distribution within the range of 2–4. In short, the uncertain reservoir parameters are considered in this paper, in which and are considered to obey the Gaussian distribution in a certain range, whereas and obey the uniform distribution.

Four production wells in the corners of the model at constant bottom hole pressure (BHP) and a central water injector at constant injection rate are considered in this paper. The bottom hole pressure (), injection rate (Q), and distance (L) between the injection wells and production wells are important components of injection-production parameters and will directly affect the development and ultimate recovery from the reservoir. Therefore, injection-production parameters are considered in the design optimization of oil field development waterflooding in this paper. When considering reservoir uncertainty, NPV is determined based on the design parameters and uncertainty parameters . The NPV function is computed as follows [44]:where , , and are the flow rates of oil, water produced, and water injected for well at the output interval , respectively, and represents the length (in days) of each of time steps in the simulation. The oil price and the cost of water produced and injected are denoted by , , and , respectively. Furthermore, is the annual discount rate.

2.2. Construction of Surrogate Model Using PCE

Considering the uncertainty of reservoir parameters, the uncertain and injection-production parameters constitute the input variables of reservoir’s numerical simulations. A large number of input parameters can lead to unbearable and huge computational load. The PCE can be regarded as a mathematical algorithm to construct a model response in the form of a high dimension polynomial to bypass the reservoir’s fluid flow simulator for quantifying uncertainties, which ultimately provides fast evaluation of NPV in an uncertain environment. Conventional PCE using Hermite polynomials [45], which are orthogonal with respect to standard normal distribution, is applicable in solving uncertain parameters following normal distributions. Due to the limited applicability of these methods, the arbitrary polynomial chaos was developed, which can be applied to arbitrary distributions with known moments. The procedure for the construction of surrogate model of NPV using PCE relates to the uncertainty in input parameters and is presented in Figure 2.

2.2.1. Polynomial Chaos Expansion of NPV

According to the PCE method [37], NPV associated with uncertain input parameters is represented in a series form as follows:where values are the unknown coefficients of the PCE; is the multi-dimensional random vector, in which each element is associated with uncertain input parameters; and are the univariate orthogonal polynomials and/or multivariate orthogonal polynomials, which are the product of one-dimensional orthogonal polynomials with different variables. For the convenience of calculations, the input parameters can be linearly transformed into random variables with the mean value of 0 and variance of 1, as follows:where is the expectation of and is the variance of . Conversely speaking, if is known, can be obtained using the following equation:

Furthermore, in order to estimate the uncertainty of output, the polynomial chaos expansion of NPV in equation (6) can be truncated to a certain order. The d-order truncated polynomial chaos expansion of NPV with respect to n-dimensional random variable is rewritten in the form, as shown in the following equation:with . Because the parameters are independent of each other, the parameters obtained through linear transformation of equation (7) are also independent of each other. Therefore, can be decomposed into a product of several one-dimensional orthogonal polynomials as depicted in the following equation:

In order to determine the polynomial chaos expansion of NPV, it is necessary to determine the one-dimensional orthogonal polynomial and polynomial chaos expansion coefficients .

2.2.2. Selection of Orthogonal Polynomial Chaos Bases

The construction of orthogonal polynomial bases is an important part of polynomial chaos expansion. Based on the type of distribution of input parameters, the orthogonal polynomial basis of the surrogate model of NPV is determined. Askey and Wilson [46] provided the orthogonal polynomials for different distributions, and the results are presented in Table 1.

Based on the work of Oladyshkin and Nowak [35], a method for determining the base of one-dimensional orthogonal polynomial is given when is arbitrarily distributed.

The k-order origin moment of can be defined as follows:

The coefficients of orthogonal polynomial of the order can be determined by solving the following matrix equation [35]:

According to equation (12), statistical moments are the only form of information required about the input distributions and it is necessary to obtain the -order origin moments of random parameters to determine the orthogonal base of d-order polynomial chaos expansion. In order to avoid divergence of polynomial values, a normalized basis has more useful properties and can be written as follows:

Therefore, the normalized orthogonal base of is as follows:

2.2.3. Determination of Unknown Coefficients of PCE

The polynomial chaos coefficients determine the accuracy of the output NPV estimation. Based on the work of Sudret [47], the regression nonintrusive method is chosen to calculate the unknown coefficients of PCE for NPV in this work. After samples are determined, the input parameters can be obtained using equation (8). Then, the samples of input parameter are brought into reservoir numerical simulator to obtain samples . Once the samples and are determined, the polynomial chaos expansion coefficient can be calculated using the following equation:

Let

Then, that makes equation (15) valid should be

The input samples are extracted in turn according to the above order until the matrix becomes reversible. When the matrix is ill-conditioned, the singular value decomposition (SVD) algorithm [48] can be used to determine .

2.2.4. Evaluation of Convergence of the PCE Model of NPV

In order to ensure accuracy, it is necessary to evaluate the convergence of PCE of NPV. The uncertainty between the PCE and higher-order PCE can be compared to evaluate the convergence of PCE. Once the PCE of NPV is obtained, the expected value (), variance (), probability density function (PDF), and cumulative distribution function (CDF) of NPV can be estimated. According to the definition of expectation and orthogonality of PCE, the expectation of NPV is given as follows:

Because of the normality of and the independence of , the variance of NPV can be calculated as follows:

Because the PCE of NPV is the sum of orthogonal polynomials, its computational cost of Monte Carlo simulation can be neglected. Therefore, sampling-based statistical methods can be used to estimate the PDF and CDF of NPV, which are calculated using PCE. In this paper, the convergence of PCE is adjudged based upon the difference between CDFs of NPV calculated using PCE of adjacent order. If the CDF of NPV, estimated using PCE of order n, has little deviation from that estimated using PCE of order n + 1, then the PCE of the NPV of order n is regarded as convergent. The difference between the CDFs can be measured using the root mean square deviation (RMSD). If RMSD is less than 5%, the PCE of NPV is considered convergent. Otherwise, a higher-order PCE needs to be constructed until it converges. The RMSD is calculated using the following equation [44]:

2.2.5. PCE Surrogate Model

The NPV of oil field development is affected by both the uncertain parameters and injection-production parameters. Equation (5) is rewritten as follows:where is the set of injection-production parameters, is the number of injection-production parameters, is the set of uncertain parameters, and is the number of uncertain parameters.

The two-stage nested PCE developed consists of two parts: external circulation, which sampled injection and production design parameters through Latin hypercube sampling (LHS) [49] to obtain the injection-production parameters set samples , and internal circulation, which sampled uncertain parameters through LHS to obtain the uncertain parameter set samples . According to equation (14), the set of injection-production parameters and uncertain parameters is transformed to , and . The uncertain parameter set U and injection-production parameters are incorporated into reservoir numerical simulator to obtain NPV samples . According to the theory of PCE, can be expanded as

Once the random base is selected, the coefficient will become a definite value. The value of depends only on , due to which, can be expressed as . Based on response surface model (RSM), can be approximated using the following equation:

In equation (23), is the base of RSM. are the coefficients of RSM. represents the transpose operator. Changing equation (23) into equation (22), can be expressed as follows:

The RSM can also be developed using LHS, and the number of selected sampling points is at least [50]. Assuming that the number of sampling points is , the coefficient vector of RSM can be obtained through a least square method as follows:where

This can be represented in the matrix form aswhere

2.3. Method of Optimizing Design under Uncertainty

After constructing the NPV surrogate model by PCE with response surface as equation (24), LFPSO algorithm combined with reliability optimization under uncertainty was utilized to solve the model.

When there are uncertain variables in the model, the optimization problem becomes an uncertain optimization problem. The optimization model is as follows:where stands for the probability of conforming to constraint and represents the reliability index. The closer the is to 1, the higher the reliability of the system. In this paper, the object function is the surrogate model of PCE with response surface model expressed as equation (24).

2.3.1. LFPSO

LFPSO is used to solve the optimization problems as equation (24) to obtain the optimal injection-production parameters with the best economic benefits. LFPSO algorithm [51] embeds Levy flight into PSO algorithm and overcomes the problem of premature convergence and insufficient global search ability of PSO. Levy flight is a non-Gaussian stochastic process, which is a random walk method extracted from Levy stable distribution. The method ensures that PSO can search globally and more effectively without falling into local stagnation. The LFPSO algorithm sets a threshold for every particle, whereas the particle velocity and its position within the threshold are updated as follows:

When the threshold is exceeded, the particle performs the Levy flight operation. The positions of particles are then redistributed:

The basic flow chart of LFPSO algorithm is shown in Figure 3. In Figure 3, denotes the current number of trials and denotes the threshold set for each particle. It is worth noticing that if the updated particles are better, is reset to 0; otherwise, is added to 1.

2.3.2. Steps and Procedure of Optimization Design

The steps and procedure employed in optimizing the injection-production parameters based on the converged PCE of NPV with response surface methodology and LFPSO are discussed in detail in this section. The procedure for determining the optimal injection-production parameters under uncertainty is shown in Figure 4. The simplification steps of the method are as follows:Step 1. To select uncertain and design variables and determine the value space and distribution type according to the design requirementsStep 2. To arrange the test design combination table by the LHS method and call the reservoir numerical simulator to calculate the response function values (NPV) of each test design combinationStep 3. To construct the PCE model of NPV and uncertain variablesStep 4. To inspect the accuracy of PCE model and rearrange the test design until the accuracy requirements are metStep 5. To construct the PCE with RSM model of NPV and design variables on the basis of the PCE model satisfying the accuracy testStep 6. To analyze the uncertainty of the surrogate model and optimize the design variables based on the surrogate model and LFPSO algorithm

3. Case Study and Analysis

3.1. Example Description

According to the analysis provided in Section 2, oil-water well distance, injection rate, and bottom hole pressure were the design parameters. These design parameters were optimized under uncertain reservoir porosity, permeability, and relative permeability as characterized by the Corey model with uncertain oil-phase and water-phase indices. In this section, a three-dimensional two-phase Eclipse black-oil simulator was used to compute the oil and water production from the production well and the NPV was computed from the fluid production profiles generated from the simulation run associated with corresponding injection-production scheme. A synthetic five-spot 2D reservoir model established by Eclipse and consisting of 100 × 100 × 1 grid cells for a total field size of 1000 × 1000 × 10 m was examined. The reservoir contained two-phase fluid of oil and water. The production mainly relied on the elastic expansion of rock pores resulting from the pressure drops in the formation and the external energy supply from the injected water.

In the synthetic reservoir, the average porosity and permeability, which were both assumed to obey truncated Gauss distribution, were 0.15 and 75 mD, respectively. Because of the boundedness of porosity and permeability, their ranges were given according to realistic values in oil fields. In addition, the water-phase and oil-phase indices for calculating oil-water relative permeability were assumed to obey uniform distribution within the range of 2–4. Meanwhile, these uncertain parameters were assumed to be independent. The distribution and range of uncertain reservoir parameters used in this example are presented in Table 2.

The values of parameters for computing NPV and the remaining reservoir system properties are presented in Table 3.

According to the principle of reservoir engineering and the experience of oil field production and development, the optimum range of injection and production parameters is presented in Table 4.

3.2. Results and Discussion

The order of PCE model should be determined when the PCE is established using internal circulation sampling. Taking the external design variable , the 2nd, 3rd, and 4th order PCE were constructed by using 21, 56, and 126 samples (respectively) to determine the surrogate model of NPV in this example. The CDFs of NPV, which were estimated by coupling the 2nd, 3rd, and 4th order PCE with LHS, are shown in Figure 5. The RMSD between the CDFs of NPV estimated from the 2nd and 3rd order PCE was 7.23%, while based on the 3rd and 4th order PCE, the corresponding value was 3.42% (less than 5%), which showed that the convergence of the 3rd order PCE was satisfactory.

In the external circulation, the surrogate model of NPV was established by taking 32 sample points in the design variable space. In order to verify the accuracy of the 3rd order PCE with the response surface model, 1000 samples of NPV estimated by the surrogate model and calculated using the reservoir numerical simulator with the same sample input parameters are compared in Figure 6. The results showed that the surrogate model of the 3rd order PCE with response surface model is considerably accurate. Therefore, the surrogate model of NPV can be used as the equivalent evaluation model for numerical simulation when using the optimization algorithm to design parameters of oil-water well distance, injection rate, and bottom hole pressure.

In this study, the reliability value of equation (31) is 0.97. The LFPSO method is used to maximize the 3rd order PCE with response surface model of NPV to find the optimal injection-production parameters in the design variable space. After 20 iterations, the optimization results tend to be stable, as shown in Figures 7 and 8.

The optimization results of LFPSO combined with PCE with response surface model are presented in Table 5. In order to verify the effectiveness of the proposed method, two of the three parameters, namely, the bottom hole pressure (), injection rate (Q), and distance between the injection and production wells (L), were taken from Table 5. The other parameter was sampled 1000 times using the MCS method to obtain the relationship between NPV value and design variables, as shown in Figure 9. The results showed that the optimized results obtained using the method proposed in this paper were reliable. In addition, the optimization time of LFPSO combined with PCE is 6563 s, which was much lower than 30217 s required by the MCS method.

Under the optimal injection and production parameters, the uncertainty variables were sampled by using the MCS method and the probability density function curve of NPV was obtained, as shown in Figure 10. It can be seen that, under the optimal injection and production parameters, the expected NPV value was high, though the variance was small. This showed that, under the influence of uncertain parameters, the proposed method can be effectively applied to design optimization of injection-production parameters in waterflooding.

4. Conclusions

In order to optimize the injection-production parameters of waterflooding at a reasonable computational cost and at an acceptable reliability level, a surrogate model of NPV constructed by PCE with RSM is used to solve the problem of complex nonlinear implicit relationship between NPV and design variables considering the random uncertainty parameters of reservoirs. In addition to solving the problem of uncertainty description and quantification between NPV and design parameters, the surrogate model has been applied to the optimization of a synthetic five-spot 2D reservoir combined with LFPSO algorithm. It improves the efficiency of the method to analyze the optimum injection-production parameters of reservoirs, which correspond to the reliability of NPV arbitrary values. The method proposed in this paper can also be applied to other oil field optimization problems considering any uncertain factors. Compared with the traditional technology, the proposed method can effectively deal with the uncertainty of the development of oil field at a low computational cost.

Data Availability

The data used to support the findings of this study are currently under embargo, while the research findings are commercialized. Requests for data, 12 months after publication of this article, will be considered by the corresponding author.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research project was financially supported by the National Science and Technology Major Project (Grant number 2017ZX05009-005) and the National Natural Science Foundation of China (no. 51774255). The authors would also like to thank Lian Wang and XiLing Wang for sharing their meaningful points and discussions.