Estimation of Breakthrough Time for Water Coning in Fractured Systems: Experimental Study and Connectionist Modeling Sohrab Zendehboudi, Ali Elkamel, and Ioannis Chatzis Dept. of Chemical Engineering, University of Waterloo, Waterloo, ON N2L 3G1, Canada Mohammad Ali Ahmadi Faculty of Petroleum Engineering, Petroleum University of Technology, Ahwaz, Khuzestan, Iran Alireza Bahadori School of Environment, Science and Engineering, Southern Cross University, Lismore, NSW,Australia Ali Lohi Dept. of Chemical Engineering, Ryerson University, Toronto, ON M5B 2K3, Canada DOI 10.1002/aic.14365 Published online February 2, 2014 in Wiley Online Library (wileyonlinelibrary.com) Water coning in petroleum reservoirs leads to lower well productivity and higher operational costs. Adequate knowledge of coning phenomena and breakthrough time is essential to overcome this issue. A series of experiments using fractured porous media models were conducted to investigate the effects of production process and pore structure characteristics on water coning. In addition, a hybrid artificial neural network (ANN) with particle swarm optimization (PSO) algorithm was applied to predict breakthrough time of water coning as a function of production rate and physical model properties. Data from the literature combined with experimental data generated in this study were used to develop and verify the ANN-PSO model. A good correlation was found between the predicted and real data sets having an absolute maximum error percentage less than 9%. The developed ANN-PSO model is able to estimate breakthrough time and critical production rate with higher accuracy compared to the conventional or back propagation (BP) ANN (ANN-BP) and common correlations. The presence of vertical fractures was found to accelerate considerably the water coning phenomena during oil production. Results of this study using combined data suggest the potential application of ANN-PSO in predicting the water breakthrough time and critical production rate that are critical in designing and evaluating production strategies for naturally C 2014 American Institute of Chemical Engineers AIChE J, 60: 1905–1919, 2014 fractured reservoirs. V Keywords: water coning, fractured media, breakthrough time, production rate, artificial neural network-particle swarm optimization Introduction A high recovery factor is achievable in oil reservoirs with active aquifer support.1–5 However, when high oil production rates are used, this causes water production through upward flow merged with oil known as “water coning” where the deformation of water–oil contact happens.1–5 Water coning is a common problem which reduces oil production and increases water production treatment requirements, and operational expenses. Water coning is caused due to the pressure gradients close to the production well and the imbalance between gravity and viscous forces. Viscous forces push the oil phase into the production well with a pressure drawdown.1–5 The dynamic pressure distribution near a production well causes *Mohammad Ali Ahmadi is the second author of the article in term of contribution. Correspondence concerning this article should be addressed to S. Zendehboudi at szendehb@uwaterloo.ca. C 2014 American Institute of Chemical Engineers V AIChE Journal the lift up of the macroscopic water–oil interface just below the bottom part of a vertical well. In contrast, gravity forces cause the water to remain below the oil saturated region of the reservoir. When the dynamic forces at the wellbore surpass the gravitational forces, the water cone develops toward the perforated part of the well (Figure 1) and finally breaks into the wellbore in the form of an oil/water mixture.1–5 Water coning has been studied by many researchers worldwide.1–10 A number of important aspects such as the critical pumping rate, breakthrough time (BT), and water oil ratio after breakthrough have been well addressed in the literature, often based on laboratory tests and/or simulations. A critical rate corresponding to the maximum rate of oil flow which is free of water is also reported.1–10 If the oil production rate becomes higher than the critical rate, then water breakthrough happens, normally leading to lower oil recovery and lower oil production rate.5–10 The models used to predict the critical rate generally fall into two categories. In the first category, the critical rate is May 2014 Vol. 60, No. 5 1905 Figure 1. Schematic representation of water coning. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] calculated through analytical techniques that use the equilibrium between viscous and gravitational forces by developing an oil potential function. However, the methodologies to determine oil potential differ for various models.4,11–16 For instance, the Laplace equation was used by Muskat and Wyckoff11 for single-phase flow; whereas, Chaney et al.12 and Chierici et al.13 used potentiometric models. In addition, Wheatly’s method14 used an oil potential equation and considered the effect of cone shape. Meyer and Garder4 suggested an inexact function assuming a radial flow regime where they used a simple model to estimate the critical rate at which the stable cone contacts the well bottom.4 An analytical solution was developed by Abbas and Bass15 to compute the critical rate at steady state and pseudosteady state conditions. In their study, a 2-D radial flow model was used to solve the corresponding governing equations through average pressure concept.15 The influence of limited wellbore penetration on the oil productivity was investigated by Guo and Lee16 while determining critical rate. Their method was based on a radial/spherical/combined 3-D flow system where low conductivity is assumed.16 Empirical models are the second group of predictive tools for the critical rate. Schols17 obtained a model through conducting a series of flow experiments. Also, Hoyland et al.18 derived a correlation to estimate critical rate in anisotropic reservoirs. Their model was obtained based on several simulation runs through a numerical simulation technique.18 Oil fields are typically operated at flow rates greater than the critical flow rates that are determined using common correlations; otherwise, the field would be uneconomic. If a well produces above its critical rate, the water cone will breakthrough after a certain time period known as “breakthrough time.”1,2,19–22 There are some correlations in the literature such as Yang and Wattenbarger23 and Rhecam et al.24 to determine this specific time.23–28 Two of the most common methods used are the Sobocinski-Cornelius26 and BournazelJeanson27 techniques. Both approaches were based on the laboratory data and modeling results. Also, they correlated the BT with two parameters, namely a dimensionless time and a dimensionless cone height, respectively.26,27 Correlations developed for homogeneous reservoirs consider the effect of vertical permeability in the form of vertical permeability to horizontal permeability ratio (Kv/Kh). Presence of shale seals and compaction phenomenon lower the vertical per1906 DOI 10.1002/aic meability.1,2,28–30 However, there are cases such as naturally fractured reservoirs (NFRs) where the Kv/Kh is generally higher than one.1,2,28–30 High values of vertical permeability due to presence of fractures accelerate the water coning process that leads to reduction in the critical rate and also BT.1,2,28–30 Empirical models normally include adjustable parameters which represent a simple 1-D system to predict the critical rate and production performance of the displacement.1–5,28–30 The dilemma here is that such correlations have been normally obtained for homogeneous porous media and are not based on comprehensive physics of fluid flow. By contrast, analytical and numerical predictive techniques contain sound physical concepts of multiphase flow in porous media. However, both methods appear to be complex tasks that require fundamental knowledge of the mechanisms contributing in the water coning phenomena. In addition, knowing the mathematical principles of numerical and analytical solutions is necessary to develop equations or compute the target variables. For instance, the number of grid blocks used in the simulation runs should be fairly high to attain accurate predictions. This means that the computational effort will increase and this may lead to divergence issues. Besides developing predictive tools, core flooding tests are common for small-scale laboratory modes to investigate microscopic oil recovery performance and determine optimum production rate.1–3 Experimental data are then used to evaluate the effectiveness of the enhanced oil recovery (EOR) technique and predict the reservoir behavior during the recovery process before any production strategies are field tested. Noteworthy is the fact that core flood trials are usually expensive and time consuming. Some methods have been introduced to diminish the effect of undesired water in production wells such as (a) maintaining recovery rate lower than the critical value; (b) perforating farther a well from the original water–oil interface; and/or (c) injecting cross-linking polymers or gels to make a water blocking zone around the wellbore.1–5,28–30 However, none of these techniques can effectively solve the BT impasse. The artificial neural network (ANN) is considered a robust smart tool which has been broadly implemented across several different disciplines to cope with the uncertainties associated with the nature of materials or processes.31–35 The ANN is capable of modeling ill-defined, uncertain, and nonlinear phenomena. A notable benefit of using this technique is that it does not rely on theoretical acquaintance and/or human experience during the training and testing stages. All unknown relations are predicted with neural network (NN) rather than using conventional relationships. A number of evolutionary algorithms are available such as genetic algorithm (GA), particle swarm optimization (PSO), imperialist competitive algorithm (ICA), and pruning algorithm (PA) that can be joined with ANN to find the proper network structures.35–39 In this article, a feed-forward ANN optimized by PSO was developed to predict the BT and critical rate of water coning in NFRs. An experimental investigation was also performed to understand the physics of water coning in NFRs and the effects of fractured porous media characteristics on water coning. In addition, the experimental phase and the literature survey about the BT help to determine important factors affecting BT and critical rate. The PSO algorithm is used to determine the initial parameter weighting involved in the NN. The ANN-PSO model developed in this article is constructed based on the real data (e.g., 90 laboratory data Published on behalf of the AIChE May 2014 Vol. 60, No. 5 AIChE Journal points, 40 pilot plant-scale data points, and 30 data points from the actual reservoir cases) that are acquired from the current study and from available results in the literature.15– 19,24,25,40–51 The results obtained in this article show the superiority of hybrid ANN-PSO model compared to conventional ANN model and common equations for estimating BT and critical rate. In addition, this study can help in choosing an appropriate production rate for a given petroleum reservoir. Theory Water coning High water cut is prevalent issue during oil recovery from conventional oil reservoirs which is normally caused by either the normal rise of oil/water interface or water fingering and water coning. High water cut leads to an increase in operational expense and decrease in production performance (less oil is produced).1–5 The majority of predictive techniques for water coning provide an expression and/or procedure for obtaining “critical rate” at which a stable cone exists from the oil/water interface to the nearby perforations. However, the theories behind the critical rate are usually unable to forecast when breakthrough will take place. A considerable amount of the world’s oil production comes from the naturally fractured reserves located mostly in the Middle East and USA.45,50,52–54 This type of reservoirs is usually characterized by fractures with high permeability and low porosity.45,50,52–54 High oil production rates make the bottom water to push the oil phase to form a cone causing an early breakthrough into the production well.45–50 Thus, an improved management over the production process appears to be vital to avoid water coning and also enhance the oil production. Although water coning has been studied comprehensively for homogeneous reservoirs, just a few research works address key aspects of water coning in fractured reservoirs.41–45 Hence, it seems beneficial to perform experimental and theoretical research to obtain the critical production rate and BT for fractured porous systems. According to Muskat’s methodology,2 the water BT due to water coning in a porous medium is determined as follows BT 5 0:33/e d 2 lo lnð2d=rw Þ KðPe 2Pw Þ (1) where BT refers to the breakthrough time, /e represents the effective porosity, d is the distance between the drainage boundary and the wellbore, rw is defined as the wellbore radius, K is the symbol for permeability, and Pe and Pw are the reservoir pressure and the wellbore pressure, respectively. Equation 1 was originally obtained for homogeneous porous media; however, it can be used for heterogeneous or fractured reservoirs if the permeability and porosity are modified for these types of porous media.2 For fractured porous systems, Saad et al.42 obtained Eq. 2 for prediction of BT in terms of production rate (q), fracture aperture (b), and the distance between the gas/liquid interface, and the production well (de) through modifying Eq. 1 BT 5 0:33pbde2 q (2) There are two main drawbacks of Eq. 2: (1) only one vertical fracture in the experimental study has been used and (2) the equation does not offer reliable results for oils with moderate to high viscosity. According to Saad’s results,42 AIChE Journal May 2014 Vol. 60, No. 5 fairly big differences were found between the experimental data and predicted values. In addition, a numerical modeling simulation was performed by Zamonsky et al.51 to investigate the behavior of water coning with respect to formation characteristics. Based on the results obtained from several simulation runs, an equation to compute the water BT was developed.51 It is important to note that Eqs. 1 and 2 have not been used to develop the hybrid ANN in this study. The main reasons to address these two equations are as follows: 1. Available correlations (e.g., empirical, analytical, and statistical) are not proper to predict BT in fractured porous media as they have been obtained for particular cases (homogeneous or/and fractured media with vertical fractures). Needless to mention that in the case of some modifications, they can be used to determine a rough estimate of BT in oil fractured reservoirs, but the magnitude of error would be at fairly high level. 2. Effects of input or independent parameters have not been considered well while deriving the equations. The two points listed above convey the message that developing a systematic and accurate method for prediction of BT and critical production rate is inevitable as the new technique (e.g., smart tools) presented in our articles captures the nonlinearity nature of these parameters with respect to process (or phenomenon) variables. Artificial neural network The first wave of interest in ANNs appeared after McCulloch and Pitts introduced simplified neurons in 1943.33–37,55 ANNs represent computing systems whose main idea has arisen from the analogy of biological NNs. They stand for very simplified mathematical systems of biological NNs. ANNs have the capability to learn and find the relationships between independent and dependent variable(s) from sample patterns to generate consequential solutions to problems even when the input data include errors or are shortened. Also, ANNs can process information quickly and adjust solutions in case if circumstances or system conditions change.33–37,55 In general, ANNs are effective models to simulate biological NNs by establishing numerous simple neurons in an efficient communication manner. The neuron receives inputs from a single or various sources and generates outputs through simple computation using nonlinear functions.33–37,55 The network structure contains a series of neurons linked together which are typically arranged in a number of layers. Each node in a layer obtains and processes input information from the last layer and transfers its output to neurons in the next layer by links. Each link has a certain weight, which indicates connection power. The summation of weights entered to a node is translated to an output using a transfer function.33–37,55 The majority of ANNs include three layers or more; an input layer consisting of input data, an output layer which produces a suitable answer to the specified input, and one or more hidden or middle layers that operate as a group of parameters detectors. A general schematic of an ANN system is demonstrated in Figure 2. Obtaining the proper network structure is one of the most vital aspects, but one of the most complicated tasks in constructing an ANN model, as well. There are various NN categories, depending on the specific characteristics. One of the common classifications is based on the training algorithm such that supervised and unsupervised networks are two major categories.33–37,55 Unsupervised Published on behalf of the AIChE DOI 10.1002/aic 1907 According to the fundamentals of the PSO algorithm, particle swarm contains “n” particles, and the location of each particle represents the potential solution in the multidimensional space. The particle updates its condition to: (a) maintain its inertia; (b) alter its condition based on its optimum status (individual); and (c) adjust the condition on the basis of the swarm optimum position (social).35,65–68 Following the above procedure as shown in Figure 3, the velocity of each particle is calculated as follows65–68 vk11 5w vki 1c1 k1 ðpBest 2ski Þ1c2 k2 ðgBest 2ski Þ i i Figure 2. A schematic of ANN structure containing three layers. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] networks are also recognized as self-organizing ANN plans. In the unsupervised ANN, no feedback is given during the training process. However, both input and output data are provided for the network in the supervised training algorithm and the network is trained on the basis of a feedback.33–37,55 ANN models need to be trained and then tested. Links weights are regulated to attain the desired outputs in the training stage. In the testing process, the precision of forecasted output data is assessed using results not used during the training process. The training and testing stages are iterated through adjusting weighting parameters until the magnitude of error is minimized.33–37,55 The commonly used error function is the mean squared error (MSE) as follows56–58 MSE 5 G X m h i2 1X ypj ðkÞ2yaj ðkÞ 2 k51 j51 (4) where vki is the velocity of particle i at iteration k, w presents the weighting function, cj is the weighting factor, kj is the random number in the range of 0–1, ski is the current position of particle i at iteration k, and pBest and gBest are the individi ual best position of particle i and global best position for the surroundings or group, respectively. The weighting function normally used in the PSO algorithm is65–68 w5wmax 2 ðwmax 2wmin Þiter iter max (5) Here, “wmax” presents the initial weight, “wmin” is the last weight, “itermax” is the maximum number of iteration, and “iter” is the current iteration number. The current position of the particle is updated as follows65–68 sk11 5ski 1vk1 i i (6) Further information regarding the PSO algorithm can be found elsewhere.35,62–68 (3) where m is the number of output links, G is the number of training data points, ypj represents the predicted output, and yaj is defined as the actual output. The reader is encouraged to see these references for more details on the ANN technique.33–37,55,58–64 Particle swarm optimization PSO is a global optimization technique introduced by Kennedy and Eberhart.35,65–68 PSO has been built up from swarm intelligence, on the basis of studies dealt with the behavior of bird and fish flock movement. Due to its simplicity and high performance, the algorithm can be used extensively in various fields such as function optimization, neutral network training, signal procession, process engineering, reaction engineering, and reservoir engineering.35,65–68 Each individual in the PSO algorithm is labeled as a “particle,” and is subject to a movement in a multidimensional space corresponding to the belief space. Individuals contain memory. Therefore, they keep a part of their earlier conditions. Particles have no limitations to maintain the same characteristics in the hyperspace; however, their individual properties remain unchanged in this case. Movement of each particle depends on the initial random velocity and two random weights based on individual and social effects. The respectively individual influence causes the particle to revisit its own best last position and the social effect leads the particle toward the surrounding’s best previous position.35,65–68 1908 DOI 10.1002/aic Figure 3. Schematic of the PSO algorithm. Published on behalf of the AIChE [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] May 2014 Vol. 60, No. 5 AIChE Journal constructed by milling grooves on Plexiglas strips of known lengths with a thickness which is equal to that of the model. Each fracture is covered by a wire mesh to prevent glass beads from entering into the fracture space. The fractures are used in various orientations and patterns in a fracture network embedded in a matrix of surrounding glass beads. The matrix porosity and permeability are in the ranges of 33– 40% and 20–100 Darcy, respectively. The fracture aperture is varied between 0.5 and 2 mm. In addition, the physical models used have dimensions ranging from 40 to 80 cm (height) 3 80–180 cm (width) 3 1–10 cm (thickness). It should be noted that the presence of fractures is indicated with black strips in Figure 4. Test liquids Figure 4. Schematic diagram of the experimental setup. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] Experimental Work Experimental procedure Several experimental trials were conducted to cover the lower end of the reservoir spectrum and also explore more aspects (e.g., physics and effect of fracture) of water coning in fractured porous media (Figure 4). The data produced in these experiments in conjunction with the results from 100 additional water coning tests (e.g., 30 laboratory data, 40 pilot plant data, and 30 actual reservoir data) reported in literature were used in this study.15–19,24,25,40–51 An experimental study using a particular set-up was performed by the authors of this article to further investigate important aspects (e.g., physical concept and effect of fracture) of critical rate and BT during production, prior to implementation of the hybrid ANN model. When the ANN model was formulated, all the data borrowed from open literature plus the experimental dataset of the current research were used in the modeling approach as input data. It is important to note that the data taken from other works include laboratory scale and large world scale cases. Also, the physical models in these studies are different in terms of geometry, properties, and dimensions. These differences exist in the form of various values for permeabilities, critical rate, and distance between drainage boundary and the production well. This novelty of the ANN model is in its potential to train with combining wide ranges of data from the real world and experiments. The data used in model training cover whole spectrum of the main parameters. Experimental setup The experimental setup is made of the following components (see Figure 4): (a) a rectangular porous system; (b) glass beads to act as the porous medium; (c) certain networks of fractures; (d) high-definition video-recording and high-resolution digital-imaging facilities; (e) a peristaltic pump with variable injection and production rates; (f) a digital balance to weigh the withdrawn liquid; and (g) a vacuum apparatus to remove dissolved gas from the test liquids. The setup is schematically shown in Figure 4. The porous media system is composed of transparent to flow porous media created by placing glass beads in a rectangular in shape box made of glass plates. The physical model has three holes in its bottom for water injection and one hole on the top as the production well. The fractures are AIChE Journal Immiscible oil and water are used to study the behavior of water coning in fractured media produced through vertical wells. Two types of oil with different properties are used in the experiments as shown in Table 1. May 2014 Vol. 60, No. 5 Four ports (each two in one side) are used to provide oil from a large oil storage tank for the physical model to simulate an infinite fractured reservoir at steady state condition. This design also removes the boundary limitation introduced by the physical dimensions anyways. Therefore, all four valves connected to the oil tank are kept open during the production process. Moreover, three ports at the bottom of the model are used to sustain continuous water supply from a water container at steady state condition. It should be noted here that experimental runs were performed at the temperature of 28 C. The experimental procedure is as follows: 1. First, the porous model equipped with a particular fracture network is filled with glass beads, and then the porous system is flooded with water from the bottom side. Then, the oil is injected from both right and left sides of the physical model. The initial thickness of oil layer and the water/oil interface position can be simply varied within a wide range through changing the rates of oil and water injected to the fractured model. 2. The outlet valve of the model is opened to produce oil at a constant recovery rate. The peristaltic pump is used to produce the oil from the fractured models. A broad range of production rates is applied to each certain fractured medium during the experiments. 3. The rise in the horizontal oil/water interface is observed as the oil is produced from the production well. To noticeably see the water–oil contact (WOC) movements during the experimental trials, the water is colored with food dye. 4. The point at which gas/liquid interface commences to form a cone is determined. The distance between the producing well and the cone apex corresponds to the critical radius (rc). 5. For each trial, the water BT is recorded. Observations of water cone development and measurements of oil Table 1. Properties of Test Fluids Used in the Experiments Type of Fluid Viscosity (cp) Density (gr/cc) 1.95 0.46 1.01 0.80 0.72 0.99 Kerosene Gasoline Water Published on behalf of the AIChE DOI 10.1002/aic 1909 production are conducted every 15 s until the water cut is approximately 95%. 6. The above steps are repeated for various rates, fracture aperture sizes, and fluid types. Although, the experimental study reported in this manuscript was carried out by the authors in the ambient (or laboratory) conditions, the ranges of pressure and temperature that the majority of the data (taken from other studies) have been obtained at are 1.0–800.0 atm and 20–160 C, respectively. These intervals are typical for most of conventional light oil reservoirs in the world. seven (7) hidden neurons, and one (1) output neuron. ANN model was trained using back propagation (BP) network through the Levenberg-Marquardt method to predict the amount of BT (one output) from three input parameters namely, the height (H), the fracture to matrix permeability ratio (Kf/Km), and the flow rate (q). The transfer functions are sigmoid and linear in the hidden and output layers, respectively. Further details about the procedure of ANN modeling along with the MATLAB code are given in Appendix A. Results and Discussion ANN-PSO Model Prior to implementing the ANN-PSO modeling, the production rate (q), model height or oil zone thickness (H), and the ratio of fracture to matrix permeability (Kf/Km) are selected as the input variables for the ANN network, based on the experimental investigation, literature review, and determination of key physical parameters of the water coning process in fractured media. Normally, the ANN consists of three important steps, namely training, selection, and testing. The training process is used to train the network, the second stage to stop overtraining process, and the testing phase is used to ensure that the system results after network selection are generalized appropriately. The available data (both real world and experimental) are divided into two series including, a training data set and a testing data set. During the training process, the aim is to determine the optimum values for the weights of the links in the network layers. If the number of weights is greater than the number of available data, the magnitude of error in fitting the nontrained data primarily lowers; however, it increases as the network turns out to be overtrained. On the contrary, when the number of weights is lower than the number of data used, the possibility of over fitting issue does not exist. The size of data in the training phase should be optimized to remove pointless data for the cases where lowering the data points does not considerably affect the precision of predictions. Optimization of the number of hidden neurons is also carried out to drop off the time required for the ANN to forecast the objective variable. The PSO is developed to optimize the NN and the MSE is used to measure the fit of the optimization. The objective in the proposed hybrid ANN-PSO is to minimize the error in the specific algorithm proposed in this study. Before starting the ANN-PSO implementation, the parameters are normalized to be able to account for differences in dimensions and type (e.g., production rate and height). The following expression is used to normalize the data values so they fall between 21 and 11 Xi 5 2ðxi 2xmin Þ 21 ðxmax 2xmin Þ (7) where Xi represents the normalized value of the parameter “xi”, and xmax and xmin are the maximum and minimum magnitudes in all observations for that study. Statistical parameters such as the coefficient of determination (R2), MSE, minimum absolute percentage error (MIPE), and maximum absolute percentage error (MAPE) are used in this study to evaluate the performance of the ANN models. The best ANN architecture for the case under study (e.g., BT) was determined to be 3-7-1: three (3) input neurons, 1910 DOI 10.1002/aic In this article, an ANN optimized with PSO was used to develop a predictive tool to find the value of BT during production process in fractured media if water coning occurs. Data used in this study were obtained from various reports and/or articles available in the literature (100 data points) and from the experimental work (60 data points) conducted to cover whole spectrum in implementing a comprehensive ANN-PSO modeling.15–19,24,25,40–51 These data have been obtained from the trials that include the porous media with different fracture patterns. A part of the water coning results is presented in Appendix B. Gathering the experimental data from dissimilar sources, this study envelops wide ranges of fluids properties, and porous media with different characteristics. Experimental work The experimental trials were performed with different fracture and matrix permeabilities. As a sample, the production data of the porous systems (one homogeneous and two fractured systems) with matrix permeability of 30 Darcy and matrix porosity of 37% are provided in this article. The fractured porous media have Kf/Km of 10 and 40, respectively. It is observed in all experiments that oil/water interface goes upward gradually, in the form of a uniform horizontal plane, and then at a specific distance or radius from the production point, this contact begins to moves toward the production well that leads to formation of a cone. The coning process exacerbates and rapidly arrives at the producing well, resulting in production of oil and water jointly. It is also concluded that the critical radius is not dependent on the location of the primary oil/water interface. The BT vs. production rate for the homogeneous and fractured porous systems are illustrated in Figure 5. The comparison of these porous systems with respect to BT shows that the presence of fractures causes earlier breakthrough as it increases the effective vertical permeability, leading to formation of an easy path flow for the water and oil phases. It is also found that higher fracture permeability results in lower value of BT, meaning that lower recovery factor at breakthrough is attainable. Figure 6 also describes the production behavior of homogeneous and fractured porous models with respect to BT at various oil zone thicknesses during the oil recovery. According to this figure, increasing oil zone thickness increases the BT; however, this increase is lower in the fractured porous medium compared to the unfractured one. The ANN-PSO technique Data selection is an essential stage in the training process to enhance the performance of an ANN model. It is important to choose good training data that cover all the possible conditions of water coning phenomenon. System complexity is a central feature in determining the required quantity of Published on behalf of the AIChE May 2014 Vol. 60, No. 5 AIChE Journal Figure 5. BT in terms of recovery rate for porous systems with different fracture permeability (Kf; Km 5 30 D, Matrix porosity 5 37%, H 5 100 cm). [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] data for ANN modeling.35,62–68 ANN performance strongly depends on the quality and precision of the data used. Therefore, the selected data for ANN training should have the maximum accuracy and the minimum improbabilities. Choosing a representative data series including all the complexities and nonlinearities of the system behavior under study will help assure prediction accuracy.35,65–68 Here, it is presumed that the sample data, combination of the real world and experimental data, selected for the ANN model statistically represent the system. Parameters affecting dataset size in the ANN consist of the problem type, complexity, possible relationships between factors, user experience, precision, network type, and time required to solve the problem. In ANN applications, data selection is normally performed randomly. If the size of training data is very small, the results of network may not be reliable. Therefore, a proper number of data can guarantee the validity of the both training and testing stages. Data not previously used should be used in testing phase. In this Figure 6. BT for fractured and unfractured models vs. oil thickness (Matrix permeability: 30 Darcy, Matrix porosity: 37%). [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] AIChE Journal May 2014 Vol. 60, No. 5 study, we used 120 data points (40 experiments and remaining from the literature) for training and 40 data for testing. The ranges of oil zone thickness, oil viscosity, production rate, ratio of fracture permeability to matrix permeability, and density difference are 0.3–65 m, 0.4–35 cp, 0.5–8000 STB/day, 2–150, and 0.01, 0.55 g/cm3, respectively. In this study, both feedforward ANN system and combined ANN system were examined to predict BT. According to the results obtained, although the latter model attained better forecasts for the training stage, the improvement was really low. On the other hand, finding the optimum ANN structure and conducting training stage is more difficult and also lengthy in terms of time, compared to the similar procedure in feedforward ANN model using BP. In addition, the joined feedforward-feedback (or recurrent) ANN exhibited a little lower predictive performance (about 1.5%) than the feedforward system during the testing and validation stages. The reason might be that the combined system should convey the information from back to forward and vice versa, leading to be confused and/or unstable when the target function (e.g., BT) has multiple local minima. Given the above justification, the feedforward ANN model combined with an evolutionary algorithm (PSO) was chosen for the purpose of BT prediction. As explained previously, the PSO was considered in this study to optimize the developed NN system. The weights of the network training phase were chosen as variables of an optimization problem. The MSE was implemented as a cost function in the PSO algorithm. BP is a gradient descent algorithm on the error space which most likely gets trapped into a local minimum making it entirely dependent on initial settings (weights). This shortcoming can be removed if an evolutionary algorithm such as PSO that has global searching ability is used. The main goal to apply the PSO algorithm was minimizing the cost function. Every weight in the network was initially set in the range of (21, 1). The ANNPSO model was trained by 50 generations, followed by a BP training technique. The learning coefficient and momentum correction factor were set at 0.7 and 0.001, respectively, for the BP training algorithm. The ANN-PSO model was tested with 3, 5, 6, 7, and 8 neurons in the hidden layer. The results achieved while using various numbers of hidden neurons are tabulated in Table 2. The accuracy degree of the predictions was improved when the number of hidden neurons was increased from 3 to 7. It was found that the network does not show enough degrees of freedom to learn the process perfectly, when 3, 5, and 6 hidden neurons were assigned to the hidden layer. However, a decline in the performance of ANN-PSO was observed at hidden neurons of 8. Thus, the optimum value for the hidden neurons is determined 7, as it takes a long time for the ANN-PSO model at 8 hidden neurons to learn and probably the data are over fitted. The same procedure was followed for other parameters of the ANN algorithm. Using the trial and error technique, the optimum network of the ANN-PSO model has 20 swarms. The c1 and c2 were obtained as 2 each. Also, 0.45 and 0.55 were considered for k1 and k2 parameters, respectively. The maximum value for the inertia weight factor was taken as 0.9; whereas, the final (minimum) value of the velocity inertia was set at 0.4. To evaluate the capability and dominances of the hybrid ANN-PSO algorithm, a back propagation neural network Published on behalf of the AIChE DOI 10.1002/aic 1911 Table 2. Performance of the ANN-PSO in terms of the Number of Hidden Neurons Number of Hidden Neurons 3 5 6 7 8 Training R 2 0.8136 0.8792 0.9443 0.9785 0.9351 Testing 2 MSE MIPE (%) MAPE (%) R MSE MIPE (%) MAPE (%) 0.0941 0.0672 0.0415 0.0287 0.0463 7.4568 6.6679 5.7437 3.5042 5.7894 14.5176 12.8634 9.5273 7.4052 9.9867 0.7425 0.8693 0.9531 0.9207 0.8796 0.1523 0.0742 0.0495 0.0266 0.0499 10.7892 9.1963 8.4321 6.7592 8.9356 15.5186 13.4951 11.1764 9.8552 12.1113 (ANN-BP) was built with the same data used in the ANNPSO model. The performances of the training and testing phases were examined by plotting real data and predicted BT vs. number of samples as shown in Figures 7 and 8, respectively, for the ANN-BP and ANN-PSO models. As clearly seen, the model predictions obtained from the training process are in a good agreement with the real water coning data for both networks (Figures 7 and 8), though ANNBP does not exhibit a reasonable performance in estimating BT over the testing stage. This implies that training the ANN system by PSO can results in better outputs, compared to the conventional ANN (Figures 7 and 8). A comparison of predicted BT using ANN-BP and the observed BT is presented in Figure 9. The predictions that match real values should be placed on the diagonal line (Y 5 X). Although, the training stage shows an acceptable match between predicted and real results (R2 5 0.9438), the testing performance of ANN-BP to estimate of the objective variable in this study is not satisfactory (R2 5 0.9062; Figure 9). Figure 10 in the form of a scatter diagram shows ANN-PSO outputs vs. the real data. A high correlation between observed and predicted magnitudes of BT is observed in Figure 10 as the correlation coefficient (R2) values are 0.9913 and 0.9888 for the ANN-PSO training and testing phases, respectively. As seen in Figure 10, all the data lie in the line of Y 5 X which proves the precision of the ANN-PSO model. It was found that the ANN-PSO has the potential of not being trapped in the local optima. The main reason is that the ANN-PSO model has global searching ability of the PSO as well as local searching ability of the BP. The high value of R2 gives initial impression that ANNPSO technique is constructive and applicable to predict behavior of water coning in fractured porous media. An ideal prediction has an R2 of one (1). However, a correlation coefficient of one does not essentially show a perfect prediction, though R2 is a good indicator of an ANN performance in practice.56,57 It also offers an easy means to compare the performances of various ANN models. However, other statistical parameters such as MSE, MIPE, and MAPE are required to examine the effectiveness of an ANN system.56,57 To have Figure 7. Measured vs. predicted BT for the ANN-BP model: (a) training and (b) testing. Figure 8. Measured vs. predicted BT for the ANN-PSO model: (a) training and (b) testing. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] 1912 DOI 10.1002/aic Published on behalf of the AIChE May 2014 Vol. 60, No. 5 AIChE Journal Figure 9. Performance of the ANN-BP model in estimating BT in terms of R2: (a) training and (b) testing. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] apposite comparison, Table 3 lists R2, MSE, MIPE, and MAPE to validate the developed NNs. The parameters R2 5 0.9888, MSE 5 0.02397, MIPE 5 3.4%, and MAPE 5 8.9% for the ANN-PSO are obtained in contrast to R2 5 0.9062, MSE 5 0.03657, MIPE 5 6.3%, and MAPE 511.2% for the ANN-BP, showing superiority of the ANN-PSO model in prediction of BT. Thus, it is concluded that using an evolutionary optimization strategy like PSO demonstrates an excellent performance in terms of convergence rate and achieving global optima. The plots of MSE against number of epochs are shown in Figures 11 and 12 to demonstrate the performances of the ANN-PSO and ANN-BP models, respectively. The curves within the figures show the relationship between the testing, training, and validation phases initiated for estimation of BT during oil production from fractured systems. As designated with brown circles on Figures 11 and 12, the best efficiency for the validation stage was obtained at MSE 5 0.03657 which was characterized with epoch of 11 for the ANN-BP model, while the ANN-PSO model experienced the paramount performance (MSE  0.02397) for the validation process when the number of epochs is set at 8. According to Figures 11 and 12, the convergence rate of the proposed algorithm (ANNPSO) is noticeably greater than that of the ANN-BP model. A sensitivity analysis was carried out for the ANN-PSO model through the analysis of variance (ANOVA) method56,57 by testing the influence of the independent variaAIChE Journal May 2014 Vol. 60, No. 5 Figure 10. Performance of the ANN-PSO model in estimating BT in terms of R2: (a) training and (b) testing. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] bles [e.g., flow rate (q) and height (H)] on the dependent variable (e.g., BT). The outcome of this analysis is summarized in Figure 13. The greater dependence between any input parameter and the output parameter implies superior importance of that variable on the extent of the dependent factor. As shown in Figure 13, the production rate has the largest effect on the value of BT. As the production rate increases, the breakthrough occurs earlier, leading to the reduction of the BT. This is consistent with often published results.1–10 Increasing fracture to matrix permeability (Kf/Km) also results in lowering the value of BT. When the reservoir height or oil zone thickness increases, it takes longer time for the water to breakthrough. The three observations and results are physically logical. The same methodology was followed to accomplish ANN modeling for critical production rate. The variables in the input layers are oil zone thickness or height (H), fracture permeability to matrix permeability ratio (Kf/Km), density Table 3. Statistical Performance of the ANN-PSO Model vs. ANN-BP Model Parameters ANN-PSO ANN-BP R2 MSE MIPE (%) MAPE (%) 0.99 0.024 3.4 8.9 0.91 0.037 6.3 11.2 Published on behalf of the AIChE DOI 10.1002/aic 1913 Figure 11. MSE vs. number of Epochs to check validation of ANN-BP model. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] difference (qw – qo) and oil viscosity (lo). A satisfactory match (R250.992; MSE 5 0.048) was observed between the estimated values and the actual data when the hybrid ANN model was used. To keep the manuscript in an acceptable size, only the results attained from ANN-PSO and ANN-BP models in the testing stage are presented in Figure 14. Again, superiority of PSO-ANN with respect to the conventional ANN is confirmed. Figure 14 shows the effectiveness of the hybrid smart technique in terms of statistical analysis very well and the results of the training phase do not seem necessary for the purpose of comparison and performance evaluation. To determine the relative significance of the main variables affecting the critical rate, the ANOVA technique was used. The importance weight of the input parameters includes height or thickness of oil zone (37%), permeability ratio (28%), density difference (20%), and viscosity (15%) Although, the main objective of this work was to demonstrate the potential of ANN to learn the functional relationship between input and output parameters affecting fluid displacement during water coning in fractured systems, this work launches a foundation for implementation of NNs in prediction of production performance of heterogeneous reser- Figure 13. Relative effect of important independent variables on BT. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] voirs containing vugs and fractures, which will be a part of our future study. In production engineering, determination of critical production rate and BT appears to be important while designing and/or operating a production well. Having accurate values for these two important parameters is inevitable to obtain efficient oil flow rate and recovery factor in terms of Figure 12. MSE vs. number of Epochs to check validation of ANN-PSO model. Figure 14. Predictive performance of the ANN systems used for estimation of critical rate in the testing stage: (a) ANN-BP and (b) ANN-PSO. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] 1914 DOI 10.1002/aic Published on behalf of the AIChE May 2014 Vol. 60, No. 5 AIChE Journal technical and economic perspectives. Such an investigation and the predictive connectionist modeling study provide vital information for chemical (or process) engineers while working with petroleum engineers on the engineering design of surface and subsurface facilities involved in oil production operations. Given further clarifications, production rate significantly affects the power required for pumps, size of pipeline, design of surface equipment, and safety regulations, resulting in considerable variations in capital and operational expenses. In general, this study covers various practical applications as well as research areas such as production engineering, process engineering and transport phenomena in porous media. Example Given an oil fractured reservoir with the following information. Calculate the BT and the critical production rate. Data: Oil density (qo) 5 0.541 g/cm3; water density (qw) 5 0.995 g/cm3; oil viscosity (lo) 5 1.2 cp; matrix permeability (Km) 5 500 mD; fracture permeability (Kf) 5 2000 mD; oil zone thickness (H) 5 65 m; production rate (q) 5 3500 STB/day; temperature (T) 5 55 C; pressure (P) 5 120 atm. Note: STB stands for the standard barrel. Solution. As discussed in the text, the BT is a function of oil thickness, production rate, and ratio of fracture permeability to matrix permeability. The normalized values of these parameters are obtained using Eq. 7 as follows i 2xmin Þ 215 2ð6520:3Þ Oil thickness (H): X1 5 ðx2ðx ð12020:3Þ 2150:081 max 2xmin Þ Production rate (q): X2 5 2ð350020:5Þ ð800020:5Þ 21  20:125 Ratio of fracture permeability to matrix permeability 2ð422Þ (Kf/Km): X3 5 ð15022Þ 21  20:97 The critical production rate is also dependent on oil thickness, oil viscosity, ratio of fracture permeability to matrix permeability, and density difference. The variables are normalized as below i 2xmin Þ Oil thickness (H): X1 5 ðx2ðx 215 2ð6520:3Þ ð12020:3Þ 2150:081 max 2xmin Þ Conclusions Oil viscosity (lo): X2 5 2ð1:220:4Þ ð3520:4Þ 21  20:95 Ratio of fracture permeability to matrix permeability (Kf/Km): 2ð422Þ X3 5 ð15022Þ 21  20:97 Density difference (qw 2 qo): X4 5 2ð0:45420:01Þ ð0:5520:01Þ 21  0:65 To predict the BT and critical production rate, the above points are placed in the training part. Then we can use the MATLAB toolbox or the MATLAB code written by the authors to compute the output parameters. Using the MATLAB toolbox, the following results are obtained: BT 5 385 days Critical production rate 5 5476 STB/day Practical implications Nowadays, management of petroleum reserves needs the high-technology implements. Using such tools may lower or increase the cost of production and development of oil and gas reservoirs. However, the oil companies are always looking for strong tools with low expenses. Smart techniques such as ANN models have capability to address a number of AIChE Journal May 2014 Vol. 60, No. 5 problems such as critical production rate and BT in production engineering that conventional methods are not capable to explain through an efficient, reasonably priced, and accurate manner. Furthermore, the ANN systems supply a practical solution to the problem of converting basic data of reservoirs and production to important parameters (e.g., critical rate and gas/oil ratio). Petroleum and chemical engineers can benefit from connectionist modeling when sufficient engineering data are not available for design and analysis purposes, but where a large amount of production history is attainable. It should be also noted that shortage of required engineering data might be due to high expenses of well testing, core flooding, and so forth. To be specific, a better understanding of water breakthrough phenomenon and accurate estimation of critical production rate and BT for oil reservoirs offer the petroleum and chemical engineers suitable rules of thumb for optimization of oil production conditions. This helps in minimizing the amount of water in the production and surface facilities, resulting in reduction of cost of separation that is attainable in the case of a successful engineering design. For instance, our study can assist skilled engineers to know the most important parameters affecting the critical rate and BT, leading to regulation of the operating conditions to delay the water breakthrough in the production stream. Taking another example, it is not easy to obtain the exact value for the critical production rate based on field tests. However, reservoir properties (e.g., thickness and permeability) are usually available to determine a safe range for oil production rate. From the modeling and test results, it was found that the classic methods like the approach developed by Muskat and Wyckoffl2,3,11 are not valid where there are highly fractured reservoirs with high oil flow rates as the conventional techniques or empirical equation predict the BT and critical rate 20–25% greater than the actual values. It can be concluded that the NN model can aid in making key decisions on oil formations particularly fractured ones while producing oil. Thus, the outcomes of this research and predictive tools come into sight to be helpful in the design phase of more proficient production processes. This study presents a novel application of the ANN technique to predict water coning in fractured porous systems. The ANN-BP and ANN-PSO models were developed to forecast the output parameters after the input parameters were used to train and test the models. The BT was defined in terms of production rate, oil zone thickness, and fracture permeability to matrix permeability ratio. Based on the results of this study, the following conclusions can be made: 1. Different network structures were examined using BP and PSO to minimize the value of error. For this case, the most successful network architecture included one sigmoid hidden layer and one linear output layer. The hidden layer contained seven neurons. 2. The model predictions were compared with the real data, indicating an acceptable agreement as the absolute maximum error is below 9% for the ANN-PSO model. 3. A high density of vertical fractures characterizes the worst case scenario for a fractured reservoir during oil production due to rapid growth of water coning. This condition is associated with the porous medium whose vertical permeability is much higher than the horizontal permeability, leading to early water breakthrough. Published on behalf of the AIChE DOI 10.1002/aic 1915 4. The performance of the ANN-PSO model in estimating BT is better than that of the conventional NN model that uses only the BP technique. 5. ANN-PSO has the capability to avoid of being stuck in local optima while forecasting the BT as this developed model has both global and local searching abilities. 6. The ANOVA technique confirmed that the production rate has the most significant impact on the BT even in fractured media, while critical production rate is mostly affected by the oil zone thickness and permeabilities ratio. 7. The conventional techniques such as Muskat and Wyckoffl’s model2,3,11 are not suitable to estimate critical oil rate in fractured reservoirs because the values predicted by the classic model are 20–25% higher than the real ones. 8. The proposed NN structure was obtained through trial and error procedure. An alternative technique is required so to combine with the PSO algorithm to optimize the NN structure. Acknowledgments Financial support from the Mitacs Elevate program and the Natural Sciences and Engineering Research Council of Canada (NSERC) is gratefully acknowledged. Notation Acronyms ANN = ANOVA = BP = BT = EOR = FL = GA = ICA = MAPE = MIPE = MSE = NN = OWC = PA = PSO = STB = WOC = WOR = artificial neural network analysis of variance back propagation breakthrough time enhanced oil recovery fuzzy logic genetic algorithm imperialist competitive algorithm maximum absolute percentage error minimum absolute percentage error mean squared error neural network oil–water contact pruning algorithm particle swarm optimization standard barrel water–oil contact water oil ratio Variables b= BT = c1, c2 = d= de = G= gBest = H= itermax = iter = K= Kf = Km = m= n= Pe = Pw = = pBest i q= 1916 fracture aperture, mm or cm breakthrough time, s trust parameters distance between drainage boundary and wellbore, cm or m distance between gas/liquid interface and production well, cm or m number of training samples global best position oil zone thickness, cm or m maximum number of iterations current number of iteration permeability, Darcy or mDarcy fracture permeability, Darcy or mDarcy matrix permeability, Darcy or mDarcy number of output nodes number of samples reservoir pressure, Pa or atm wellbore pressure, Pa or atm individual best position of particle i production rate, cm3/s or m3/s DOI 10.1002/aic R2 = rw = rc = ski = vki = w= wmax = wmin = Xi = xi = xmax = xmin = ypj = yaj = coefficient of determination wellbore radius, cm or m critical radius, cm or m current position of particle i at iteration k velocity vector of particle i at iteration k, m/s or cm/s weight inertia maximum weight inertia minimum weight inertia normalized value of the parameter “xi” output or input variable maximum value of output or input variable minimum value of output or input variable predicted output actual output Greek letters ki lo qo qw /e = = = = = random number corresponded to particle i oil viscosity, cp or kg/ms oil density, g/cm3 or kg/m3 water density, g/cm3 or kg/m3 effective porosity Subscripts c= e= f= i= m= max = min = o= w= critical effective fracture particle i matrix maximum minimum oil well Superscripts a = actual p = predicted. Literature Cited 1. Ahmed T. Reservoir Engineering Handbook, 4th ed. Oxford: Gulf Professional Publishing, imprint Elsevier, 2010. 2. Muskat M. The flow of homogeneous fluids through porous media. New York: McGraw-Hill, 1937. 3. Muskat M. Physical Principles of Oil Production. New York: McGraw-Hill, 1949. 4. Meyer HI, Garder AO. Mechanics of two immiscible fluids in porous media. J Appl Phys. 1954;25(11):1400–1406. 5. Blake JR, Kueera A. Coning in oils reservoirs. Math Sci. 1988;13: 36–47. 6. Elkins LF. Fosterton field-An unusual problem of bottom water coning and volumetric water invasion efficiency. Petroleum Transactions, AIME, SPE (Society of Petroleum Engineers), 1959;216:130–137. 7. Karp JC, Lowe DK, Marusov N. Horizontal barriers for controlling water coning. JPT. 1962;14:783–790. 8. Fortunati F. Water coning at the bottom of the well. Technical note. SPE J. 1962;544:1–9. 9. Smith CR, Pirson SJ. Water coning control in oil well by fluid injection. SPE Reservoir Eng J. 1963;613:314–326. 10. Outmans HD. Effect of coning on clean production rate of well in heterogeneous reservoir. SPE 893. In: 39th Annual Fall Meeting SPE. SPE (Society of Petroleum Engineers), Houston, TX, October 11–14, 1964. 11. Muskat M, Wycoff RD. An approximate theory of water-coning in oil production. Trans AIME. 1935;114(1):144–163. 12. Chaney PE, Noble MD. How to perforate your well to prevent water and gas coning. Oil Gas J. 1956;55:108–115. 13. Chierici GL, Ciucci GM, Pizzi G. A symmetric study of gas and water coning by potentiometric models. JPT. 1964;16(8):923–929. 14. Wheatley MJ. An approximate theory of oil/water coning. SPE 14210. In: 60th Annual Technical Conference, 1985. SPE (Society of Petroleum Engineers), Las Vegas, Nevada, USA. 15. Abass HH, Bass DM. The critical production rate in water coning system. SPE 17311. In: Permian Basin Oil and Gas Recovery Conference. SPE (Society of Petroleum Engineers), Midland, TX, March 10–11, 1988. Published on behalf of the AIChE May 2014 Vol. 60, No. 5 AIChE Journal 16. Guo B, Lee RLH. A simple approach to optimisation of completion interval in oil/water coning systems. SPE 23994. SPE Reservoir Eng J. 1993;8(4):249–255. 17. Schols RS. An experimental formula for the critical oil production rate. Erdoel-Erdgas Z J. 1972;88:6–11. 18. Hïyland LA, Papatzacos P, Skjaeveland SM. Critical rate for water coning: correlation and analytical solution. SPE Reservoir Eng. 1989;4(4):495–502. 19. Ould-amer Y, Chikh S, Naji H. Attenuation of water coning using dual completion technology. J Pet Sci Eng. 2004;45:109–122. 20. Ould-amer Y, Chikh S. Transient behavior of water–oil interface in an upward flow in porous media. J Porous Media. 2003;6(2):1–12. 21. Arthur MG. Fingering and coning of water and gas in homogeneous oil sand. Trans AIME. 1944;145:184–199. 22. Chappelear JE, Hirasaki GJ. A model of oil-water coning for twodimensional, areal reservoir simulation. SPE J. 1976;16(2):65–72. 23. Yang W. An analytical solution to two-phase flow in porous media and its application to water coning. SPE Reservoir Eng J. 1992:1–38. 24. Yang W, Wattenbarger RA. Water coning calculations for vertical and horizontal wells. SPE 22931. In: 66th Annual Technical Conference and Exhibition. SPE (Society of Petroleum Engineers), Dallas, TX, October 6–9, 1991. 25. Recham R, Osisanya SO, Touami M. Effects of water coning on the performance of vertical and horizontal wells—a reservoir simulation study of Hassi R’mel Field, Algeria. In: SPE/Petroleum Society of CIM paper 65506 presented at the SPE/Petroleum Society of CIM International Conference on Horizontal Well Technology. Petroleum Society of Canada, Calgary, Alberta, Canada, November 6–8, 2000. 26. Sobocinski DP, Cornelius AJ. A correlation for predicting water coning time. J Pet Technol. 1965;17(5):594–600. 27. Bournazel C, Jeanson B. Fast water-coning evaluation method. SPE 3628. In: The SPE 46th Annual Fall Meeting. SPE (Society of Petroleum Engineers), New Orleans, October 3–6, 1971. 28. Coleman SB, Clay HB, McCurdy DG, Norris LH III. A new look at predicting gas-well load-up. JPT. 1991;43(3):329–333. 29. Menouar HK, Hakim AA. Water coning and critical rates in vertical and horizontal wells. SPE 29877. In: The SPE Middle East Oil Show. SPE (Society of Petroleum Engineers), Bahrain, March 11–14, 1995. 30. Lee SH, Tung WB. General coning correlations based on mechanistic studies. SPE 20742. In: The SPE Annual Technical Conference and Exhibition. SPE (Society of Petroleum Engineers), New Orleans, LA, September 23–26, 1990. 31. Hagan MT, Demuth HB, Beal M. Neural Network Design. Boston: PWS Publishing Company, 1966. 32. Mohaghegh S. Neural network: what it can do for petroleum engineers. JPT. 1995;47(1):42–55. 33. Hornick K, Stinchcombe M, White H. Multilayer feed forward networks are universal approximators. Neural Networks. 1989;2:359–366. 34. Mohaghegh SD. Recent developments in application of artificial intelligence in petroleum engineering. JPT. 2005;57(4):86–91. 35. Zendehboudi S, Ahmadi MA, James L, Chatzis I. Prediction of condensate-to-gas ratio for retrograde gas condensate reservoirs using artificial neural network with particle swarm optimization. Energy Fuels. 2012;26(6):3432–3447. 36. Brown M, Harris C. Neural Fuzzy Adaptive Modeling and Control, 1st ed. Englewood Cliffs, NJ: Prentice-Hall, 1994. 37. Holland JH. Adaptation in Natural and Artificial Systems, 2nd ed. Cambridge, MA: MIT Press, 1975. 38. Zahedi G, Fazlali AR, Hussein SM, Pazuki GR, Sheikhattar L. Prediction of asphaltene precipitation in crude oil. J Pet Sci Eng. 2009; 68:218–222. 39. Ahmadi MA, Zendehboudi S, Lohi A, Elkamel A, Chatzis I. Reservoir permeability prediction by neural networks combined with hybrid genetic algorithm and particle swarm optimization. J Geophys Prospecting. 2013;61(3):582–598. 40. Perez-Martınez E, Rodrıguez-de la Garza F, PEMEX E&P, Samaniego-Verduzco F. Water coning in naturally fractured carbonate heavy oil reservoir—a simulation study. In: SPE Latin American and Caribbean Petroleum Engineering Conference. SPE (Society of Petroleum Engineers), Mexico City, Mexico, April 16–18, 2012. 41. Beattie DR, Roberts BE. Water coning in naturally fractured gas reservoirs. In: The Gas Technology Conference. SPE (Society of Petroleum Engineers), Calgary, Alberta, Canada, April 28 to May 1, 1996. 42. Saad SD, Darwich TD, Yousri A. Water coning in fractured basement reservoirs. SPE 29808. In: SPE Middle East Oil Show. SPE (Society of Petroleum Engineers), Bahrain, March 11–14, 1995. 43. van Golf-Racht TD, Fernand S. Water-coning in a fractured reservoir. SPE 28572. In: 69th Annual Technical Conference and Exhibi- AIChE Journal May 2014 Vol. 60, No. 5 44. 45. 46. 47. 48. 49. 50. 51. 52. 53. 54. 55. 56. 57. 58. 59. 60. 61. 62. 63. 64. 65. 66. 67. 68. tion. SPE (Society of Petroleum Engineers), New Orleans, LA, September 25–28, 1994. Hidalgo OJ, Brito LE, Garrido R, Munoz J, Paz F, Flamenco F, Aguilera R. Critical oil rates in naturally fractured reservoirs to minimize gas and water coning: case history of a mexican carbonate reservoir. SPE 121743. In: Latin American and Caribbean Petroleum Engineering Conference. SPE (Society of Petroleum Engineers), Cartagena, Colombia, May 31 to June 3, 2009. Sahimi M. Flow and Transport in Porous Media and Fractured Rock, 2nd ed. Germany: Wiley VCH, 1995. Aggour MA, Kandil AA. Experimental study of horizontal well performance in fractured reservoirs with bottom-water drive. Pet Sci Technol. 2001;19(7&8):933–947. Gunning J, Paterson L, Poliak B. Coning in dual completed systems. J Pet Sci Eng. 1999;23:27–39. Bahadori A, Nouri A. Prediction of critical oil rate for bottom water coning in anisotropic and homogeneous formations. J Pet Sci Eng. 2012;82–83:125–129. Khan AR. A scaled model study of water coning. AIME Trans. 1970;249:771–776. Saidi AM. Reservoir Engineering of Fractured Reservoir-Fundamental and Practical Aspects. Paris: Total Edition Press, 1987. Zamonsky G, Lacentre PE, Larreteguy AE. Towards better correlations for water production prediction using sensitivity analysis and numerical simulation models. SPE 94457. In: SPE Europe/EAGE Annual Conference. SPE (Society of Petroleum Engineers), Madrid, Spain, June 13–16, 2005. Zendehboudi S, Chatzis I, Shafiei A, Dusseault MB. Empirical modeling of gravity drainage in fractured porous media. Energy Fuels. 2011;25(3):1229–1241. Zendehboudi S, Rezaei N, Chatzis I. Effects of fracture properties on the behavior of free-fall and controlled gravity drainage processes. J Porous Media. 2012;15(4):343–369. Zendehboudi S, Chatzis I, Mohsenipour AA, Elkamel A. Dimensional analysis and scale-up of immiscible two-phase flow displacement in the fractured porous media under controlled gravity drainage. Energy Fuels. 2011;25(4):1731–1750. Hornik K, Stinchcombe M, White H. Universal approximation of an unknown mapping and its derivatives using multilayer feed forward networks. Neural Networks. 1990;3(5):551–600. Montgomery DC. Introduction to Statistical Quality Control, 3rd ed. New York: John Wiley and Sons, 2008. Montgomery DC, Runger GC. Applied Statistics and Probability for Engineers, Student Solutions Manual, 3rd ed. New York: John Wiley and Sons, 2006. Haykin S. Neural Networks: A Comprehensive Foundation, 2nd ed. NJ: Publisher Prentice Hall PTR Upper Saddle River, 1998. Demuth H, Beale M. Neural Network Toolbox User’s Guide. MATLAB Version 4, The MathWorks, Inc, 2001. Simpson PK. Artificial Neural System. New York: Pergmon Press Elmsford, 1989. Aminzadeh F, Jamshidi M. Soft Computing: Fuzzy Logic, Neural Networks, and Distributed Artificial Intelligence. New York: Prentice-Hall, 2008. Zendehboudi S, Ahmadi MA, Bahadori A, Shafiei A, Babadagli T. A developed smart technique to predict minimum miscible pressure—EOR implications. Can J Chem Eng. 2013;91(7):1325–1337. Shafiei A, Dusseault MB, Zendehboudi S, Chatzis I. A new screening tool for evaluation of steamflooding performance in naturally fractured carbonate reservoirs. Fuel. 2013;108:502–514. Zendehboudi S, Shafiei A, Bahadori A, James L, Lohi A, Elkamel A. Asphaltene precipitation and deposition in oil reservoirs–Technical aspects, experimental and hybrid neural network predictive tools. Chem Eng Res Des. In press. Available at: http://dx.doi.org/10.1016/ j.cherd.2013.08.001 (accessed Oct. 15, 2013). Eberhart RC, Kennedy J. A new optimizer using particle swarm theory. In: Proceedings of IEEE the 6th International Symposium on Micro Machine and Human Science. IEEE (Institute of Electrical and Electronics Engineers), Nagoya, 1995:39–43. Eberhart RC, Simpson PK, Dobbins RW. Computational intelligence PC tools. Academic Press Professional: Boston; 1996. Shi Y, Eberhart RC. A modified particle swarm optimizer. In: Proceedings of IEEE International Conference on Evolutionary Computation, 1988:69–73, IEEE (Institute of Electrical and Electronics Engineers), Anchorage, AK. Kennedy J. The particle swarm: social adaptation of knowledge. In: Proceedings of the 1997 IEEE International Conference on Evolutionary Computation ICEC’97. IEEE (Institute of Electrical and Electronics Engineers), Indianapolis, Indiana, 1997:303–308. Published on behalf of the AIChE DOI 10.1002/aic 1917 Appendix A: Method of ANN-PSO Modeling In several engineering and science problems, it is not easy to find the precise relationship between input variables and target function (output parameter) through statistical, analytical modeling and numerical modeling tools. The problem may raise due to high nonlinearity of the phenomenon (or process) under study as well as lack of enough knowledge about process physics and mathematical methods. On the other hand, the corresponding tests to measure the target parameters might be costly and time consuming. Thus, developing a fast and reliable technique to estimate target function(s) is inevitable. Nowadays, ANN models are recognized as strong predictive tools to forecast important parameters in engineering and science disciplines. The ANN system mainly includes two stages, namely training and testing. It also contains three types of layers, including input, hidden, and output. The layers have various numbers of processing neurons. Each neuron in the input layer is associated to entire neurons in the hidden and output layers through weights. The numbers of input and output neurons are usually determined based on the problem physics, available data, and targets defined by researchers. A fairly large number of real data (including inputs and outputs taken from laboratory tests, pilot plant runs, and real cases) are collected from various sources in the open literature. A part of the data may also be provided by the researchers who are conducting the corresponding research work. It should be noted that the real data are generally gathered from laboratory data, pilot plant data, and real plant data (even accurate modeling outputs are sometimes used). Therefore, it is logical to cover wide ranges of input and output parameters so that the developed model is applicable to the actual cases (or processes) in the world. After collecting the required data, they are split into training and testing parts through a random manner. In general, 70–80% of the data are assigned to the training stage and the rest goes to the testing stage. Then, the training process is performed to find the appropriate relationship between inputs and output(s) through determination of the optimum magnitudes for the weights of the connections in the ANN structure layers. The ANN parameters are obtained at this stage as a reasonable accuracy is achieved in terms of statistical analysis (e.g., R2, MSE, MIPE, MAPE). After the training stage, the data assigned (just inputs, not outputs) for the testing phase are put in the model to predict the target parameter(s). To examine the performance of the ANN system, the predicted outputs are compared with the measured (or real) data. It is clear that a good agreement is noticed if the ANN model parameters are optimized well throughout the training phase. Needless to mention that if one provides the magnitudes of the input parameters for a known case (or process), the model is able to predict the target parameter(s), accurately. In the current study, we followed the same procedure as the total number of 160 data (60 data points from the tests conducted by the authors 1 100 data points taken from the open literature including real world) were used to construct the hybrid ANN model. The training phase was implemented using 120 real data points to find out the most suitable relationship between the target parameters (e.g., BT and critical rate) and input variables such as height (H), the fracture to matrix permeability ratio (Kf/Km), and flow rate (q). After the optimum values for the connection weights were determined at the training stage, the testing process was performed to examine the accuracy of the developed ANN model through comparing the estimated values with the measured (or actual) ones. At this stage, the predictive tool has been built and no more experimental works or 1918 DOI 10.1002/aic measurements are required to predict the outputs such as BT and critical production rate if the magnitudes of input variables are available for an oil reservoir. Although the program code written by the authors has a strong capability to estimate BT and critical production rate for various oil reservoirs, the MATLAB toolbox can be also used for this purpose as follows: 1. First type “nnstart” in the Command Window of MATLAB. 2. This command can be launched to “Neural Network Start.” 3. This MATLAB page includes various wizards (e.g., Input-Output and Curve Fitting Tool, and Pattern Recognition Tool) to solve different kinds of problems. 4. For our case, “Input-Output and Curve Fitting Tool” is selected. 5. Then, a page is opened for selecting input and output data from a file. We can load total data from an Excel file. 6. After loading the data, the MATLAB toolbox randomly takes 70–75% of the data for the training, 15–20% for the testing, and 10–15% for the validation phases. 7. Before starting the training stage, the characteristics of the ANN network (e.g., number of hidden layers, hidden neurons, etc.) should be put in the MATLAB window. We can use the defaults for this part, as well. If a good solution (or accuracy) is not obtained, it is possible to return the page and make changes on the magnitudes of parameters through trial and error procedure until high R2 and low error percentage are attained. In the current work, the PSO technique is used to determine the optimum values of the parameters without the trial and error technique. 8. When the training part is finished, the toolbox gives you the opportunity to plot performance figures and have the statistical parameters to examine the model accuracy. Table B1. A Part of Measured Data Used in ANN-PSO15–19,24,25,40–51 H (cm) 5 10 15 20 25 30 35 40 15 20 25 30 35 40 10 20 30 40 10 20 30 40 10 20 30 Published on behalf of the AIChE q (cm3/s) Kf/Km BT (s) 0.8 0.8 1 0.8 1 0.7 1 0.7 6 6 6 5 5 5 0.42 0.42 0.45 0.6 1.5 1.5 1.5 1.5 7.9 7.9 7.9 15 15 15 15 15 15 15 15 15 15 15 15 15 15 25 25 25 25 25 25 25 25 25 25 25 10 29 47 60 80 110 140 180 9 15 21 28 40 56 6.5 15 33 51 5 10 25 46 2.5 3.8 7.5 May 2014 Vol. 60, No. 5 AIChE Journal 9. Then, the testing and validation stages are run. This part provides the outputs for corresponding input variables. If an acceptable precision is attained, it means that the ANN model is properly built. 10. Now, if just input variables from a certain reservoir are given, the model can predict the BT and critical production rate without performing relevant tests. The method is the same, just new data points are added to the testing part and the program is run. If the new input data are within the range of old previous data, AIChE Journal May 2014 Vol. 60, No. 5 the ANN model estimates the output parameters with high confidence and accuracy. Also, engineers and researchers can use the MATLAB code developed by the authors, upon request. Appendix B: Real Data Table B1 presents a part of the measured data that were used in the ANN-PSO model. Manuscript received Nov. 5, 2012, and revision received Nov. 29, 2013. Published on behalf of the AIChE DOI 10.1002/aic 1919