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

Authors: Liang Zhang; ZhiPing Li; Hong Li; Caspar Daniel Adenutsi; FengPeng Lai; KongJie Wang; Sen Yang
Journal: Mathematical Problems in Engineering
 (2020)

---

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


---

## Body

## 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 [4–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–15]. 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) [17–21] and proper orthogonal decomposition (POD) method [22–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]. There are also numerous applications of PCE in petroleum engineering problems [30–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 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 Section2.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, 40–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 two-phase flow in reservoir system can be given as(1)ϕ∂Sw∂t−∇⋅K⋅krwμw∇p−ρwg−qw=0,−ϕ∂So∂t−∇⋅K⋅kroμo∇p−ρog−qo=0,which is constrained by the equation(2)Sw+So=1.The subscriptsw and o stand for water phase and oil phase, respectively; Φ is the porosity; Sw and So are the water-phase saturation and oil saturation, respectively; K is the absolute permeability; krw and kro are the relative permeability of water and oil, respectively, which are functions of Sw and So; ρw and ρo are the water density and oil density, respectively; g is the gravity vector;t stands for time; and qw and qo 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:(3)Krw=KrwSorSwdnw,Kro=KroSwi1−Swdno,where(4)Swd=Sw−Swi1−Swi−Sor,where KrwSor and KroSwi denote the water and oil relative permeability, Sor is the residual oil saturation, and Swi is the connate water saturation. Furthermore,Sw is the water saturation, Swd is the normalized water saturation considering the irreducible saturations of both phases, and nw and no 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 ofSwi, Sor, KrwSor, KroSwi, nw, and no. In this paper, only nw and no 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 nw and no that satisfies the uniform distribution within the range of 2–4. In short, the uncertain reservoir parameters u=ϕ,K,nw,no are considered in this paper, in which ϕ and K are considered to obey the Gaussian distribution in a certain range, whereas nw and no obey the uniform distribution.Figure 1
Oil-water relative permeability curves with uncertainnw and no.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 (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. Therefore, injection-production 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,nw,no. The NPV function is computed as follows [44]:(5)NPVd,u=∑k=1Nt∑j=1Npcoqoj,kd,uΔtk−∑j=1Npcwpqwpj,kd,uΔtk−∑j=1Ntcwiqwij,kd,uΔtk1+DkΔtk,where qoj,k, qwpj,k, and qwij,k are the flow rates of oil, water produced, and water injected for well j at the output interval k, respectively, and Δtk represents the length (in days) of each of Nt time steps in the simulation. The oil price and the cost of water produced and injected are denoted by co, cwp, and cwi, respectively. Furthermore, D 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.Figure 2
Procedure for the construction of PCE of NPV.


#### 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:(6)NPV=∑i=0∞αiΨiξ,ξ=ξ1,ξ2,…,ξn,where αii=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=x1,x2,…,xn can be linearly transformed into random variables ξ=ξ1,ξ2,…,ξn with the mean value of 0 and variance of 1, as follows:(7)ξi=xi−ExiVarxi,i=1,2,…,n,where Exi is the expectation of xi and Varxi is the variance of xi. Conversely speaking, if ξi is known, xi can be obtained using the following equation:(8)xi=Varxiξi+Exi,i=1,2,…,n.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:(9)NPV=∑i=0p−1αiΨiξ,ξ=ξ1,ξ2,…,ξn,with p=n+d!/n!d!. Because the parameters x1,x2,…,xn are independent of each other, the parameters ξ1,ξ2,…,ξn obtained through linear transformation of equation (7) are also independent of each other. Therefore, Ψiξ can be decomposed into a product of several one-dimensional orthogonal polynomials as depicted in the following equation:(10)Ψiξ1,ξ2,…,ξn=∏k=1nφklkξk,i=0,1,…,p−1;lk=0,1,…,d;∑k=1nlk≤d.In order to determine the polynomial chaos expansion of NPV, it is necessary to determine the one-dimensional orthogonal polynomial∏k=1nφklkξkk=1,2,…,n;lk=0,1,…,d and polynomial chaos expansion coefficients αii=0,1,…,p−1.


#### 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.Table 1
Random variables and orthogonal polynomial bases.
Distribution
Density functionρω
Orth. polynomial
Interval
Gaussian
1
/
2
π
exp
−
1
/
2
ω
2
Hermite
−
∞
,
+
∞
Gamma
exp
−
ω
ω
α
−
1
1
/
Γ
α
Laguerre
0
,
+
∞
Beta
Γ
α
+
β
/
Γ
α
Γ
β
2
α
+
β
−
1
ω
+
1
α
−
1
1
−
ω
β
−
1
Jacobi
−
1,1
Uniform
ω
+
1
/
2
Legendre
−
1,1
Poisson
λ
ω
/
ω
!
e
−
λ
Charlier
0,1,2
,
…Based on the work of Oladyshkin and Nowak [35], a method for determining the base of one-dimensional orthogonal polynomial φiliξili=0,1,…,d is given when ξi is arbitrarily distributed.The k-order origin moment ofξi can be defined as follows:(11)uk=∫Ωξikpξidξi.The coefficients of orthogonal polynomialφijiξi of the order ji can be determined by solving the following matrix equation [35]:(12)μ0μ1⋯μjiμ1μ2⋯μji+1⋮⋮⋮⋮μji−1μji⋯μ2ji−100⋯1β0jiβ1ji⋮βji−1jiβjiji=00⋮01.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:(13)φijiξi2=∫Ωφijiξi2pξidξi.Therefore, the normalized orthogonal baseφ¯ijiξi of φijiξi is as follows:(14)φ¯ijiξi=1φijiξi∑k=0jiβkjiξik.


#### 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 ξ1,ξ2,…,ξN are determined, the input parameters x1,x2,…,xN can be obtained using equation (8). Then, the samples of input parameter are brought into reservoir numerical simulator to obtain samples NPV=NPV1,NPV2,…,NPVN. Once the samples ξ and NPV are determined, the polynomial chaos expansion coefficient α can be calculated using the following equation:(15)α^=argminα1N∑j=1Nfxξj−∑i=0P−1αiΨiξj.Let(16)Z=Ψ0ξ1Ψ1ξ1⋯ΨP−1ξ1Ψ0ξ2Ψ1ξ2⋯ΨP−1ξ2⋮⋮⋱⋮Ψ0ξNΨ1ξN⋯ΨP−1ξN.Then,α^ that makes equation (15) valid should be(17)α^=ZTZ−1ZTNPV.The input samples are extracted in turn according to the above order until the matrix becomes reversible. When the matrixZTZ 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 (ENPV), variance (VarNPV), 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:(18)ENPV=∫Ω∑i=0P−1αiΨiξpξdξ=α0.Because of the normality ofΨiξi=0,1,…,P−1 and the independence of ξii=1,2,…,n, the variance of NPV can be calculated as follows:(19)VarNPV=ENPV−ENPV2=ENPV2−ENPV2=∑i=1P−1αi2.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 ordern, 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]:(20)RMSD=∑i=1mNPVi−NPV^i2/m∑i=1mNPVi/m.


#### 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:(21)NPV=fd,u,where d=d1,d2,…,dnd is the set of injection-production parameters, nd is the number of injection-production parameters, u=u1,u2,…,unu is the set of uncertain parameters, and nu 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 D samples di=d1i,d2i,…,dndi,i=1,2,…,M, and internal circulation, which sampled uncertain parameters through LHS to obtain the uncertain parameter set U samples uji=u1ji,u2ji,…,unuji,j=1,2,…,N. According to equation (14), the set of injection-production parameters and uncertain parameters is transformed to ηi=η1i,η2i,…,ηndi,i=1,2,…,M, and ξji=ξ1ji,ξ2ji,…,ξnuji,j=1,2,…,N. The uncertain parameter set U and injection-production parameters di are incorporated into reservoir numerical simulator to obtain NPV samples NPVdi,U=NPV1i,NPV2i,…,NPVNi. According to the theory of PCE, NPVdi,U can be expanded as(22)NPVdi,U=∑k=0s−1ukiΨkξi.Once the random baseΨξi is selected, the coefficient uki will become a definite value. The value of uki depends only on ηi, due to which, uki can be expressed as ukηi. Based on response surface model (RSM), ukη can be approximated using the following equation:(23)ukη≈∑l=0t−1bk,lφlη=bkTφlη.In equation (23), φη=φ0η,φ1η,…,φt−1ηT is the base of RSM. bk=bk,0,bk,1,…,bk,t−1 are the coefficients of RSM. T represents the transpose operator. Changing equation (23) into equation (22), NPVd,U can be expressed as follows:(24)NPVd,U=∑k=0s−1ukiΨkξi=∑k=0s−1∑l=0t−1bk,lφlηΨkξi.The RSM can also be developed using LHS, and the number of selected sampling pointsdi is at least 2t [50]. Assuming that the number of sampling points is M, the coefficient vector of RSM bk can be obtained through a least square method as follows:(25)bk=FηTFη−1FηTukη,where(26)Fη=φ0η1φ1η1⋯φP−1η1φ0η2φ1η2⋯φP−1η2⋮⋮⋱⋮φ0ηMφ1ηM⋯φP−1ηM.This can be represented in the matrix form as(27)B=FηTFη−1FηTNPVD,UZξ,where(28)NPVD,U=NPVd1,u1NPVd1,u2⋯NPVd1,uNNPVd2,u2NPVd2,u2⋯NPVd2,uN⋮⋮⋱⋮NPVdM,u1NPVdM,u2⋯NPVdM,uN.


### 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:(29)maxxfx,us.t.Pgcx,u≤0≥αcc=1,2,…,C,x∈xmin,xmax,where P stands for the probability of conforming to constraint gcx,u≤0 and αc represents the c_th reliability index. The closer the αc is to 1, the higher the reliability of the system. In this paper, the object function fx,u 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:(30)Vi,dt+1=ωVi,dt+c1rand1pbesti,dt−Xi,dt+c2rand2gbestdt−Xi,dt,Xi,dt+1=Xi,dt+Vi,dt+1.When the threshold is exceeded, the particle performs the Levy flight operation. The positions of particles are then redistributed:(31)Xt+1=Xt+randomsizen⊕Levyβ.The basic flow chart of LFPSO algorithm is shown in Figure3. In Figure 3, triali 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,triali is reset to 0; otherwise, triali is added to 1.Figure 3
The flow chart of the LFPSO algorithm.


#### 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 Figure4. 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 algorithmFigure 4
Procedure for determining the optimal injection-production parameters under uncertainty.


## 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, 40–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 two-phase flow in reservoir system can be given as(1)ϕ∂Sw∂t−∇⋅K⋅krwμw∇p−ρwg−qw=0,−ϕ∂So∂t−∇⋅K⋅kroμo∇p−ρog−qo=0,which is constrained by the equation(2)Sw+So=1.The subscriptsw and o stand for water phase and oil phase, respectively; Φ is the porosity; Sw and So are the water-phase saturation and oil saturation, respectively; K is the absolute permeability; krw and kro are the relative permeability of water and oil, respectively, which are functions of Sw and So; ρw and ρo are the water density and oil density, respectively; g is the gravity vector;t stands for time; and qw and qo 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:(3)Krw=KrwSorSwdnw,Kro=KroSwi1−Swdno,where(4)Swd=Sw−Swi1−Swi−Sor,where KrwSor and KroSwi denote the water and oil relative permeability, Sor is the residual oil saturation, and Swi is the connate water saturation. Furthermore,Sw is the water saturation, Swd is the normalized water saturation considering the irreducible saturations of both phases, and nw and no 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 ofSwi, Sor, KrwSor, KroSwi, nw, and no. In this paper, only nw and no 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 nw and no that satisfies the uniform distribution within the range of 2–4. In short, the uncertain reservoir parameters u=ϕ,K,nw,no are considered in this paper, in which ϕ and K are considered to obey the Gaussian distribution in a certain range, whereas nw and no obey the uniform distribution.Figure 1
Oil-water relative permeability curves with uncertainnw and no.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 (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. Therefore, injection-production 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,nw,no. The NPV function is computed as follows [44]:(5)NPVd,u=∑k=1Nt∑j=1Npcoqoj,kd,uΔtk−∑j=1Npcwpqwpj,kd,uΔtk−∑j=1Ntcwiqwij,kd,uΔtk1+DkΔtk,where qoj,k, qwpj,k, and qwij,k are the flow rates of oil, water produced, and water injected for well j at the output interval k, respectively, and Δtk represents the length (in days) of each of Nt time steps in the simulation. The oil price and the cost of water produced and injected are denoted by co, cwp, and cwi, respectively. Furthermore, D 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.Figure 2
Procedure for the construction of PCE of NPV.


### 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:(6)NPV=∑i=0∞αiΨiξ,ξ=ξ1,ξ2,…,ξn,where αii=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=x1,x2,…,xn can be linearly transformed into random variables ξ=ξ1,ξ2,…,ξn with the mean value of 0 and variance of 1, as follows:(7)ξi=xi−ExiVarxi,i=1,2,…,n,where Exi is the expectation of xi and Varxi is the variance of xi. Conversely speaking, if ξi is known, xi can be obtained using the following equation:(8)xi=Varxiξi+Exi,i=1,2,…,n.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:(9)NPV=∑i=0p−1αiΨiξ,ξ=ξ1,ξ2,…,ξn,with p=n+d!/n!d!. Because the parameters x1,x2,…,xn are independent of each other, the parameters ξ1,ξ2,…,ξn obtained through linear transformation of equation (7) are also independent of each other. Therefore, Ψiξ can be decomposed into a product of several one-dimensional orthogonal polynomials as depicted in the following equation:(10)Ψiξ1,ξ2,…,ξn=∏k=1nφklkξk,i=0,1,…,p−1;lk=0,1,…,d;∑k=1nlk≤d.In order to determine the polynomial chaos expansion of NPV, it is necessary to determine the one-dimensional orthogonal polynomial∏k=1nφklkξkk=1,2,…,n;lk=0,1,…,d and polynomial chaos expansion coefficients αii=0,1,…,p−1.


### 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.Table 1
Random variables and orthogonal polynomial bases.
Distribution
Density functionρω
Orth. polynomial
Interval
Gaussian
1
/
2
π
exp
−
1
/
2
ω
2
Hermite
−
∞
,
+
∞
Gamma
exp
−
ω
ω
α
−
1
1
/
Γ
α
Laguerre
0
,
+
∞
Beta
Γ
α
+
β
/
Γ
α
Γ
β
2
α
+
β
−
1
ω
+
1
α
−
1
1
−
ω
β
−
1
Jacobi
−
1,1
Uniform
ω
+
1
/
2
Legendre
−
1,1
Poisson
λ
ω
/
ω
!
e
−
λ
Charlier
0,1,2
,
…Based on the work of Oladyshkin and Nowak [35], a method for determining the base of one-dimensional orthogonal polynomial φiliξili=0,1,…,d is given when ξi is arbitrarily distributed.The k-order origin moment ofξi can be defined as follows:(11)uk=∫Ωξikpξidξi.The coefficients of orthogonal polynomialφijiξi of the order ji can be determined by solving the following matrix equation [35]:(12)μ0μ1⋯μjiμ1μ2⋯μji+1⋮⋮⋮⋮μji−1μji⋯μ2ji−100⋯1β0jiβ1ji⋮βji−1jiβjiji=00⋮01.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:(13)φijiξi2=∫Ωφijiξi2pξidξi.Therefore, the normalized orthogonal baseφ¯ijiξi of φijiξi is as follows:(14)φ¯ijiξi=1φijiξi∑k=0jiβkjiξik.


### 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 ξ1,ξ2,…,ξN are determined, the input parameters x1,x2,…,xN can be obtained using equation (8). Then, the samples of input parameter are brought into reservoir numerical simulator to obtain samples NPV=NPV1,NPV2,…,NPVN. Once the samples ξ and NPV are determined, the polynomial chaos expansion coefficient α can be calculated using the following equation:(15)α^=argminα1N∑j=1Nfxξj−∑i=0P−1αiΨiξj.Let(16)Z=Ψ0ξ1Ψ1ξ1⋯ΨP−1ξ1Ψ0ξ2Ψ1ξ2⋯ΨP−1ξ2⋮⋮⋱⋮Ψ0ξNΨ1ξN⋯ΨP−1ξN.Then,α^ that makes equation (15) valid should be(17)α^=ZTZ−1ZTNPV.The input samples are extracted in turn according to the above order until the matrix becomes reversible. When the matrixZTZ 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 (ENPV), variance (VarNPV), 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:(18)ENPV=∫Ω∑i=0P−1αiΨiξpξdξ=α0.Because of the normality ofΨiξi=0,1,…,P−1 and the independence of ξii=1,2,…,n, the variance of NPV can be calculated as follows:(19)VarNPV=ENPV−ENPV2=ENPV2−ENPV2=∑i=1P−1αi2.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 ordern, 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]:(20)RMSD=∑i=1mNPVi−NPV^i2/m∑i=1mNPVi/m.


### 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:(21)NPV=fd,u,where d=d1,d2,…,dnd is the set of injection-production parameters, nd is the number of injection-production parameters, u=u1,u2,…,unu is the set of uncertain parameters, and nu 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 D samples di=d1i,d2i,…,dndi,i=1,2,…,M, and internal circulation, which sampled uncertain parameters through LHS to obtain the uncertain parameter set U samples uji=u1ji,u2ji,…,unuji,j=1,2,…,N. According to equation (14), the set of injection-production parameters and uncertain parameters is transformed to ηi=η1i,η2i,…,ηndi,i=1,2,…,M, and ξji=ξ1ji,ξ2ji,…,ξnuji,j=1,2,…,N. The uncertain parameter set U and injection-production parameters di are incorporated into reservoir numerical simulator to obtain NPV samples NPVdi,U=NPV1i,NPV2i,…,NPVNi. According to the theory of PCE, NPVdi,U can be expanded as(22)NPVdi,U=∑k=0s−1ukiΨkξi.Once the random baseΨξi is selected, the coefficient uki will become a definite value. The value of uki depends only on ηi, due to which, uki can be expressed as ukηi. Based on response surface model (RSM), ukη can be approximated using the following equation:(23)ukη≈∑l=0t−1bk,lφlη=bkTφlη.In equation (23), φη=φ0η,φ1η,…,φt−1ηT is the base of RSM. bk=bk,0,bk,1,…,bk,t−1 are the coefficients of RSM. T represents the transpose operator. Changing equation (23) into equation (22), NPVd,U can be expressed as follows:(24)NPVd,U=∑k=0s−1ukiΨkξi=∑k=0s−1∑l=0t−1bk,lφlηΨkξi.The RSM can also be developed using LHS, and the number of selected sampling pointsdi is at least 2t [50]. Assuming that the number of sampling points is M, the coefficient vector of RSM bk can be obtained through a least square method as follows:(25)bk=FηTFη−1FηTukη,where(26)Fη=φ0η1φ1η1⋯φP−1η1φ0η2φ1η2⋯φP−1η2⋮⋮⋱⋮φ0ηMφ1ηM⋯φP−1ηM.This can be represented in the matrix form as(27)B=FηTFη−1FηTNPVD,UZξ,where(28)NPVD,U=NPVd1,u1NPVd1,u2⋯NPVd1,uNNPVd2,u2NPVd2,u2⋯NPVd2,uN⋮⋮⋱⋮NPVdM,u1NPVdM,u2⋯NPVdM,uN.


## 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:(6)NPV=∑i=0∞αiΨiξ,ξ=ξ1,ξ2,…,ξn,where αii=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=x1,x2,…,xn can be linearly transformed into random variables ξ=ξ1,ξ2,…,ξn with the mean value of 0 and variance of 1, as follows:(7)ξi=xi−ExiVarxi,i=1,2,…,n,where Exi is the expectation of xi and Varxi is the variance of xi. Conversely speaking, if ξi is known, xi can be obtained using the following equation:(8)xi=Varxiξi+Exi,i=1,2,…,n.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:(9)NPV=∑i=0p−1αiΨiξ,ξ=ξ1,ξ2,…,ξn,with p=n+d!/n!d!. Because the parameters x1,x2,…,xn are independent of each other, the parameters ξ1,ξ2,…,ξn obtained through linear transformation of equation (7) are also independent of each other. Therefore, Ψiξ can be decomposed into a product of several one-dimensional orthogonal polynomials as depicted in the following equation:(10)Ψiξ1,ξ2,…,ξn=∏k=1nφklkξk,i=0,1,…,p−1;lk=0,1,…,d;∑k=1nlk≤d.In order to determine the polynomial chaos expansion of NPV, it is necessary to determine the one-dimensional orthogonal polynomial∏k=1nφklkξkk=1,2,…,n;lk=0,1,…,d and polynomial chaos expansion coefficients αii=0,1,…,p−1.


## 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.Table 1
Random variables and orthogonal polynomial bases.
Distribution
Density functionρω
Orth. polynomial
Interval
Gaussian
1
/
2
π
exp
−
1
/
2
ω
2
Hermite
−
∞
,
+
∞
Gamma
exp
−
ω
ω
α
−
1
1
/
Γ
α
Laguerre
0
,
+
∞
Beta
Γ
α
+
β
/
Γ
α
Γ
β
2
α
+
β
−
1
ω
+
1
α
−
1
1
−
ω
β
−
1
Jacobi
−
1,1
Uniform
ω
+
1
/
2
Legendre
−
1,1
Poisson
λ
ω
/
ω
!
e
−
λ
Charlier
0,1,2
,
…Based on the work of Oladyshkin and Nowak [35], a method for determining the base of one-dimensional orthogonal polynomial φiliξili=0,1,…,d is given when ξi is arbitrarily distributed.The k-order origin moment ofξi can be defined as follows:(11)uk=∫Ωξikpξidξi.The coefficients of orthogonal polynomialφijiξi of the order ji can be determined by solving the following matrix equation [35]:(12)μ0μ1⋯μjiμ1μ2⋯μji+1⋮⋮⋮⋮μji−1μji⋯μ2ji−100⋯1β0jiβ1ji⋮βji−1jiβjiji=00⋮01.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:(13)φijiξi2=∫Ωφijiξi2pξidξi.Therefore, the normalized orthogonal baseφ¯ijiξi of φijiξi is as follows:(14)φ¯ijiξi=1φijiξi∑k=0jiβkjiξik.


## 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 ξ1,ξ2,…,ξN are determined, the input parameters x1,x2,…,xN can be obtained using equation (8). Then, the samples of input parameter are brought into reservoir numerical simulator to obtain samples NPV=NPV1,NPV2,…,NPVN. Once the samples ξ and NPV are determined, the polynomial chaos expansion coefficient α can be calculated using the following equation:(15)α^=argminα1N∑j=1Nfxξj−∑i=0P−1αiΨiξj.Let(16)Z=Ψ0ξ1Ψ1ξ1⋯ΨP−1ξ1Ψ0ξ2Ψ1ξ2⋯ΨP−1ξ2⋮⋮⋱⋮Ψ0ξNΨ1ξN⋯ΨP−1ξN.Then,α^ that makes equation (15) valid should be(17)α^=ZTZ−1ZTNPV.The input samples are extracted in turn according to the above order until the matrix becomes reversible. When the matrixZTZ 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 (ENPV), variance (VarNPV), 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:(18)ENPV=∫Ω∑i=0P−1αiΨiξpξdξ=α0.Because of the normality ofΨiξi=0,1,…,P−1 and the independence of ξii=1,2,…,n, the variance of NPV can be calculated as follows:(19)VarNPV=ENPV−ENPV2=ENPV2−ENPV2=∑i=1P−1αi2.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 ordern, 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]:(20)RMSD=∑i=1mNPVi−NPV^i2/m∑i=1mNPVi/m.


## 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:(21)NPV=fd,u,where d=d1,d2,…,dnd is the set of injection-production parameters, nd is the number of injection-production parameters, u=u1,u2,…,unu is the set of uncertain parameters, and nu 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 D samples di=d1i,d2i,…,dndi,i=1,2,…,M, and internal circulation, which sampled uncertain parameters through LHS to obtain the uncertain parameter set U samples uji=u1ji,u2ji,…,unuji,j=1,2,…,N. According to equation (14), the set of injection-production parameters and uncertain parameters is transformed to ηi=η1i,η2i,…,ηndi,i=1,2,…,M, and ξji=ξ1ji,ξ2ji,…,ξnuji,j=1,2,…,N. The uncertain parameter set U and injection-production parameters di are incorporated into reservoir numerical simulator to obtain NPV samples NPVdi,U=NPV1i,NPV2i,…,NPVNi. According to the theory of PCE, NPVdi,U can be expanded as(22)NPVdi,U=∑k=0s−1ukiΨkξi.Once the random baseΨξi is selected, the coefficient uki will become a definite value. The value of uki depends only on ηi, due to which, uki can be expressed as ukηi. Based on response surface model (RSM), ukη can be approximated using the following equation:(23)ukη≈∑l=0t−1bk,lφlη=bkTφlη.In equation (23), φη=φ0η,φ1η,…,φt−1ηT is the base of RSM. bk=bk,0,bk,1,…,bk,t−1 are the coefficients of RSM. T represents the transpose operator. Changing equation (23) into equation (22), NPVd,U can be expressed as follows:(24)NPVd,U=∑k=0s−1ukiΨkξi=∑k=0s−1∑l=0t−1bk,lφlηΨkξi.The RSM can also be developed using LHS, and the number of selected sampling pointsdi is at least 2t [50]. Assuming that the number of sampling points is M, the coefficient vector of RSM bk can be obtained through a least square method as follows:(25)bk=FηTFη−1FηTukη,where(26)Fη=φ0η1φ1η1⋯φP−1η1φ0η2φ1η2⋯φP−1η2⋮⋮⋱⋮φ0ηMφ1ηM⋯φP−1ηM.This can be represented in the matrix form as(27)B=FηTFη−1FηTNPVD,UZξ,where(28)NPVD,U=NPVd1,u1NPVd1,u2⋯NPVd1,uNNPVd2,u2NPVd2,u2⋯NPVd2,uN⋮⋮⋱⋮NPVdM,u1NPVdM,u2⋯NPVdM,uN.


## 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:(29)maxxfx,us.t.Pgcx,u≤0≥αcc=1,2,…,C,x∈xmin,xmax,where P stands for the probability of conforming to constraint gcx,u≤0 and αc represents the c_th reliability index. The closer the αc is to 1, the higher the reliability of the system. In this paper, the object function fx,u 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:(30)Vi,dt+1=ωVi,dt+c1rand1pbesti,dt−Xi,dt+c2rand2gbestdt−Xi,dt,Xi,dt+1=Xi,dt+Vi,dt+1.When the threshold is exceeded, the particle performs the Levy flight operation. The positions of particles are then redistributed:(31)Xt+1=Xt+randomsizen⊕Levyβ.The basic flow chart of LFPSO algorithm is shown in Figure3. In Figure 3, triali 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,triali is reset to 0; otherwise, triali is added to 1.Figure 3
The flow chart of the LFPSO algorithm.


### 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 Figure4. 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 algorithmFigure 4
Procedure for determining the optimal injection-production parameters under uncertainty.


## 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:(30)Vi,dt+1=ωVi,dt+c1rand1pbesti,dt−Xi,dt+c2rand2gbestdt−Xi,dt,Xi,dt+1=Xi,dt+Vi,dt+1.When the threshold is exceeded, the particle performs the Levy flight operation. The positions of particles are then redistributed:(31)Xt+1=Xt+randomsizen⊕Levyβ.The basic flow chart of LFPSO algorithm is shown in Figure3. In Figure 3, triali 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,triali is reset to 0; otherwise, triali is added to 1.Figure 3
The flow chart of the LFPSO algorithm.


## 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 Figure4. 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 algorithmFigure 4
Procedure for determining the optimal injection-production parameters under uncertainty.


## 3. Case Study and Analysis

### 3.1. Example Description

According to the analysis provided in Section2, 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 Table2.Table 2
Distribution and range of uncertain parameters for the case.
Uncertain parameters
Distribution
Mean
Range
Porosity
Truncated Gaussian
15
10∼20
Permeability (mD)
Truncated Gaussian
75
50∼100
Water-phase index
Uniform
3
2∼4
Oil-phase index
Uniform
3
2∼4The values of parameters for computing NPV and the remaining reservoir system properties are presented in Table3.Table 3
System properties.
Properties
Value
Properties
Value
Price of oilCo
60$/STB
Water densityρw
1014
kg
⋅
m
3
Costs of water disposalCw
6$ (m3)
Oil densityρo
850
kg
⋅
m
3
The annual discount rated
0.05
Water compressibilitycw
10
−
10
Pa
−
1
Water viscosityμw
10
−
3
Pa
⋅
s
Oil compressibilityco
10
−
10
Pa
−
1
Oil viscosityμo
8
×
10
−
3
Pa
⋅
s
Rock compressibilitycr
1.8
×
10
−
10
Pa
−
1
Costs of water injectionCwi
2$ (m3)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 Table4.Table 4
Range of injection and production parameters for the case.
Optimized parameters
Range
Distance of oil and water wells (L)
50∼500 m
Injection rate (Q)
50∼500 m³/d
Bottom hole pressure (P)
5∼18 MPa


### 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 variabled=15,300,100, 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.Figure 5
CDFs of NPV associated with the three cases based on PC.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 Figure6. 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.Figure 6
Comparison of NPV estimated by the 3rd order PCE and calculated by reservoir numerical simulation.In this study, the reliability valueαc 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.Figure 7
NPV reliability optimization curve under uncertainty (LFPSO- PCE).Figure 8
Reliability optimization curve of injection-production parameters under uncertainty.The optimization results of LFPSO combined with PCE with response surface model are presented in Table5. 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. 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.Table 5
Optimization results of LFPSO combined with PCE and RSM.
Optimized parameters
Values
Distance of oil and water wells (L)
254.94 m
Injection velocity (Q)
358.70 m³/d
Bottom hole pressure (P)
12.75 MPaFigure 9
Diagram of the relationship between design parameters and NPV.
(a)
(b)
(c)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 Figure10. 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.Figure 10
Probability density of NPV with optimal injection-production parameters.


## 3.1. Example Description

According to the analysis provided in Section2, 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 Table2.Table 2
Distribution and range of uncertain parameters for the case.
Uncertain parameters
Distribution
Mean
Range
Porosity
Truncated Gaussian
15
10∼20
Permeability (mD)
Truncated Gaussian
75
50∼100
Water-phase index
Uniform
3
2∼4
Oil-phase index
Uniform
3
2∼4The values of parameters for computing NPV and the remaining reservoir system properties are presented in Table3.Table 3
System properties.
Properties
Value
Properties
Value
Price of oilCo
60$/STB
Water densityρw
1014
kg
⋅
m
3
Costs of water disposalCw
6$ (m3)
Oil densityρo
850
kg
⋅
m
3
The annual discount rated
0.05
Water compressibilitycw
10
−
10
Pa
−
1
Water viscosityμw
10
−
3
Pa
⋅
s
Oil compressibilityco
10
−
10
Pa
−
1
Oil viscosityμo
8
×
10
−
3
Pa
⋅
s
Rock compressibilitycr
1.8
×
10
−
10
Pa
−
1
Costs of water injectionCwi
2$ (m3)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 Table4.Table 4
Range of injection and production parameters for the case.
Optimized parameters
Range
Distance of oil and water wells (L)
50∼500 m
Injection rate (Q)
50∼500 m³/d
Bottom hole pressure (P)
5∼18 MPa


## 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 variabled=15,300,100, 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.Figure 5
CDFs of NPV associated with the three cases based on PC.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 Figure6. 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.Figure 6
Comparison of NPV estimated by the 3rd order PCE and calculated by reservoir numerical simulation.In this study, the reliability valueαc 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.Figure 7
NPV reliability optimization curve under uncertainty (LFPSO- PCE).Figure 8
Reliability optimization curve of injection-production parameters under uncertainty.The optimization results of LFPSO combined with PCE with response surface model are presented in Table5. 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. 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.Table 5
Optimization results of LFPSO combined with PCE and RSM.
Optimized parameters
Values
Distance of oil and water wells (L)
254.94 m
Injection velocity (Q)
358.70 m³/d
Bottom hole pressure (P)
12.75 MPaFigure 9
Diagram of the relationship between design parameters and NPV.
(a)
(b)
(c)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 Figure10. 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.Figure 10
Probability density of NPV with optimal injection-production parameters.


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



---

*Source: 5374523-2020-04-21.xml*