Nonlinear risk optimization approach to water drive gas reservoir production optimization using DOE and artificial intelligence

Nonlinear risk optimization approach to water drive gas reservoir production optimization using DOE and artificial intelligence

Accepted Manuscript Nonlinear Risk Optimization Approach to Water Drive Gas Reservoir Production Optimization Using DOE and Artificial Intelligence Me...

2MB Sizes 0 Downloads 50 Views

Accepted Manuscript Nonlinear Risk Optimization Approach to Water Drive Gas Reservoir Production Optimization Using DOE and Artificial Intelligence Meysam Naderi, Ehsan Khamehchi PII:

S1875-5100(16)30176-7

DOI:

10.1016/j.jngse.2016.03.069

Reference:

JNGSE 1381

To appear in:

Journal of Natural Gas Science and Engineering

Received Date: 25 December 2015 Revised Date:

21 March 2016

Accepted Date: 22 March 2016

Please cite this article as: Naderi, M., Khamehchi, E., Nonlinear Risk Optimization Approach to Water Drive Gas Reservoir Production Optimization Using DOE and Artificial Intelligence, Journal of Natural Gas Science & Engineering (2016), doi: 10.1016/j.jngse.2016.03.069. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Nonlinear Risk Optimization Approach to Water Drive Gas Reservoir Production Optimization Using DOE and Artificial Intelligence Meysam Naderi, Ehsan Khamehchi∗ Faculty of Petroleum Engineering, Amirkabir University of Technology (Tehran Polytechnic), Hafez Avenue, Tehran, Iran

RI PT

Highlights

A new method of uncertainty reduction for water drive gas reservoirs was proposed. Design of experiments and relative variation factor were used in this work.

The method combines conventional uncertainty assessment and genetic algorithm.

SC

The proposed model has minimum relative error among different methods.

M AN U

Abstract

The risk of investment in petroleum industry is usually high due to uncertainties associated with performance prediction of hydrocarbon reservoir. To reduce the investment risk, the uncertainties should be estimated accurately. Statistical methods of risk and uncertainty assessment in oil and gas industry have limitations because of underlying principles.

TE D

This study proposes an improved strategy to reduce uncertainty associated with different engineering and geologic factors. The method optimizes conventional uncertainty assessment methods using evolutionary computing algorithm. In this regard, first; gas recovery factor from water drive gas reservoirs is predicted using five

EP

different methods including Box-Behnken design (BBD), Full Factorial design (FFD), Central Composite design (CCD), Uniform design (UD) and relative variation factor

AC C

(RVF). Next, genetic algorithm (GA) is applied to optimally integrate all conventional uncertainty assessment methods in one equation as a new approach. Part of water drive gas reservoir based on actual data is simulated to extract necessary information in order to be able to quantify the uncertainty of ultimate recovery factor. The effect of average reservoir permeability (KG), permeability anisotropy (Kv/Kh), aquifer size (Vaq), production rate (Qg), perforated thickness (Hp) and lower limit of tubing head pressure (THP) on reservoir performance has been studied.



Corresponding author. Tel.: 98 21 6454 5154 E-mail addresses: [email protected]

Page 1 of 24

ACCEPTED MANUSCRIPT The results show that average reservoir permeability, tubing head pressure, permeability anisotropy, aquifer size, perforated thickness, and production rate have the greatest impact on recovery factor uncertainty, respectively. The results also indicate the higher performance of proposed method over the conventional uncertainty assessment methods to predict the most likely recovery

RI PT

factor. The prediction error of most likely recovery factor by applying GA, BBD, FFD, CCD, UD, and RVF is 0.11, 2.2, 3.15, 2.78, 7.22, and 11.4% respectively.

Keywords: Box-Behnken Design, Full Factorial Design, Central Composite Design,

SC

Uniform Design, Relative Variation Factor, Genetic Algorithm 1) Introduction

Gas recovery from water drive reservoirs depends on many uncertain factors. The

M AN U

uncertainty associated with engineering, geologic, and economic variables influences the risk of reservoir development plans. Therefore, it is required to study the effect of various uncertain factors on reservoir performance in order to reduce the risk of exploration and exploitation plans. In this regard, uncertainty assessment techniques provide useful information for investors to invest with minimum risk.

TE D

A literature review indicates that several authors have investigated the performance of gas reservoirs with an active aquifer. The experimental study of residual gas saturations under water drive showed that prediction of recovery factor could be uncertain because of large trapped gas volume at abandonment pressure (Geffen et al., 1952). Performance

EP

evaluation of real water drive fields showed that the errors in the past performance data and aquifer characteristics cause uncertain predictions (Simons and Spain, 1958;

AC C

Chierici and Ciucci, 1967; Zee Ma, 2010). The simulation study by Cronquist (1973) revealed that areal permeability variation, production rate and well pattern influences recovery efficiency from partial water drive gas reservoirs. McVay and Dossary (2014) studied the value of assessing uncertainty, and the monetary impact of overconfidence on portfolio performance. Also different uncertainty quantification techniques have been used to estimate uncertainty in predictive reservoir simulations including gradient-based techniques (Lepine et al., 1999), maximum likelihood model (Subbey et al., 2004), response surface method (Fetel and Caumon, 2008), Markov Chain Monte Carlo (Liu and McVay, 2010), iterative procedures (Maschio and Schiozer, 2013) and desirability concept (Naderi et al., 2014). Page 2 of 24

ACCEPTED MANUSCRIPT In order to achieve better prediction accuracy, this paper aims to introduce a new useful methodology to reduce uncertainty associated with recovery factor from water drive gas reservoirs. In this regard, we applied design of experiment (DOE) techniques and relative variation factor (RVF) for uncertainty assessment. We used Box-Behnken design (BBD), Full Factorial design (FFD), Central Composite design (CCD) and Uniform

RI PT

design (UD) in DOE section. To develop the recovery factor response function, the simultaneous effect of six uncertain factors including average reservoir permeability (KG), permeability anisotropy (Kv/Kh), aquifer size (Vaq), production rate (Qg),

of gas recovery factor has been studied.

SC

perforated thickness (Hp) and lower limit of tubing head pressure (THP) on distribution

Eventually, genetic algorithm (GA) was applied to combine all five methods optimally in

AC C

EP

TE D

M AN U

order to improve the prediction performance.

Page 3 of 24

ACCEPTED MANUSCRIPT 2) Model description The purpose of this paper is to present a methodology to reduce the uncertainty associated with performance prediction of water drive gas reservoirs. In this regard, we used three dimensional and heterogeneous reservoirs with the well located at the center. The aquifer (in blue) is attached to the bottom of reservoir (in red). pressure. Fig. 1 shows the numerical simulated model.

RI PT

The well is capable of producing either at constant gas rate or at constant bottomhole

This model totally consists of 40 layers in the vertical direction. The reservoir has 30 layers with total thickness of 300 ft. The aquifer has 10 layers. The first 5 layers of

SC

aquifer have thickness of 50 ft. However, the other 5 layers of aquifer have variable thickness in order to be able to define the aquifer strength. The ratio of the aquifer pore volume to the gas pore volume is called aquifer strength or aquifer size. Aquifer

M AN U

strength determines the amount of reservoir energy that can be provided by water drive mechanism. The areal dimension of reservoir and aquifer is 3000 ft × 3000 ft and 6000 ft × 6000 ft respectively.

The flow is simulated using a gridding scheme that is locally refined around the well and coarsened away from the well to allow accurate modeling of the pressure drop in the

TE D

near wellbore region. The grid refinement is necessary because a large gridblock with dimension of about 165 ft cannot resolve near wellbore effects and consequently overestimate well productivities (Malachowski et al., 1995). Dykstra-Parson (1950) coefficient was used to characterize the heterogeneity level of % −  .% %

(1)

AC C

 =

EP

model. The Dykstra-Parsons permeability variation is defined by the following equation:

Where K50% and K84.1% are corresponding permeability values at 50% and 84.1% of cumulative sample thickness respectively. We used a heterogeneous model with Dykstra-Parson coefficient of 0.66. Reservoir and aquifer porosity is constant in all simulation models and it was set to 25%. The initial reservoir pressure and temperature is 5290 psi and 220°F at datum depth of 8530 ft. For simulation models, we used a vertical well with 7-inch internal – tubing diameter. Gas is produced through the vertical well at tubing head pressure control mode with two production constraints. The constraints are economic production rate and bottom

Page 4 of 24

ACCEPTED MANUSCRIPT hole pressure. The economic production rate is set to 10% of maximum initial gas production rate, and minimum bottom hole pressure is 500 psi in all simulations. The correlation developed by Lee et al. (1966) was used to predict reservoir gas viscosity. In addition, gas deviation factor was estimated using Dranchuk et al. (1974) correlations.

RI PT

The relative permeability data and capillary pressure is identical to the data sets measured for gas – water systems in the experimental work of Chierici et al. (1963). Table 1 and 2 presents fluid properties and the relative permeability data respectively. The Petalas and Aziz (2000) mechanistic model is used for calculating vertical flow

SC

performance curves.

It is noteworthy that the original gas in place and recovery factor of base case model is 1.44 × 108 (MSCF) and 75.46% respectively. Six factors of geometric average of

M AN U

reservoir permeability, permeability anisotropy, aquifer size, production rate, perforated thickness, and tubing head pressure were selected due to their expected impact on water drive gas reservoirs performance for uncertainty analysis. The following section describes five conventional methods of uncertainty assessment and

3) Methodology

TE D

our approach by applying genetic algorithm.

There are various statistical methods of uncertainty analysis, which have been used in the petroleum industry for many years. Monte Carlo simulation (MCS), design of

EP

experiment and response surface methodology, multiple realization tree and relative variation factor are typical methods of quantifying uncertainties with production

AC C

forecast.

Design of experiment and response surface methodology have been used in various fields of petroleum engineering including screening and sensitivity studies in recoverable reserves (White et al., 2001; Corre et al., 2000), production forecasting and estimating ultimate recovery curves (Damsleth et al., 1992; Friedmann et al., 2003; Venkataraman, 2000) and field development optimization (Aanonsen et al., 1995; Dejean and Blanc, 1999). Relative variation factor has been used to assess recovery factor and reserve uncertainty (Steagall and Schiozer, 2001; Salomao and Grell, 2001). The following section introduces design of experiment and response surface

Page 5 of 24

ACCEPTED MANUSCRIPT methodology application for uncertainty analysis of water drive gas reservoirs using BBD, FFD, CCD, UD and MCS. 3.1)

Design of experiment and response surface methodology

Design of experiment is a useful method to obtain maximum information at the minimum experimental cost by varying all the uncertain parameters simultaneously.

for empirical model building.

RI PT

Response surface methodology is a collection of mathematical and statistical techniques

Uncertainty assessment by using DOE and RSM consists of the following steps: 1) Build a reservoir simulation model

SC

2) Perform sensitivity study and determine the most effective uncertain factors 3) Design simulation tests for various levels of factors based on DOE 4) Run simulator for designed models and extract all necessary information

M AN U

5) Apply RSM to develop response functions

6) Assign appropriate distribution function for uncertain factors based on available information, statistical principles, expert opinion and historical data 7) Determine final distribution of response by applying Monte Carlo simulation The first step of uncertainty assessment using DOE and RSM is building a simulation

TE D

model of reservoir. We used a synthetic reservoir model in current study. However, actual data has been used in the construction of simulated models. The second step is determination of the most effective factors. Usually, there are many number of reservoir factors that affect the recovery from water drive gas reservoirs. Because of expensive

EP

and time-consuming nature of reservoir simulation, and rapid increasing of required simulation tests with the number of parameters, only the most effective factors should

AC C

be selected for uncertainty analysis. The most influential reservoir parameters can be determined from sensitivity analysis. In this regard, we used tornado chart of recovery factor. To construct the tornado diagram of recovery factor for n reservoir parameters, it is required to perform 2n+1 reservoir simulation runs: one run for all the reservoir parameters at their most likely value and two runs for each parameter; one at the pessimistic value and the other at the optimistic value for each parameter. Then, the most-likely, pessimistic, and optimistic recovery factor values are calculated with the corresponding reservoir parameter value. For example, for average reservoir permeability, the optimistic recovery factor is calculated with the optimistic average reservoir permeability value; the most likely Page 6 of 24

ACCEPTED MANUSCRIPT recovery factor is calculated with the most likely average reservoir permeability value, and the pessimistic recovery factor value is calculated with the pessimistic value. Finally, a recovery factor range is obtained from the pessimistic and optimistic recovery factor values, and a tornado diagram is formed. Fig. 2 shows the tornado diagram of recovery factor for six factors. As can be seen from

RI PT

tornado diagram, average reservoir permeability, tubing head pressure, permeability anisotropy, aquifer size, perforated thickness, and production rate have the greatest impact on recovery respectively. Therefore, we selected four factors of average reservoir permeability, tubing head pressure, permeability anisotropy and aquifer size

SC

for uncertainty analysis based on their relative importance.

We used four designs of BBD, FFD, CCD, and UD to perform simulation prediction tests. Based on BBD, FFD, CCD, and UD, it is required to design and simulate 25, 81, 25, and 30

M AN U

models respectively. Therefore, the total number of required simulations is 161. A full factorial design is one type of DOE that allows the experimenter to study the effects of all possible combinations of the several factors levels on the response including main and interaction effects. Also by using full factorial design, it is possible to test the model curvature (Box, 1980).

TE D

A central composite design is one of the main types of response surface designs that contains a factorial or fractional factorial design with center points, augmented with a group of axial points. This design is capable to estimate the main and square power effect of factors with curvature efficiently. Circumscribed, inscribed and face centered

EP

are three different types of central composite designs (Box-Wilson, 1951). A Box-Behnken design is another type of response surface design that does not contain

AC C

an embedded factorial or fractional factorial design. This design like central composite designs allows efficient estimation of first and second order terms in the model. Relative to central composite designs, Box-Behnken designs usually have fewer design points. Therefore, for the same number of factors performing Box-Behnken design can be less time consuming and expensive (Box- Behnken, 1960). The uniform design is one type of space filling designs that seeks its design points to be uniformly scattered on the domain of experiment. This design like other types of designs can be used for situations when the underlying model is unknown (Fang, 1980). Table 3 shows the range of variables examined in this survey with their appropriate transfer functions. In RSM, transfer functions are used to normalize range of data Page 7 of 24

ACCEPTED MANUSCRIPT between -1 and 1, and also ease of computation. In this table, X1, X2, X3 and X4 are coded variables between -1 and 1 for KG, Kv/Kh, Vaq and THP respectively. For example, range of average reservoir permeability (KG) from (1, 10, 100) can be normalized to (-1, 0, 1) by using transfer function of Log (KG)-1. It is worth to note that the production rate is set to 90 MMSCFD and perforated thickness is 60% of total gas column thickness in all

RI PT

simulations. The forth step of uncertainty analysis is to perform flow simulations based on DOE.

Table 4 only shows the result of performing flow simulation based on FFD. It is noteworthy that all flow simulations were undertaken using the Eclipse® 100 black oil

SC

simulator.

The fifth step is determination of response function using RSM. We analyzed the simulator output using response surface methodology. There are three main

M AN U

approaches to develop response functions including forward selection, backward elimination, and bidirectional elimination (Efroymson, 1960). In this study, we used backward elimination method, and the following equations developed for gas recovery factor (RF) based on BBD, FFD, CCD and UD:

 = 75.4 + 15.08 − 1.28 + 1.45 − 8.77 − 10.35 − 2.79

(2)

(3)

!! = 74.9 + 15.57 + 0.81 + 1.96 − 8.08 − 10.34 − 3.86 + 3.15 − 1.62 

(4)

EP

TE D

 = 74.1 + 15.6 + 0.0065  + 1.39 − 8.43 − 9.57 − 2.01 + 1.01 − 0.764   − 1.81  + 0.753 

(5)

AC C

" = 73.6 + 8.89  − 12.4  − 7.08  − 5.33  − 2.78  + 2.4   + 5.66   + 5.29   + 1.96  

Where X1, X2, X3 and X4 are coded variables for mentioned factors in Table 3. We have different responses for recovery factor because of different combinations of design points based on BBD, FFD, CCD, and UD. We have used 25, 81, 25, and 30 simulations for developing the equations based on BBD, FFD, CCD, and UD respectively. Table 5 shows the coefficient of determination (R2) and adjusted coefficient of determination (R2adj) for derived response functions. We used R2 and R2adj to show the goodness of fit for the regressions. It is worth to note that R2 explains the percentage of variability in the process by the fitted model. The R2adj values show the goodness of fit for the regressions adjusted for number of terms. The low difference between R2adj and Page 8 of 24

ACCEPTED MANUSCRIPT R2 indicates that the necessary terms have been included in the model (Myers and Montgomery, 1995). Fig. 3 shows the validity of proxy equations with a plot of predicted versus actual experiments. In addition, the validity of developed equations has been tested by performing 50 blind simulation tests outside the training group (Table 7).

RI PT

The significance of derived response functions for recovery factor was tested by analysis of variance. In addition, the regressions so obtained were analyzed for relative importance of the factors main effects, quadratic and interaction effects. The result of analysis of variance shows that all terms included in the above equations are significant

SC

at the 10% significance level for prediction at 90% confidence interval.

Finally, the cumulative probability distribution of recovery factor can be obtained by

performing Monte Carlo simulation. 3.2)

M AN U

assigning appropriate probability distribution function for uncertain factors and

Relative variation factor method

This method of uncertainty assessment consists of the following steps: 1) Build a reservoir simulation model

TE D

2) Select the most effective uncertain factors using sensitivity analysis 3) Calculate the most likely RF by setting all factors at their most likely values 4) Determine corresponding pessimistic and optimistic RF for all factors 5) Divide the resultant three RF by the most likely RF to obtain three real numbers

EP

of: less than one, one and greater than one for each factor 6) Assign triangular probability density function for factors based on step 5

AC C

7) Perform Monte Carlo simulation to obtain final distribution of RF The first and second step of uncertainty analysis using relative variation factor is same as DOE and RSM method.

The third step is to obtain the most likely recovery factor. The most likely recovery factor is determined by setting all factors at their most likely values and performing flow simulation. The resultant recovery factor is called a base case value. The pessimistic and optimistic recovery factor for individual factors can be obtained using flow simulation while the other factors are kept constant at their most likely values. The next step is dividing pessimistic, most likely and optimistic recovery factor by the base case recovery factor to obtain three real numbers of less than one, one, and Page 9 of 24

ACCEPTED MANUSCRIPT greater than one. The resultant three real numbers are used for construction of triangular distribution function. These three numbers are the lower, most likely and upper bound of triangular distribution function. Finally, the distribution of recovery factor can be obtained using Monte Carlo simulation and following equation: (6)

RI PT

RFRVF = RF ML × TPDF X1 × TPDF X2 × TPDF X3 × TPDF X4

Where RFML is the most likely recovery factor, TPDFXi is the triangular probability density function for uncertain factors of Xi, and RFRVF is the final distribution of recovery factor by RVF method. Table 6 shows the calculated lower and upper bound of

SC

triangular distribution for average reservoir permeability, permeability anisotropy,

3.3)

Genetic algorithm

M AN U

aquifer size and tubing head pressure.

Genetic algorithm (GA) is a heuristic search and an effective optimization technique that mimics the process of natural selection (Holland, 1975). Genetic algorithms have been successfully applied to a wide range of complex problems in the petroleum and natural gas industry. Genetic algorithms were used for optimum design of gas transmission

TE D

lines (Goldberg, 1983), reservoir characterization and modeling (Guerreiro et al., 1997 and Sen et al., 1992), petrophysics and petroleum geology (Fang et al., 1992), well test analysis (Yin et al., 1998), hydraulic fracturing design (Mohaghegh et al., 1999), optimum determination of all effective infiltration parameters in furrow irrigation

EP

(Valipour and Montazar, 2012), and sensitivity analysis of optimized infiltration parameters in SWDC model (Valipour and Montazar, 2012). In addition, GA was applied

AC C

for sand production onset (Khamehchi et al., 2014), gas-lift performance prediction and long-term gas-lift allocation optimization (Rasouli et al., 2015), history matching (Firoozjaee and Khamehchi, 2015), and pressure drop estimation in vertical multiphase flow (Ebrahimi and Khamehchi, 2015). The optimization process by using genetic algorithm might be as follow: 1) Generate an initial population of chromosomes randomly 2) Evaluate fitness of all chromosomes in the population using objective function 3) Select chromosomes with highest fitness to produce the next generation 4) Create a new population of chromosomes by using genetic operators such as crossover and mutation Page 10 of 24

ACCEPTED MANUSCRIPT 5) Replace the initial population with the new population 6) Continue the process until the stopping criteria have been met In genetic algorithm, chromosomes are string representations of solutions, a population is a set of chromosomes, and fitness function is a representation of a particular problem. More details about genetic algorithm optimization are given in Gen and Cheng (1997)

RI PT

and Hajela (1990). In this study, we applied GA to combine all five methods optimally in an integrated equation as a new approach for uncertainty assessment. We used GA for optimization instead of conventional methods (like least square method) because GA searches large

SC

numbers of candidate solutions simultaneously. In addition, the main advantage is that the GA gives the best set of coefficients to predict the recovery factor, and the result of least square method is one of the candidate solutions by GA.

M AN U

The final distribution of recovery factor is determined using following equation: NRFGA = A (NRFBBD) + B (NRFFFD) + C (NRFCCD) + D (NRFUD) + E (NRFRVF)

(7)

Where NRFBBD, NRFFFD, NRFCCD, NRFUD and NRFRVF are the distributions of normalized recovery factor using BBD, FFD, CCD, UD and RVF respectively. NRFGA is the final distribution of normalized recovery factor obtained by GA optimization. Also A, B, C, D,

new approach.

TE D

and E are optimal coefficients of each method. Fig. 4 shows the schematic diagram of

In order to obtain the optimal coefficients of each method, we used mean squared error (MSE) as the fitness function that should be minimized by GA. The MSE can be 0

#$% = &

(()* − (+,- ) /

(8)

AC C

12

EP

calculated using the following equation:

Where n is the number of new blind simulation tests, NRFSIM is the normalized actual recovery factor by simulation. We performed 50 new blind simulation tests to calculate MSE for GA optimization. Table 7 shows the result of these new blind tests. Parameter settings for GA optimization is given in Table 8. In this table, initial population size is 75 which specifies the number of individuals in each generation, and the initial range is [-2, 2] which specifies the range of vectors in the initial population. The crossover function is scattered with fraction of one. Mutation function is adaptive feasible. After running the GA, we applied optimized coefficients in equation 7 to produce the final distribution of recovery factor. Page 11 of 24

ACCEPTED MANUSCRIPT 4) Results and discussion The performance of all methods of uncertainty assessment illustrated in Table 9 shows that the GA optimized, BBD, FFD, CCD, RVF and UD have the minimum mean squared error respectively. In addition, GA optimized has maximum coefficient of determination among different methods. Therefore, GA optimized method was used for estimation of

RI PT

recovery factor. As can be seen from Fig. 5, the GA has reached to the best fitness value of 0.0038456 for normalized MSE after 51 generations. In addition, Fig. 5 shows best, worst and mean scores within 51 generations. The GA optimized coefficients for NRFBBD, NRFFFD, NRFCCD,

SC

NRFUD and NRFRVF are 1.598, -0.925, 0.327, 0.124 and -0.114 respectively. The estimated coefficients by GA yield the best fit with minimum prediction error. The recovery factor

M AN U

can be calculated using the following equation:

NRFGA = 1.598 × NRFBBD – 0.925 × NRFFFD + 0.327 × NRFCCD + 0.124 × NRFUD – 0.114 × NRFRVF

(9)

The final distribution of recovery factor is determined by applying Monte Carlo simulation. It is required to assign appropriate probability distribution function for all independent variables prior to performing Monte Carlo simulation. Therefore, based on statistical and previous studies (Venkataraman, 2000;Peake et al., 2005; Naderi et al.,

TE D

2014; Naderi et al., 2015), we assigned log-normal distribution for three factors of average reservoir permeability, permeability anisotropy and aquifer size; and triangular distribution for tubing head pressure.

Fig. 6 shows the result of Monte Carlo simulation by considering Latin Hypercube

EP

Sampling (LHS) and generating 1000 random numbers from entire range of factors. As can be seen, the relative variation factor and uniform design generate the smallest

AC C

uncertainty range (P90%-P10%) respectively while full factorial design predicts largest range for recovery factor. We summarized the findings from Fig. 6 in Table 10. Table 10 shows the minimum, P10%, P50%, P90%, maximum value of recovery factor predicted using six different methods, and also the prediction error of methods for base case (most likely) recovery factor. The most likely recovery factor is 75.46%, which was determined by setting all factors at their most likely values (middle levels) and performing flow simulation. As can be seen, GA optimized method has provided the smallest prediction error (0.11%). The relative error of BBD, CCD, FFD, UD, and RVF is 2.2, 2.78, 3.15, 7.22, and 11.4% respectively. Applying GA method for prediction of recovery factor has provided minimum MSE, maximum R2 and smallest relative error. Page 12 of 24

ACCEPTED MANUSCRIPT This indicates that GA optimized method performs better than individual conventional methods for prediction of recovery factor. In order to investigate the effect of different probability distributions on the final distribution of recovery factor; we considered two types of triangular and uniform distributions for all factors.

RI PT

Table 11 shows the effect of triangular and uniform distribution on final distribution of RFBBD, RFFFD, RFCCD, and RFUD. Results indicate that assigning triangular distribution for all factors decreases the prediction error of most likely recovery factor regardless of DOE type.

SC

However, relative error increases by considering uniform distribution. On the other hand, in case of assigning uniform distribution for all factors, construction of recovery factor distribution using CCD, BBD, FFD and UD respectively generates less prediction

M AN U

error. However, in case of triangular distribution, construction of recovery factor distribution using BBD, CCD, FFD and UD respectively provides better predictions. 5) Conclusions and recommendations 5.1) Conclusions

TE D

In this paper, we studied the uncertainty associated with gas recovery from water drive gas reservoir by applying design of experiment and relative variation factor method. Four designs of Box-Behnken, full factorial, central composite and uniform design have been used for developing recovery response function.

EP

Next, genetic algorithm was applied to optimally combine different DOE and RVF for prediction of final recovery factor distribution. In this regard, part of heterogeneous

AC C

water drive gas reservoir has been selected and simulated with properties described in section 2. It is important to note that the results and conclusions brought as an insight of this study are valid under the specific model and assumed properties and its ranges. In addition to review the most effective parameters, we aimed to introduce the useful methodology for reduction of recovery factor uncertainty. The following conclusions could be drawn from this investigation: 1) The uncertainty associated with ultimate recovery factor depends mostly on average reservoir permeability, permeability anisotropy, tubing head pressure and aquifer size respectively.

Page 13 of 24

ACCEPTED MANUSCRIPT 2) GA optimized and RVF methods have the smallest and largest prediction error for most likely recovery factor, respectively. 3) RVF and UD generate the smallest uncertainty range (P90%-P10%) for recovery factor while the FFD gives the largest range. 4) Among different DOEs and without considering two methods of GA and RVF;

RI PT

applying BBD, CCD, FFD, and UD generate better predictions respectively. 5) Prediction error decreases by assigning triangular for all factors and increases by uniform distributions. In case of uniform distribution, construction of recovery factor distribution using CCD generates less prediction error, and in case of

SC

triangular distribution, BBD has higher prediction performance. 5.2) Recommendations

M AN U

The following items are proposed to improve the use of methodology for uncertainty assessment and optimum integration of different methods.

1) Global and local optimization algorithms other than the GA considered here should be explored for robust optimum combination of methods including particle swarm optimization (PSO), gene expression programming (GEP), bat algorithm (BA), and hybrid models.

TE D

2) DOEs other than BBD, FFD, CCD and UD should be used for generating proxy models in order to understand the effects of proxy performance on the result of proposed methodology.

EP

3) Further blind tests should be performed to improve the framework efficiency. 4) The method developed in this work should be tested on real field applications.

AC C

Nomenclature BBD FFD CCD UD RVF GA KG Kv/Kh Vaq Qg Hp THP VDP P Bg µg

Box-Behnken design Full Factorial design Central Composite design Uniform design Relative variation factor Genetic algorithm Average reservoir permeability Permeability anisotropy Aquifer size Production rate Perforated thickness Tubing head pressure Dykstra-Parson coefficient Pressure Gas formation volume factor Gas viscosity

Page 14 of 24

Gas saturation Gas relative permeability Water relative permeability Capillary pressure Monte Carlo simulation Design of experiment Response surface methodology Coded variable for average reservoir permeability Coded variable for permeability anisotropy Coded variable for aquifer size Coded variable for tubing head pressure Recovery factor Triangular probability density function Coefficient of determination Adjusted coefficient of determination Mean squared error Number of factors or blind simulation tests

SC

Sg Krg Krw Pc MCS DOE RSM X1 X2 X3 X4 RF TPDF R2 R2adj MSE n

[5]

[6] [7] [8]

[9] [10]

[11]

[12] [13] [14] [15] [16] [17]

TE D

[4]

EP

[3]

Anonsen S. I., Eide A. L., and Holden L., Aasen, J. O. 1995. Optimizing Reservoir Performance Under Uncertainty with Application to Well Location, Paper SPE 30710 presented at the Annual Technical Conference & Exhibition held in Dallas, U. S. A., 22-25, October. Box, G. E. P. and Behnken, D. W. 1960. Some New Three Level Designs for the Study of Quantitative Variables, Technometrics, 2 (4): 455-475. Box, G.E.P and Wilson, K.B. 1951. On the experimental attainment of optimum conditions, Journal of the Royal Statistical Society, Series B, 13, 1-45. Box, JF, 1980. R. A. Fisher and the Design of Experiments, 1922-1926. The American Statistician, JSTOR2682986, 34 (1): 1–7. Chierici, G. L., Ciucci, G. M. and Long, G. 1963. Experimental Research on Gas Saturation Behind the Water Front in Gas Reservoirs Subjected to Water Drive, WPC-10134, Proceeding of the 6th World Petroleum Congress, Frankfurt am Main, Germany, 25 June. Chierici G.L. and Ciucci G.M. 1967. Water Drive Gas Reservoirs: Uncertainty in Reserves Evaluation From Past History, J. Pet. Tech. 19 (2): 237-244. SPE-1480-PA. Corre, B. et al. 2000. Integrated Uncertainty Assessment for Project Evaluation and Risk Analysis, paper SPE 65205 presented at the SPE European Petroleum Conference, Paris, 24-25 October. Cronquist C. 1973. Effects of Permeability Variation and Production Rate on Recovery from Partial Water Drive Gas Reservoirs, Paper SPE-4635-MS presented at the SPE Annual Fall Meeting held in Las Vegas, Nevada30 September-3 October. Damsleth, E., Hage, A., and Volden, R. 1992. Maximum Information at Minimum Cost: A North Sea Field Development Study With an Experimental Design, JPT Dejean, J. P. and Blanc, G. 1999. Managing Uncertainties on Production Predictions Using Integrated Statistical Methods, Paper SPE 56696 presented at the Annual Technical Conference and Exhibition held in Houston, Texas, 3-6 October. Dranchuk, P. M., Purvis, R. A. and Robison, D. B. 1974. Computer Calculations of Natural Gas Compressibility Factors Using the Standing and Katz Correlation, Institute of Petroleum Technical Series, No. IP 74-008. Dykstra, H., and Parsons, R. L., 1950. The Prediction of Oil Recovery by Water Flood, In Secondary Recovery of Oil in the United States, 2nd ed., pp. 160–174. API. Ebrahimi A., Khamehchi E., 2015. A robust model for computing pressure drop in vertical multiphase flow, JNGSE, 26: 1306-1316 Efroymson, MA. 1960. Multiple regression analysis. In Ralston, A. and Wilf, HS, editors, Mathematical Methods for Digital Computers. Wiley. Fang, J.H., et. al. 1992. Genetic Algorithm and Its Application to Petrophysics, SPE 26208. Fang, K. T. 1980. The uniform design: application of number-theoretic methods in experimental design. Acta Math. Appl. Sinica 3, 363-372 Fetel E. and CaumonG. 2008. Reservoir flow uncertainty assessment using response surface

AC C

[2]

M AN U

References [1]

RI PT

ACCEPTED MANUSCRIPT

Page 15 of 24

ACCEPTED MANUSCRIPT

[26] [27] [28]

[29] [30]

[31] [32] [33]

[34]

[35] [36] [37]

[38] [39]

[40]

[41]

[42] [43]

RI PT

[24] [25]

SC

[23]

M AN U

[21] [22]

TE D

[20]

EP

[19]

constrained by secondary information, JPSE, 60 (3-4): 170–182. Firoozjaee R. A.and Khamehchi E., 2015. A Novel Approach to Assist History Matching Using Artificial Intelligence, Chemical Engineering Communications, 202 (4): 513-514. Friedmann, F., Chawathé, A., and Larue, D.K. 2003. Assessing Uncertainty in Channelized Reservoirs Using Experimental Designs , SPEREE. Geffen, T.M., Parrish, D.R., Haynes, G.W. and Morse, R.A. 1952. Efficiency of Gas Displacement from Porous Media by Liquid Flooding, J. Pet. Tech. 4 (2): 29-38. SPE-952029-G. Gen, M., Cheng, R., 1997. Genetic algorithms and engineering design, John wiley& Sons, Inc. Goldberg, D. E., 1983. Computer Aided Gas Pipeline Operation Using Genetic Algorithms and Rule Learning, Ph.D. dissertation, University of Michigan, Ann Arbor, Michigan. Guerreiro, J.N.C., et. al., 1997. Identification of Reservoir Heterogeneities Using Tracer Breakthrough Profiles and Genetic Algorithms, SPE 39066, Latin American and Caribbean Petroleum Engineering Conference and Exhibition held in Rio de Janeiro, Brazil, 30 August-3 September. Hajela, P. 1990. Genetic Search-Optimization Problem, AIAA. 28 (7): 1205- 1210. Holland, J.H., 1975. Adaptation in Natural and Artificial Systems. University of Michigan Press, Ann Arbor. Khamehchi E., RahimzadehKivi I., Akbari M., 2014. A novel approach to sand production prediction using artificial intelligence, JPSE, 123: 147-154. Lee, A. L., Gonzalez, M. H. and Eakin, B. E. 1966. The Viscosity of Natural Gases, J. Pet. Tech. 18 (8): 997-1000. SPE-1340-PA. LepineO.J.,BissellR.C.,Aanonsen S.I., Pallister I.C. and Barker J.W. 1999. Uncertainty Analysis in Predictive Reservoir Simulation Using Gradient Information, SPE J. 4(3): 251 – 259. SPE-57594PA. Liu C. and McVay D.A. 2010. Continuous Reservoir-Simulation-Model Updating and Forecasting Improves Uncertainty Quantification, SPE J. 13(4): 626 – 637. SPE-119197-PA. Malachowski, M. A., Yanosik, J. L., Saldana, M. A. 1995. Simulation of Well Productivity Losses Due to Near Well Condensate Accumulation in Filed Scale Simulations, Paper SPE30715 presented at the SPE Annual Technical Conference & Exhibition, Dallas. Maschio C. and Schiozer D.J. 2013. A new procedure to reduce uncertainties in reservoir models using statistical inference and observed data, JPSE. 110: 7–21 October 2013. McVay D.A. and Dossary M.N. 2014. The Value of Assessing Uncertainty, SPE J. 6(2):100-110. Mohaghegh, S., Popa, A. S. and Ameri, S. 1999. Intelligent Systems Can Design Optimum Fracturing Jobs, SPE 57433, Proceedings, SPE Eastern Regional Conference and Exhibition, October 21-22, Charleston, West Virginia. Myers, R.H. and Montgomery, D.C. 1995. Response Surface Methodology: process improvement with steepest ascent, the analysis of response Surfaces, experimental designs for fitting response surfaces, 183-351. New York: John Wiley and Sons, Inc. Naderi M., Rostami B. and Khosravi M. 2014. Optimizing production from water drive gas reservoirs based on desirability concept, JNGSE, 21: 260-269. Naderi M., Rostami B. and Khosravi M. 2015. Effect of heterogeneity on the productivity of vertical, deviated and horizontal wells in water drive gas reservoirs, JNGSE, 23: 481-491. Peake, W.T., Abadah, M. and Skander, L. 2005. Uncertainty Assessment Using Experimental Design: Minagish Oolite Reservoir, Paper SPE 91820 presented at the Reservoir Simulation Symposium, Houston, Texas U.S.A., 31 January – 2 February. Petalas, N., and Aziz, K. 2000. A Mechanistic Model for Multiphase Flow in Pipes, J. Cdn. Pet. Tech. 39 (6): 43-55. Rasouli H., RashidiF., Karimi and Khamehchi E., 2015. A Surrogate Integrated Production Modeling Approach to Long-Term Gas-Lift Allocation Optimization, Chemical Engineering Communications, 202 (5): 647-654. Salomao, M.C. and Grell, A.P. 2001. Uncertainty in Production Profiles on the Basis of Geostatistic Characterization and Flow Simulation, paper SPE 69477 presented at the SPE Latin American and Caribbean Petroleum Engineering Conference, Buenos Aires, Argentina, 25 – 28 March. Sen, M.K. et. al., 1992. Stochastic Reservoir Modeling Using Simulated Annealing and Genetic Algorithm, SPE 24754, SPE 67th Annual Technical Conference and Exhibition of the Society of Petroleum Engineers held in Washington, DC, October 4-7. Simons, L.H., Spain H.H. 1958.Limitations on Pressure Predictions for Water-Drive Reservoirs, J. Pet. Tech. 10 (10): 51-53. SPE-915-G Steagall, D.E. and Schiozer, D.J. 2001. Uncertainty Analysis in Reservoir Production Forecasts

AC C

[18]

Page 16 of 24

ACCEPTED MANUSCRIPT

[50]

RI PT

[49]

SC

[48]

M AN U

[47]

TE D

[46]

EP

[45]

AC C

[44]

During Appraisal And Pilot Production Phases, paper SPE 66399 presented at the SPE Reservoir Simulation Symposium, Houston, 11 – 14 February. Subbey S., Christie M. and Sambridge M. 2004. Risk Analysis Applied to Petroleum Exploration and Production Prediction under uncertainty in reservoir modeling, JPSE. 44 (1-2): 143–153 Valipour M., Montazar A. A., 2012. Optimize of all Effective Infiltration Parameters in Furrow Irrigation Using Visual Basic and Genetic Algorithm Programming, Australian Journal of Basic and Applied Sciences, 6(6): 132-137. Valipour M., Montazar A. A., 2012. Sensitive Analysis of Optimized Infiltration Parameters in SWDC model, Advances in Environmental Biology, 6(9): 2574-2581. Venkataraman R. 2000. Application of the Method of Experimental Design to Quantify Uncertainty in Production Profiles, Paper SPE 59422 presented at the SPE Asia Pacific Conference on Integrated Modelling for Asset Management, Yokohama, Japan, 25–26 April. White, C.D., Willis, B.J., Narayanan, K. and Dutton, S.P. 2001. Identifying and Estimating Significant Geologic Parameters With Experimental Design, SPEJ. Yin, Hongjun, Zhai, Yunfang, 1998. An Optimum Method of Early-Time Well Test Analysis-Genetic Algorithm, SPE 50905, International Oil & Gas Conference and Exhibition in China held in Beijing, China, 2-6 November. Zee Ma Y. 2010. Error types in reservoir characterization and management, JPSE, 72 (3-4): 290– 301.

Page 17 of 24

SC

RI PT

ACCEPTED MANUSCRIPT

M AN U

Fig. 1.The numerical simulated model. The aquifer (in blue) is attached to the bottom of reservoir (in red). This model consists of 40 layers in the vertical direction, 30 layers with total thickness of 300 ft and 10 layers with variable thickness for reservoir and aquifer section respectively. The reservoir and aquifer areal dimension is 3000 ft × 3000 ft and 6000 ft × 6000 ft respectively.

Recovery Factor (%)

45

65

75

85

95

TE D

Average reservoir permeability

55

Tubing head pressure Permeability anisotropy

EP

Aquifer size

Perforated thickness

AC C

Production rate

Fig. 2. The tornado chart of recovery factor with six different factors (upper bound of factors in red and lower bound of factors in blue)

Page 18 of 24

ACCEPTED MANUSCRIPT

Actual Recovery Factor Predicted Recovery Factor

100 90 80 70 60 50 40 30 20 10 0

100 90

100

y = 0.9997x R² = 0.9626

80 70 60

M AN U

y = 0.9989x R² = 0.9788

50

Actual Recovery Factor

B Predicted Recovery Factor

A

0

100

SC

50

y = 0.9988x R² = 0.9764

RI PT

y = 0.9987x R² = 0.9635

0

50 40 30 20 10 0

0

C

100 90 80 70 60 50 40 30 20 10 0

Predicted Recovery Factor

Predicted Recovery Factor

100 90 80 70 60 50 40 30 20 10 0

50 Actual Recovery Factor

0

100

D

50

Actual Recovery Factor

100

EP

TE D

Fig.3. Predicted versus actual experiments for recovery factor based on A) BBD, B) FFD, C) CCD and D) UD

AC C

Average reservoir permeability Permeability anisotropy Aquifer size Tubing head pressure

RFBBD RFFFD RFCCD

GA optimization

RFUD RFRVF

Fig. 4. The schematic diagram of new approach

Page 19 of 24

Final RF

SC

RI PT

ACCEPTED MANUSCRIPT

M AN U

Fig. 5. Plot showing the mean and best fitness values for normalized MSE, andbest, worst and mean scores within 51 generations

1

BBD FFD CCD UD RVF GA Optimized

TE D

0.8 0.7 0.6 0.5 0.4

EP

Cumulative Probability

0.9

0.3 0.2

AC C

0.1

0

0

20

40 60 Recovery Factor

80

Fig. 6. Result of all different uncertainty assessment methods

Page 20 of 24

100

ACCEPTED MANUSCRIPT

8.302392

0.0119

800

4.035185

0.012537

1200

2.627463

0.013346

1600

1.937011

0.014314

2000

1.5347

0.015425

2400

1.276632

0.016655

2800

1.100498

0.017973

3200

0.974824

0.019348

3600

0.881742

0.020752

4000

0.81073

0.022163

4400

0.755129

0.023566

4800

0.710644

0.024947

5200

0.674302

0.026302

5600

0.64413

0.027624

M AN U

SC

400

RI PT

Table 1 - Fluid properties P (psia) Bg (bbl/Mscf) µg (cp)

Table 2- The relative permeability data used in this simulation model Sg Krg Krw Pc (psia) 0

0

1

0

0.875

0.07

0

0.75

0.15

0.3

0.02

0.625

0.23

0.4

0.085

0.5

0.3

0.5

0.2

0.375

0.37

0.6

0.37

0.225

0.45

0.7

0.65

0.125

0.53

0.8

1

0

0.6

0.1

Table 3 - Levels of factors with transfer functions Coded Variable Level (-1) Level (0) Level (+1) X1 1 10 100 X2 0.01 0.1 1 X3 1 10 100

AC C

Factor KG(md) Kv/Kh Vaq

EP

TE D

0.2

THP (psi)

0

X4

500

1000

Page 21 of 24

1500

Transfer Function Log (KG) -1 Log (Kv/Kh) +1 Log (Vaq)-1 345 − 1000 500

ACCEPTED MANUSCRIPT Table 4- Result of 81 flow simulations based on FFD X1

X2

X3

X4

RF

Run

X1

X2

X3

X4

RF

1

1

-1

1

1

72.44

42

-1

-1

1

-1

52.07

1

-1

-1

-1

88.46

43

0

0

1

0

81.97

3

-1

-1

-1

-1

51.46

44

0

1

1

1

58.47

4

0

-1

-1

1

63.09

45

0

-1

0

1

64.77

5

1

-1

0

1

67.46

46

1

1

-1

0

77.96

6

-1

0

-1

-1

53.5

47

0

0

1

-1

87.31

7

0

1

-1

1

63.65

48

1

1

-1

-1

88.79

8

-1

0

1

1

43.46

49

-1

-1

0

0

47.38

9

-1

1

0

-1

52.87

50

0

0

0

-1

84.08

10

0

-1

0

0

75.06

51

1

-1

-1

0

77.15

11

-1

-1

-1

1

39.15

52

0

0

-1

1

63.37

12

-1

-1

1

1

39.46

53

-1

13

1

1

1

1

75.35

54

-1

14

-1

-1

0

1

40.14

55

-1

15

-1

1

0

0

16

1

0

1

-1

17

1

0

0

-1

18

0

-1

-1

-1

19

0

0

-1

0

20

0

-1

1

0

21

-1

0

1

0

22

1

1

1

23

-1

1

24

-1

1

25

1

0

26

1

-1

27

1

1

28

1

29

0

30

-1

31

SC

RI PT

2

M AN U

Run

-1

-1

0

46.13

-1

0

-1

52.6

1

1

-1

59.15

56

0

0

-1

-1

83.03

89.49

57

0

-1

0

-1

82.93

89.33

58

-1

1

1

1

46.82

81.89

59

1

-1

1

-1

90.38

74.02

60

-1

0

-1

0

48.1

75.44

61

-1

0

0

1

42.7

50.99

62

0

1

0

-1

77.72

0

82.16

63

-1

0

-1

1

40.78

-1

0

49.11

64

-1

1

1

0

49.41

-1

-1

54.96

65

1

1

0

0

79.23

1

0

82.43

66

1

0

-1

1

66.25

0

-1

89.2

67

-1

-1

1

0

46.67

0

-1

89.52

68

0

-1

-1

0

73.66

-1

1

0

82.49

69

0

1

0

0

66.92

0

0

0

75.46

70

1

-1

-1

1

65.76

0

1

-1

56.86

71

-1

0

0

0

49.94

0

1

0

1

58.19

72

1

0

-1

-1

88.58

32

1

0

0

0

79.07

73

-1

0

0

-1

55.41

33

0

-1

1

1

65.07

74

1

0

1

1

76.07

34

1

1

1

-1

88.56

75

0

1

-1

-1

82.67

35

0

1

1

-1

79.45

76

1

1

0

1

68.09

36

1

-1

0

0

78.44

77

0

1

1

0

68.06

37

0

-1

1

-1

83.28

78

0

0

0

1

65.13

38

1

0

0

1

67.94

79

1

0

-1

0

77.79

39

0

1

-1

0

74.37

80

-1

1

0

1

42.91

40

1

1

-1

1

66.39

81

0

0

1

1

73.53

41

-1

1

-1

1

41.92

AC C

EP

TE D

49.68

Page 22 of 24

ACCEPTED MANUSCRIPT Table 5 - Response functions statistics

Response RFBBD RFFFD RFCCD RFUD

R2 96.48% 97.69% 97.92% 96.40%

R2adj R2-R2adj 95.30% 1.18% 97.36% 0.33% 96.89% 1.03% 94.80% 1.60%

Factor Average reservoir permeability Permeability anisotropy Aquifer size Tubing head pressure

Lower bound of TPDF 0.662 0.887 0.981 0.863

RI PT

Table 6- shows the calculated lower and upper bound of triangular distribution for all four factors Upper bound of TPDF 1.048 1 1.086 1.114

RFBBD 67.95 88.14 78.03 80.58 65.01 56.62 68.57 73.35 78.97 83.75 69.62 72.71 79.93 74.07 76.93 78.14 57.87 87.00 83.38 69.76 43.87 82.76 75.10 71.99 54.24 78.73 77.81 53.38 68.13 82.71 85.81 74.88 61.54 74.33 72.44 74.84 68.91 78.08 70.12 65.19 77.41 68.04 84.68 60.87 63.98 68.69 71.37 66.46 80.32 71.85

RFFFD 65.38 88.40 77.58 80.77 64.71 54.44 67.96 71.99 78.40 83.87 68.66 71.92 79.20 73.34 77.35 77.00 57.56 86.21 82.47 68.79 43.38 82.53 74.99 71.74 51.79 78.53 75.39 52.46 67.42 82.54 84.94 74.55 62.63 73.31 73.12 72.95 68.10 76.90 70.63 65.55 77.57 67.31 85.57 61.09 65.37 70.22 68.84 64.86 79.93 70.64

RFCCD 68.33 88.14 80.93 81.46 65.83 52.81 70.27 72.13 78.48 84.24 68.97 72.68 80.55 75.00 77.72 77.37 58.35 86.45 82.39 69.70 43.64 82.85 75.01 72.49 51.43 77.47 76.00 52.79 68.69 81.57 83.31 75.39 64.87 74.92 74.46 74.70 69.88 77.53 72.94 67.89 78.31 68.09 88.08 60.08 66.20 71.62 67.21 64.88 79.31 71.83

RFUD 67.10 88.98 74.50 78.80 62.28 60.52 59.64 69.07 73.16 80.14 70.76 72.56 75.98 66.34 73.56 69.77 53.48 85.00 80.07 69.77 46.26 78.97 53.79 67.65 66.01 65.45 75.64 50.45 73.36 86.34 88.29 71.39 55.94 66.77 71.14 75.41 68.94 76.39 68.49 58.98 74.01 69.34 79.27 61.67 64.74 68.52 75.69 62.69 69.53 64.10

M AN U

X4 -0.67 -0.93 -0.31 -0.87 0.83 0.03 0.99 0.25 0.14 -0.97 -0.83 -0.25 -0.02 0.60 0.16 0.28 0.89 -0.73 -0.34 0.05 0.72 -0.22 0.68 0.69 -0.47 0.23 -0.63 0.78 -0.78 -0.71 -0.88 0.35 0.40 0.62 -0.10 -0.57 0.02 -0.15 0.39 0.93 -0.42 -0.19 -0.39 -0.55 -0.06 0.52 -0.51 0.46 0.10 0.85

TE D

X3 1.00 -0.47 0.93 -0.72 -0.15 -0.28 -0.81 -0.61 -0.54 0.24 0.30 -0.39 0.40 0.61 0.09 0.35 0.42 0.19 0.02 0.12 -0.72 0.05 0.90 -0.18 -0.41 0.50 0.57 -0.77 -0.96 -0.86 -0.33 -0.02 -0.99 -0.65 0.65 0.54 -0.91 -0.31 0.77 0.82 -0.59 0.23 0.85 -0.06 0.71 0.73 -0.22 -0.11 0.44 -0.49

EP

X2 -0.50 -0.16 0.27 0.47 0.39 -1.00 0.67 -0.64 -0.36 0.81 0.71 0.31 0.14 -0.15 0.77 -0.41 -0.37 -0.22 -0.31 0.21 -0.76 0.06 -0.88 0.75 -0.45 -0.70 -0.67 -0.87 0.09 -0.82 -0.74 0.55 0.64 -0.01 0.86 0.00 -0.25 -0.08 0.50 -0.08 0.57 0.44 0.34 0.92 0.99 0.94 -0.94 -0.55 -0.60 0.17

AC C

X1 -0.67 0.54 -0.06 -0.01 -0.12 -0.73 0.91 0.09 0.71 0.18 -0.54 -0.21 0.35 0.23 0.80 0.44 -0.51 0.46 0.48 -0.28 -0.94 0.75 0.62 0.65 -0.98 0.92 -0.21 -0.61 -0.59 0.25 0.31 0.36 -0.37 0.58 -0.10 -0.34 -0.28 0.13 -0.10 -0.19 0.05 -0.43 0.85 -0.70 -0.45 0.03 -0.37 -0.26 0.82 0.97

SC

Table 7-Result of 50 new blind simulation tests for GA optimization # 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

Page 23 of 24

RFRVF 64.07 85.27 78.00 81.46 65.30 56.23 58.50 74.36 75.91 79.89 67.88 71.50 78.31 74.94 70.46 79.87 56.51 88.38 83.70 68.14 46.46 80.24 83.59 63.49 52.76 85.92 76.09 53.73 64.51 84.77 90.43 68.45 61.86 67.06 74.63 71.83 66.85 76.05 70.63 63.99 71.08 66.31 85.44 61.06 65.18 66.70 69.63 65.09 85.70 65.57

RFSIM 68.11 86.59 79.10 82.13 65.52 59.10 66.56 72.28 75.04 81.92 72.85 74.59 80.04 75.17 77.57 76.68 61.45 85.64 81.37 72.39 44.84 82.41 72.15 70.53 52.12 78.87 79.75 55.97 69.74 82.41 84.81 74.06 65.67 70.23 66.00 78.99 70.21 77.07 70.85 69.84 79.02 72.60 84.07 62.93 56.63 63.45 73.54 67.46 79.43 68.49

ACCEPTED MANUSCRIPT

RI PT

Table 8- Parameter settings for GA optimization Initial population size 75 Creation function Constraint dependent Initial range [-2, 2] Fitness scaling function Rank Selection function Stochastic uniform Crossover function Scattered Crossover fraction 1 Mutation function Adaptive feasible Stopping criterion Tolerance value of 0.000001

RFCCD RFUD

EP

RFFFD

Table 11 – Summary of the four different methods output Distribution Type Minimum P10% P50% P90% Maximum Triangular 44.98 62.54 74.23 81.98 88.29 Uniform 39.08 55.15 72.76 83.87 90.74 Triangular 47.51 61.49 73.18 81.05 89.17 Uniform 41.93 54.54 72.05 84.02 90.80 Triangular 45.71 62.26 73.90 81.82 90.97 Uniform 40.73 55.31 72.88 84.19 92.58 Triangular 48.46 62.11 71.51 79.48 87.45 Uniform 42.33 55.72 68.76 81.54 89.62

AC C

RFBBD

Table 10 –Summary of the six different methods Minimum P10% P50% P90% Maximum P90%-P10% 42.52 57.05 73.80 81.88 89.30 24.83 44.81 56.59 73.08 81.98 90.13 25.39 43.78 57.56 73.36 82.47 91.67 24.91 39.60 58.94 70.01 79.00 89.58 20.06 46.27 56.18 66.85 76.22 89.25 20.04 40.74 58.99 75.38 83.09 90.95 24.10

TE D

Method BBD FFD CCD UD RVF GA

M AN U

SC

Table 9-Performance of six methods based on MSE and R2 for 50 new blind simulation tests Rank Method MSE R2 BBD 7.11 0.910 2 FFD 9.64 0.880 3 CCD 10.86 0.864 4 UD 35.02 0.630 6 RVF 21.3 0.739 5 GA optimized 6.67 0.913 1

Page 24 of 24

P50% RE (%) 2.20 3.15 2.78 7.22 11.4 0.11

P50% RE (%) 1.63 3.58 3.02 4.52 2.07 3.42 5.23 8.88

ACCEPTED MANUSCRIPT Highlights A new method of uncertainty reduction in water drive gas reservoirs is proposed. Design of experiments and relative variation factor are used in this work. The method combines conventional uncertainty assessment and genetic algorithm.

AC C

EP

TE D

M AN U

SC

RI PT

The proposed model has minimum relative error among different methods.

1