Application of Polynomial Chaos Expansion to Optimize Injection-Production Parameters under Uncertainty

School of Energy Resources, China University of Geosciences, Beijing 100083, China Beijing Key Laboratory of Unconventional Natural Gas Geological Evaluation and Development Engineering, Beijing 100083, China School of Mechanical Engineering, Chongqing University of Technology, Chongqing 400050, China Institute of Industrial Research, Council for Scientific and Industrial Research, P.O. Box LG 576, Legon-Accra, Ghana


Introduction
e 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]. e optimal field development strategy depends on reservoir geometry and physical parameters of rocks and fluids. e 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 [4][5][6][7]. 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 [8][9][10][11][12][13][14][15]. ese 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.
erefore, 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. e 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) [17][18][19][20][21] and proper orthogonal decomposition (POD) method [22][23][24] 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]. ere are also numerous applications of PCE in petroleum engineering problems [30][31][32][33][34][35][36]. 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 injectionproduction parameters under uncertainty. e 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.

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. e paper is structured as follows: in Section 2.1, the uncertain parameters of the reservoir as well as injectionproduction 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. is 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.

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,[40][41][42]. 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 twophase flow in reservoir system can be given as which is constrained by the equation e subscripts w and o stand for water phase and oil phase, respectively; Φ is the porosity; S w and S o are the water-phase saturation and oil saturation, respectively; K is the absolute permeability; k rw and k ro are the relative permeability of water and oil, respectively, which are functions of S w and S o ; ρ w and ρ o are the water density and oil density, respectively; g is the gravity vector; t stands for time; and q w and q o 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. erefore, 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. e most common type of such a relationship is the Corey model [43], which takes the following form: where where K rw (S or ) and K ro (S wi ) denote the water and oil relative permeability, S or is the residual oil saturation, and S wi is the connate water saturation. Furthermore, S w is the water saturation, S w d is the normalized water saturation considering the irreducible saturations of both phases, and n w and n o 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 S wi , S or , K rw (S or ), K ro (S wi ), n w , and n o . In this paper, only n w and n o 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. ey are obtained using the Corey model considering n w and n o that satisfies the uniform distribution within the range of 2-4. In short, the uncertain reservoir parameters u � ϕ, K, n w , n o are considered in this paper, in which ϕ and K are considered to obey the Gaussian distribution in a certain range, whereas n w and n o 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. e bottom hole pressure (P), 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. erefore, injectionproduction parameters d � P, Q, L { } 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 d � P, Q, L { } and uncertainty parameters u � ϕ, K, n w , n o . e NPV function is computed as follows [44]: where q j,k o , q j,k wp , and q j,k wi are the flow rates of oil, water produced, and water injected for well j at the output interval k, respectively, and Δt k represents the length (in days) of each of N t time steps in the simulation. e oil price and the cost of water produced and injected are denoted by c o , c wp , and c wi , respectively. Furthermore, D is the annual discount rate.

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. e 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. e procedure for the construction of surrogate model of NPV using PCE relates to the uncertainty in input parameters and is presented in Figure 2.

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 α i (i � 0, 1, 2, . . .) values are the unknown coefficients of the PCE; ξ � ξ 1 , ξ 2 , . . . , ξ n 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 x � x 1 , x 2 , . . . , x n can be linearly transformed into random variables ξ � ξ 1 , ξ 2 , . . . , ξ n with the mean value of 0 and variance of 1, as follows: where E[x i ] is the expectation of x i and Var(x i ) is the variance of x i . Conversely speaking, if ξ i is known, x i 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. e d-order truncated polynomial chaos expansion of NPV with respect to ndimensional random variable ξ is rewritten in the form, as shown in the following equation: with p � ((n + d)!/n!d!). Because the parameters x 1 , x 2 , . . . , x n are independent of each other, the parameters ξ 1 , ξ 2 , . . . , ξ n obtained through linear transformation of equation (7) are also independent of each other. erefore, Ψ i (ξ) can be decomposed into a product of several onedimensional 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 n k�1 φ l k k (ξ k )(k � 1, 2, . . . , n; l k � 0, 1, . . . , d) and polynomial chaos expansion coefficients α i (i � 0, 1, . . . , p − 1).

Selection of Orthogonal Polynomial Chaos Bases.
e 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 φ e k-order origin moment of ξ i can be defined as follows:  Mathematical Problems in Engineering e coefficients of orthogonal polynomial φ j i i (ξ i ) of the order j i 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 2d − 1-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: erefore, the normalized orthogonal base φ

Determination of Unknown Coefficients of PCE.
e 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 ξ (1) , ξ (2) , . . . , ξ (N) are determined, the input parameters x (1) , x (2) , . . . , x (N) can be obtained using equation (8).
en, the samples of input parameter are brought into reservoir numerical simulator to obtain samples NPV � NPV (1) , NPV (2) , . . . , NPV (N) . Once the samples ξ and NPV are determined, the polynomial chaos expansion coefficient α can be calculated using the following equation: en, α that makes equation (15) valid should be e input samples are extracted in turn according to the above order until the matrix becomes reversible. When the matrix Z T Z is ill-conditioned, the singular value decomposition (SVD) algorithm [48] can be used to determine α.

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. e 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 (E[NPV]), variance (Var[NPV]), 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 Ψ i (ξ)(i � 0, 1, . . . , P − 1) and the independence of ξ i (i � 1, 2, . . . , n), 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. erefore, 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.
e 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. e RMSD is calculated using the following equation [44]:

PCE Surrogate
Model. e NPV of oil field development is affected by both the uncertain parameters and injection-production parameters. Equation (5) is rewritten as follows: where d � d 1 , d 2 , . . . , d n d is the set of injection-production parameters, n d is the number of injection-production parameters, u � u 1 , u 2 , . . . , u n u is the set of uncertain parameters, and n u is the number of uncertain parameters. e 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 D samples Once the random base Ψ(ξ i ) is selected, the coefficient u i k will become a definite value. e value of u i k depends only on η i , due to which, u i k can be expressed as u k (η i ). Based on response surface model (RSM), u k (η) can be approximated using the following equation: In equation (23), 0 , b k,1 , . . . , b k,t− 1 ] are the coefficients of RSM. T represents the transpose operator. Changing equation (23) into equation (22), NPV(d, U) can be expressed as follows: e RSM can also be developed using LHS, and the number of selected sampling points d i is at least 2t [50]. Assuming that the number of sampling points is M, the coefficient vector of RSM b k can be obtained through a least square method as follows: where is can be represented in the matrix form as where

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. e optimization model is as follows: max where P stands for the probability of conforming to constraint g c (x, u) ≤ 0 and α c represents the c th reliability index. e closer the α c is to 1, the higher the reliability of the system. In this paper, the object function f(x, u) is the surrogate model of PCE with response surface model expressed as equation (24).

LFPSO.
LFPSO is used to solve the optimization problems as equation (24) to obtain the optimal injectionproduction 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. e method ensures that PSO can search globally and more effectively without falling into local stagnation. e 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. e positions of particles are then redistributed: e basic flow chart of LFPSO algorithm is shown in Figure 3. In Figure 3, trial (i) denotes the current number of trials and limit denotes the threshold set for each particle. It is worth noticing that if the updated particles are better,trial (i) is reset to 0; otherwise, trial (i) is added to 1.

Steps and Procedure of Optimization Design.
e 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. e procedure for determining the optimal injection-production parameters under uncertainty is shown in Figure 4. e 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 requirements Step 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 combination Step 3. To construct the PCE model of NPV and uncertain variables Step 4. To inspect the accuracy of PCE model and rearrange the test design until the accuracy requirements are met Step 5. To construct the PCE with RSM model of NPV and design variables on the basis of the PCE model satisfying the accuracy test Step 6. To analyze the uncertainty of the surrogate model and optimize the design variables based on the surrogate model and LFPSO algorithm

Example Description.
According to the analysis provided in Section 2, oil-water well distance, injection rate, and bottom hole pressure were the design parameters. ese 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 fivespot 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. e reservoir contained two-phase fluid of oil and water. e 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.
e distribution and range of uncertain reservoir parameters used in this example are presented in Table 2.
e 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.

Results and Discussion.
e order of PCE model should be determined when the PCE is established using internal circulation sampling. Taking the external design variable d � 15, 300, 100 { }, the 2nd, 3 rd , and 4th order PCE were constructed by using 21, 56, and 126 samples (respectively) to determine the surrogate model of NPV in this example.
e CDFs of NPV, which were estimated by coupling the 2nd, 3rd, and 4th order PCE with LHS, are shown in Figure 5. e 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. e results showed that the surrogate model of the 3rd order PCE with response surface model is considerably accurate. erefore, 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 α c of equation (31) is 0.97. e 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.
e 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 (P), injection rate (Q), and distance between the injection and production wells (L), were taken from Table 5. e 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. e 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. is showed that, under the influence of uncertain    Figure 4: Procedure for determining the optimal injection-production parameters under uncertainty.  parameters, the proposed method can be effectively applied to design optimization of injection-production parameters in waterflooding.

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. e 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
e 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.