Development of an effective completion schedule for a petroleum reservoir with strong aquifer to control water production

A reasonable solution, to deal with oil field water problem, is to minimize the amount of water associated with oil production using effective completion lengths. This work presents an effective method to optimize wells’ completion lengths in an oil reservoir with a strong aquifer. The suggested technique is formulated as a constrained optimization problem that defines a NPV objective function and a set of existing field/facility constraints. An effective algorithm translates the completion lengths to connections number in the dynamic simulation model. In this approach, a genetic algorithm (GA), an adaptive version of simulated annealing (ASA) and a particle swarm optimization (PSO) hybridized with polytope technique are applied to maximize NPV. A comparison is given for their performances in a strong water-drive reservoir where the combinatorial effects of wells’ completion lengths (decision variables) should be addressed. Optimizing the lengths of completions leads to an increased production period, total oil production, retarding water breakthrough, reducing total water production, and finally increasing ultimate recovery. The results showed that total oil production by GA, ASA and PSO algorithm is increased by 11.0%, 2.40% and 2.22%, respectively, related to the initial case. Total water productions are decreased by GA, 9.82%, by ASA 2.11%, and by PSO 1.82% relative to the initial schedule. The best performance belongs to the GA algorithm. Moreover, the average watercut of all wells is decreased through the optimization process. Besides, based on the numerical simulation, closing the worst connections with high watercut decreases total water production, and improves oil recovery, maximum well productivity, and NPV (oil–water ratio is increased 18.2%). Most connections are placed in the layers where water coning can occur later (considering near-well-bore permeability) and slightly far from full water zone.


Introduction
Significant amount of water is produced in upstream industries, and it will have increasing trends. Associated water (formation brine) during oil and gas production is the largest source of water. This is brought up from the oil and gas formation during the production of petroleum and natural gas (Kaur et al. 2009). This large amount of water should be treated, desalinated, or injected. Therefore, it requires different surface and subsurface facilities, that imposes lots of costs and environmental problem to the development projects (Ebrahimi and Vilcáez 2018). Based on a study, when the watercut reaches to 75%, the cost of water treatment is drastically increased (Arslan 2005). Treatment of water with considering environmental legislation is highly important in the oil industry (Jiménez et al. 2018). To face with water problem, the best solution would be minimizing the amount of water produced with oil which is in direct relation to completion length (Pintor et al. 2016). Among different sources of associated water from the petroleum reservoirs, the water (water) produced from conning or aquifer movement in the reservoirs are key players. In other words, although oil wells are drilled and completed/perforated in the zones above the water-oil contact, water production from such sources is problematic. One of the efficient strategies is to have an optimized design for completion and perforation lengths in cased-holes and perforated wells. The perforated completions have several advantages relative to open-hole wells. This is since cased-hole completions bring several advantages to an oil field. First, if some layers have issues such as high watercut or sand production, it will be easy to shut-off the worst offending connections. Second, perforated completions provide a handful for selective productions from multiple productive formations. Third, perforated wells can be used for a longer time without causing any problem for the whole production system. The productivity of a perforated well depends on both perforation and formation features, e.g., length, shot density, damage, radius, permeability, and anisotropy. The precise knowledge about the relationship between pressure drop and oil production rate leads to optimum perforation design and optimized hydrocarbon production. The perforation density and length can have effect on the inflow performance and alleviating watercut of wells in the water-driven reservoir. Early water breakthrough and a rapid increase in watercut cause serious problems where the heterogeneous reservoir has high permeability completion intervals. Short production period reduces oil production, premature water breakthrough, high watercut and lower ultimate recovery are results of an unscientific perforation plan.

Completion design modeling and optimization
There are several published studies on well completion lengths modeling and optimization (Cakici et al. 2013;Liang et al. 2019). Landman and Goldthorpe (1991) proposed a mathematical model that couples the Darcy law with a 1D momentum equation, for each perforation, to optimize the number of perforations for a horizontal well of a simplified infinite reservoir neglecting the skin factor causing by formation damage. Boyun and Lee (1993) proposed an analytical approach to optimize wellbore penetration into a pay zone to maximize critical oil rate. The unstable water cone occurs when the gradient pressure of water contact is lower than the vertical gradient pressure which can highly affect the total oil recovery factor. They concluded that the best wellbore penetration into an oil zone is one-third the total pay-zone thickness; Yildiz (2000) studied the productivity of selectivity perforated vertical wells by considering nonuniform perforation size, the random phasing angle, length and formation damage around perforations. They presented a new mathematical model controlling the impact of gas or water coning in vertical wells. Their study revealed that capillary fluid level is affected by a simultaneous twophase flow in the formation. Zhou et al. (2002) suggested a variable-density perforation approach to control production profile for horizontal wells. Yeten et al. (2003) optimized well placement, wells' completion design and well controls using a genetic algorithm. Yeten et al. (2004) applied a gradient-based optimization method combined with a reservoir simulator to optimize the performance of intelligent wells in several geological model realizations. According to their work, the optimal number of perforations has an important effect on horizontal wells. Zhou (2007) introduced a novel method to optimize the number of perforations along the wellbore. Different perforation optimization procedures have been presented to optimize the number of perforations in horizontal wells by using different optimization algorithms, while the main focuses of their studies were homogenous models (Hua et al. 2010;Wang et al. 2011). By deriving a reservoir model coupled with a wellbore flow model on the basis of inflow controlling theory, Wang et al. (2010) presented a model to optimize the perforation scheme for finding the global optimum of an objective function for recovery. Their study has proposed to use a low perforation density where permeability is high and a high perforation density where permeability is low. Wang et al. (2011) coupled reservoir and wellbore model to calculate a perforation optimization model using genetic algorithm, but their study focuses on the near-well heterogeneity. However, their model cannot be applied in heterogeneous reservoirs with distributed heterogeneous permeability in the whole drainage area. Haijing et al. (2012) improved the modeling by considering near-wellbore permeability heterogeneity. They proposed an optimal perforation scheme to determine inflow profiles of horizontal wells in anisotropic parallelepiped reservoirs through an enhanced assessment approach of permeability heterogeneity for horizontal wells. Their results showed that to achieve optimal productivity, the number of perforations should be higher in the region with low permeability and be lower in the region with high permeability. Non-uniform inflow can occur in a high permeability reservoir in a horizontal well because of the variation in pressure drawdown along the well. A similar model, presented by Pang et al. (2013), considered near-wellbore inflow and wellbore conduit flow on the basis of a reservoir numerical simulation and wellbore dynamics. They used such tools to optimize perforation density for perforation design and optimization of long horizontal wells. Xu et al. (2013) published a numerical model to take the impacts of the number of perforations on the productivity of gas well into account, where the objective function is the vertical well production length, and the decision-making parameter were the perforation distribution. Their results illustrated that inflow rate profiles are successfully increased by the optimal perforation position and perforation density distribution. There are also published works about mathematical modeling of other aspects of well completion for example hydraulic fracturing , , (Zhang and Emami-Meybodi 2020a, b), (Zhang and Emami-Meybodi 2020a, b). (Zhang and Emami-Meybodi 2020a, b) and (Zhang and Emami-Meybodi 2022). In these types of studies, the specifications and configurational properties of hydraulic fractures have been in focus instead of perforations' specifications. Khoshneshin and Sadeghnejad (2018) used two heuristic optimization methods (artificial bee colony and PSO) to determine the optimum well placement and find the best completion design in an Iranian oil field. Their study has divided into two periods of the field life, i.e., the first one is infill drilling stage and the second one is the well location optimization. To find the best completion layer, they coupled geological uncertainties and practical constraints with the finding optimum infill drilling placement problem. Based on their result, artificial bee colony algorithm has a better performance than PSO. Carpenter (2020) coupled both geological and geo-mechanical models with production and pumping history to optimize well completion design in unconventional reservoirs to achieve higher oil recovery factor. Zeynolabedini and Assareh (2021) presented effective design of well completion length in a downhole water loop and assessed different completion characteristics parameters.

Contributions of this study
When a reservoir has an active aquifer and the rate of water production in other wells of the reservoir is high, optimization of perforation length is posing a challenge due to the importance of combinatorial effects. In such a case, the length of perforations should be optimized during the lifecycle of the reservoir. In this way, not only the watercut can be controlled, but also the maximized oil can be recovered. One of the important parts of perforation optimization is planning for the length of perforation and previous studies have not focused properly on the length of perforation where the reservoir has strong water derived using reservoir-wide models. In other words, this study tries to show that just by optimization of perforation length, total oil production and oil-water ratio is increased without causing additional cost such as water treatment or remediation of perforations. In fact, this analysis gives a vision for reservoir engineers to focus on the optimization of perforation length if the production mechanism of a reservoir is water-drive. Note that shutting the worst offending connection necessitates well shut in and workover operation. There is not enough study in the literature to evaluate this kind of optimization problem. Herein, this work tries to optimize the length of perforation in an oil reservoir without considering other parameters of perforations such as phase angle, borehole diameter and length of perforated spacing. This is since the length of perforation is a critical subject for high watercut. The phasing and density of perforation are defined with, connection factor, an abstract parameter for transmissibility of simulation block to well flow. The objective function is the net present value (NPV) with a discount rate applied to take the long-term and short-term gains into account. The ASA, GA and PSO hybridized with polytope are used as global search techniques, to help the optimization escape from the local minima. In fact, authors try to show that optimized length of perforation in a reservoir with strong aquifer can decrease costs of water treatment and well-shutting besides of avoiding environmental issues. In the next section, the method is explained in detail. Later, the case study results are presented to show a comparison between these algorithms and how much they can improve the total production and lower watercut.

Production optimization method
In this work, the total oil and water productions are considered in NPV for optimization. The decision variables are the lengths of perforations (connections) of each well in an oil reservoir in which there is a strong aquifer. In fact, the optimized lengths of perforations are obtained at an initial state through optimizing of both total oil and water production in NPV for production life cycle using a reservoir model. Besides, some practical constraints are considered in the model such as minimum rate of each well, a minimum limit for the length of perforation and inactive cells due to the geometry of the reservoir. The NPV, J, addresses both the long-term and short-term gains through using a discount rate and a certain reference time denominator can impact the total oil production during a certain time. The NPV is based on the cumulative oil and water production over a fixed time horizon introduced by van Essen et al. (2009) and can be defined as: K indicates the total number of time steps k and Δt k the time interval of time step k in days. The term b is the discount rate for a specific reference time t . The water injection cost is defined as: The produced water production cost is defined as: The produced oil production benefit is calculated from: where r o , r wp and r winj are the oil price, the water production, and the water injection costs are in $/m 3 , respectively. The terms q winj.k , q wp.k and q o.k represent the total flow rate of injected water, produced water, and produced oil at time step k in m 3 . N prod and N inj are the number of production and injection wells, respectively. The optimization problem with the objective function can be formulated as follows: which is subjected to: and: where u and x are the control vector (input vector) and the state vector (grid block pressures and saturations), respectively. g is a vector-valued function presenting the system equations, x 0 is a vector of the initial conditions of the reservoir. A colon in a subscript indicates a range, e.g., u 1∶K = ( u 1 ,u 2 ,…,u k ). The vector of inequality constraints g relates to the well's limitations. The mass conversation equation and Darcy's law are combined to create the governing equations being the water-oil material balance. The effects of capillary pressure and molecular diffusion are neglected, and the problem is isothermal.
ASA, GA and PSO are applied to optimize the defined problem in this study. We also make a comparison among them to regarding their performances on this optimization problem. The adaptive version of simulated annealing which uses an advanced annealing schedule and improved generation function of random numbers is implemented. The flowchart of the implemented ASA is shown in the work of Azamipour et al. (2017) in more details.
Simulated annealing, introduced by Kirkpatrick (1984), is based on the physical process of metallurgical annealing which is molecules of hot metal are organized in a way being near to global energy minimum when the temperature falls. The metropolis algorithm helps SA to escape local optimum traps through a Boltzmann probability function. In fact, if the new objective function is decreased, the set of optimization parameters will be accepted automatically; otherwise, if a random number is less than the probability created by the Boltzmann function, it is accepted for the next iteration.
The temperature T is a global time-varying parameter and the objective functions J(u k1 ) and J(u k2 ) are parameters of Boltzmann probability function by which the SA algorithm can escape from a local minimum which is worse than the global one since the probability function P should be positive even when J(u k2 ) is larger than J(u k1 ) . If the difference between two consecutive objective functions (delta ΔJ ) or delta is positive but the metropolis function accepts the solution, the temperature based on a formulation proposed by Ingber (1996) would be adjusted. In these relations, T 0i is the initial temperature and, n i and m i can be considered free parameters that help tuning ASA for specific problems and D is the dimension of the space problem. Otherwise, a new solution is generated. Finally, when the final temperature is reached, the algorithm would be terminated.
Genetic algorithms are one of the evolutionary approaches originated from natural evolution, for example, mutation, selection and crossover and can optimize complicated objective functions (Haupt and Haupt 2004). The population consists of several arbitrarily generated strings. In fact, these strings are a set of randomly produced practical solutions. The strings are individual solutions and are evaluated by a fitness function, then their fitness values are sorted. The nearness of each solution is evaluated to the optimal solution in the fitness or evaluation function and allocates each chromosome a "rank" applied to calculate the crossover probability of each chromosome. In a crossover process, two randomly chromosomes with the highest rank, called parents, are selected, and integrated to create offspring or children to exchange genetic information. The crossover rate is the ratio of the number of children in each generation to the population size. It adjusts the expected number of chromosomes used in the crossover operation. Different types of crossover have been introduced such as single point, double point, multi-point and scattered (Koza 1992). The mutation process not only helps GA escape from local optima but also brings back any lost genetic values when the population converges too fast (Ghaedi et al. 2014). The mutation rate is the percentage of the number of genes whose values are flipped in each generation to the population size. In addition, elitism is the process in which some percentages of highest ranked individuals without any modification are transformed to the next generation to replace some lowest-ranked individuals. Finally, the next generation will be created, and the cycle is repeated to reach the stopping criteria likes the difference between two subsequent objective functions or the maximum number of generations. The flowchart of the implemented GA is presented in the work of Azamipour et al. (2018).
The particle swarm optimization algorithm, introduced by Eberhart and Kennedy (1995), is inspired by the concept of movement of animal groups, such as birds and fishes, and it is capable of solving complicated mathematics problems and optimizing nonlinear continuous functions (de Moura Meneses et al. 2009). The aim of flying a group of birds is to find a place where the accessibility of food is maximized, and predators' threat is minimized. As matter of fact, the flock should search and assess several places by applying some surviving criteria to maximize the survival situations of its population. Kennedy and Eberhart (1995) introduced a classical version of particle swarm optimization to find the optimal solution by studying the social behavior of birds' flocks which help them to increase their chance of survival by identifying a safe place to land. Researchers have proposed different versions of PSO algorithm such as the constriction factor weight (Eberhart and Shi 2000), the dynamic inertia and maximum velocity reduction and the linear-decreasing inertia weight (Shi and Eberhart 1999). The main goal of optimization is to maximize or minimize function f (X) which is represented by the position vector X = [x 1, x 2 , x 3 , … , x n ] in PSO algorithm. The position vector X denotes a variable model and it is dimensions' vector as the number of variables in a problem. Moreover, assessing the objective function f (X) is performed by evaluating how good or bad a position is based on landing points and survival criteria. The position vector is updated according to the following equations: and where P is the total particles of the swarm, i = 1, 2 … P is a component of the vector, j = 1, 2, … , n is the dimension of equations. There are three different terms in Eq. 11 representing a particle's movement in an iteration. The first term is a product between parameter w and particle's previous velocity which is why it shows a particles' previous motion into the current one. The parameter w is the inertia weight constant and a positive constant value for the classical PSO version. This parameter is applied as a balancing between global search, called exploration, and local search knows as exploitation. For instance, if w = 1 , the particle's motion is governed by its previous motion in a similar direction. If , more regions in the searching domain are explored by the swarm which increases the chance of finding a global optimum, however, the optimization procedure will be more time consuming . The second term of Eq. 11, being the individual cognition term, is determined by the difference between it is the particle's own best position, pbest ij , and its current position X t ij . When the current position of the particle is far from pbest ij position, which means the amount of pbest ij − X t ij get larger, the particle will be attracted to its best own position. The parameter c 1 , known as an individual-cognition parameter, is a positive constant, also it denoted the importance of particle's own previous experiences. Equation 11 uses r 1 , a random value with [0,1] range, to avoid premature convergences, increasing the most probable optima. The social learning term is the third term by which all particles of the swarm can share the information of best point attained nonetheless of which particle had found it, e.g., gbest j . The difference gbest j − X t ij represents the amount of attraction of particles to the best point until found at some t iteration. The parameter c 2 , known as the social learning parameter, shows the importance of the global learning of the swarm. r 2 also is a random value like r 1 . Equation 12 updates the particle's position.
Finally, the polytope approach speeds up the optimization methods by finding the best solution among the previous iterations and selecting the best one for the next guess. The polytope algorithm, being a hill-climbing algorithm, does not require derivative information. This algorithm, presented by Nelder and Mead (1965), tries to find the direction that increases the value of the objective function.
To simulate the water movement process, a coupled system of nonlinear equations system for oil and water material balance for each simulation cell should be solved.
where pressure ( P ) and water saturation ( S w ) are used as the primary variables to solve these two equations for all the simulation blocks ( x ). In this relation, ij is the block pair transmissibility between central block, i , and neighbor block, j . Kr o and Kr w are oil and water relative permeability and o and w are oil and water viscosity, respectively. PV i is the pore volume of the grid block. B o and B w are the oil and water formation volume factor correspondingly. The pressure potential difference for the oil phase (driving force) is defined as: The same is valid for water: The main aim of this approach is to find the optimum length of perforation. To do that the top of perforation for all the wells are considered constant, and the lowest one is optimized by considering some practical limitations. At first, the largest possible length of perforation for all the wells are applied. Then, the model is simulated coupled with the optimization algorithm. The reservoir flow is determined by commercial reservoir numerical software by which reservoir static information and dynamics can be captured, utilized, and predicted. The new solution goes to the polytope algorithm to compare it with the older solution to select the best one. Finally, when the algorithm reaches the stopping criteria, the process will be finished.
In this research, ASA, GA and PSO algorithms are used as the optimization method in separate cases to compare the performance of each of them on this problem. In the section of initializing, all the perforations and optimization parameters are inputs, and the simulator is run, and the results of the simulation are read, which contains yearly oil and water production, therefore the objective function is calculated. In the next section, being the main section of optimization, new perforations are generated based on the best solution of the NPV by the random generator. The simulator is run again, and the objective function is determined. The two recent objective functions are compared, and ASA algorithm decides which objective function should be used. If the new objective function is lower than the previous one, then the new objective function will be accepted. Otherwise, the metropolis algorithm decides to keep the objective function or not. Then, polytope search finds the best solution among all the solutions and select ΔΦ w = Φ wj − Φ wi = P wj − P wi − 1 2 wj + wi gΔd ij the best NPV by which the best perforations will be used in the stage of generating new perforations. Finally, if the temperature is down enough, the optimization will be stopped, and the result will be shown. The flowchart of optimization-simulation coupling using GA is started by initializing optimization and simulation parameters and first guess for the perforation. The NPV is calculated after the simulator is run and the oil and water produced are calculated annually. The perforations of each well are turned into the chromosome and new perforations are generated and are used as initial data for the simulator. After the simulation, the NPV is calculated and stored and ranked as a fitness function. The polytope algorithm helps the optimizations to initialize based on the best solution that has been found. If the generation reaches the criteria to terminate the search, the result will be presented. Figure 1 shows the flowchart of optimization-simulation coupling using PSO algorithm. In the first step of optimization, initial perforations, PSO parameters randomly, particle position and velocity are initialized. The simulator is used to calculate the NPV based on determining the amount of oil and water produced annually. The perforations are considered as positions of particles and new perforations are produced by the random generator of PSO. Again, the objective function is calculated after simulating with new perforations. The fitness function for each particle is evaluated and local and global solutions are computed. The velocities and swarm positions are updated, and the particles move with new PSO velocity. The new searching direction is calculated. Polytope search assists the algorithm  to benefit a better guess for generating new perforations. If the objective function is satisfied or the optimization is reached termination criteria, the optimization is stopped and the result are displayed; otherwise, the optimization process is repeated.

Reservoir simulation model
The optimization procedure was evaluated using a 3D oil reservoir model with two-phase flow of water and oil. The model consists of 10,800 (30 × 90 × 40) simulation grid cells and nine production wells with average properties presented in Table 1. The supplementary material file for this article contains the reservoir model information and 3D distribution. The 3D grid is presented in Figure 1 in supplementary material. The initial and bubble point pressures are 1645 and 1397 psi, respectively. The porosity and permeability maps are presented in Figure 2 in supplementary material. The anisotropy ratio in x, y and z direction equals unity. The model has four Fetkovich aquifers at the edge of the reservoir that the main properties of the aquifers are presented in Table 2.
Fetkovich is a successful and apparently accurate aquifer model which uses an equation describing flowrate from aquifer and an aquifer material balance. The first parameter of Fetkovich aquifer model is a pseudo-steady state productivity index. The other crucial parameter describes a relationship based on the material balance between the cumulative influx and aquifer pressure. Since smaller aquifers may reach a pseudo-steady state condition rapidly, the Fetkovich model is best appropriate for them. In this approach, the influx of water is modeled equivalent to the flow of oil from a reservoir into a well. The aquifer inflow is modeled by the equation: where Q ai is the inflow rate from the aquifer to the connecting grid block i . W ai is the cumulative influx from the aquifer to grid block i . J is the specified productivity index of the aquifer. i is the area fraction for connection to grid block i . The p a , p c and p i are the pressure in the aquifer at time t, the capillary pressure, and the water pressure in a connecting grid block i , respectively. The is water density in the aquifer. d i is the depth of grid block i and d a is the datum depth of the aquifer. The area fraction for each grid block connection is given by: In this relation, A i is the area of block face communicating with the aquifer and m i is an aquifer influx coefficient multiplier. The pressure response in the aquifer is calculated by the material balance equation: In this relation, W a is cumulative total influx from the aquifer, C t is the total (rock and water) compressibility of the aquifer. V w0 is the initial volume of water in the aquifer and p a0 is the initial pressure of water in the aquifer.
The aquifer time constant, having the dimension of time, and the productivity index has great roles on the aquifer performance. The aquifer time constant is given by: If the connecting grid blocks have uniform reservoir pressure, by integrating Eqs. 17 and 19, the average influx rate over the time interval Δt can be defined as: Which is the form applied by the simulator to calculate the influx rates. The aquifer's cumulative total influx is increased, and its pressure is updated using Eq. 19 at the end of each time step. Production is controlled by the production target rate. The length of the numerical simulation time is 27,375 days. All the remaining fluid and geological properties being used in the model are given in Table 3. Figure 3 and Figure 4 in the supplementary material represent the movement of oil saturation in different cross sections of the reservoir model with its wells and their perforations before starting simulation and after that. These figures show the upward movement of water saturation in a cross section of the reservoir model before and after simulation which represents the effects of production and the power of aquifer to push up water. Hence, it is important to choose optimal perforations intervals to avoid high watercut production and postpone abandoning wells as far as possible.

Results and discussion
The optimized scenarios are chosen between 450 scenarios created in the ASA, GA and PSO algorithms. The NPV is written in C# programming language with parameters r o = 38.0 $/bbl and r w =5.0 $/bbl. The discount factor is 0.25. The minimum oil production was defined as 10.00 $/bbl as a practical constraint in the optimization problem; if the oil production rate is under this limit, the well should be shut in. There are 40 layers in the model in which some cells in different layers are inactive. The wells begin to produce oil with a maximum length of perforation. The lowest perforations of wells (connection positions) are selected as decision variables, and the top perforation of all wells remained fixed at the first layer. The minimum length of perforation for all wells are considered as in 10 th layer as another practical constraint which means that the minimum number of connections is 10. The maximum temperature, initial temperature and quench factor are 200, 0.1, and 1.0, respectively, in ASA algorithm. In the GA algorithm, the probabilities of crossover and mutation are 0.70 and 0.09, respectively, and the percentage of elitism is 5%. The type of crossover and selection are double point and fitness proportionate, correspondingly. The population size which is the number of simulations runs in each generation was set to 30. In the PSO algorithm, the number of particles, max epochs, inertia weight, local weight and global weight are 30, 20, 0.729, 1.49445 and 1.49445, correspondingly. Figure 2 describes all the steps of the proposed optimization procedure. It is started with inputting the reservoir model and initial parameters of the optimization algorithm. The algorithm reads initial perforations of each well which are input into the reservoir model and stores them into a list. For the first iteration of optimization, the simulator is called, and the objective function is computed by calculating the amount of oil and water produced annually during the time of the simulation. Otherwise, the random generator of the optimization algorithm creates new perforations for each well based on the best solution that is achieved. If new perforations meet the constraints, they will be accepted and will be input (translated to the number of connections) into the simulator, otherwise, the random generator generates new perforation. The simulator is used to compute a new objective function with new perforations of wells. Then, the optimization algorithm, i.e., ASA, GA or PSO is applied to optimize the solution. Polytope algorithm searches among all the solutions and finds the best NPV. The algorithm selects the perforations of the best NPV as an initial guess to use in the random generator to create new perforations. If the algorithm reaches termination criteria, the optimal solution will be shown, and the optimization process will be stopped. The results of the optimization process are shown in Tables 4 and 5. In Table 4, the initial and best perforations of all wells for each algorithm, the total water and oil production and watercut for each well are presented. As it is shown in Table 5, GA algorithm has better performance than ASA and PSO algorithms and initial case based on the values of the objective function which are 336,832, 332,559, 332, 498 and 327,266 $, respectively. Total oil production by GA, ASA and PSO algorithm is increased by 11.0%, 2.40% and 2.22%, respectively, related to the initial case. Total water productions are by GA, 9.82%, by ASA 2.11%, and by PSO 1.82% decreased related to the initial case. As another result, the average watercut of all wells also are decreased through the optimization process; however, the best performance belongs to the GA algorithm. Figure 3 (top) represents total oil production versus time for ASA, GA and PSO algorithm and initial case. As it can be seen in the plot, GA (densely dashed) has a higher performance. The total oil productions optimized by both ASA and PSO algorithm have similar profiles which are higher than the base case (base completion design). In Fig. 3 (bottom), a comparison of total water cut between ASA, GA, PSO algorithm and the initial case is displayed. The total watercut of genetic algorithm (densely dashed) is lower than the one in other two cases; however, there are some breaks in the profile that is occurred when a well with a high watercut is shut in. Therefore, the amount of watercut is dropped suddenly. PSO and ASA have the same profiles which are in the interval of GA and the base case. The reason that the drops in the profiles do not occur simultaneously is that well/wells which are closed are not similar.
The performance of optimization algorithms (i.e., ASA, GA and PSO) and the initial (base) case for two wells of P1 and P5 are presented in Fig. 4 For well P1, not only oil production rate of genetic algorithm is higher than other cases, but also watercut of optimization with genetic algorithm is Fig. 2 Flowchart of the optimization process considering translation of perforation length to number of connections highly lower than others. The initial case has the worst performance among other cases. Due to geological properties, well P5 is prone to produce high water, therefore, the reasonable scenario is the lower amount of oil production rate. Watercut of GA is considerably lower than other methods during production since the oil production rate of optimization using GA is low. This plot clearly shows that using the optimization approach helps watercut of wells which are susceptible to produce a high amount of water is controlled by optimizing their perforation lengths (connection intervals). Figure 5 to 7 in the supplementary material displays 3D distributions of the reservoir model which represents perforations for wells in the genetic algorithm approach. Most perforations are placed in the layers where water coning can occur later and slightly far from full water saturated zone. Hence, watercut will be lower during production and the reservoir is able to produce more oil. Figure 5 shows the objective function in each iteration in GA optimization as an example. Since each alteration in the length of the perforation well has a significant effect on the objective function, the polytope algorithm cannot create an ascending trend on the objective function and that is the reason that the objective function fluctuated highly during the optimization process. Across section of the reservoir showing the rise of the water in the reservoir is presented in Figure 8 in the supplementary materials.
To evaluate the performance of optimized solutions, two cases are performed (Case 1 and Case 2) in which initial guesses for the perforations are determined based on the performance of each well. These two cases are just arbitrary unoptimized possible solutions from middle GA generations. They have been selected since they have high oil and low water production (even though they are unoptimized). The results are shown in Tables 6 and 7. Although in these two cases, the total oil, and water production and average watercut are better compared relative to the optimized solution in GA algorithm, due to discount factor and economical consideration, the objective functions of these cases are lower than the objective function of GA algorithm. Hence, it verifies that, however, there are cases with lower total water production and/or higher oil production, the optimized solution has a higher objective function in which is considered both long and short gains.
Another scenario, called GA*, was designed to find the optimized solution of GA algorithm in which a constraint on the watercut is applied to the wells. If a well violates the limitation value for watercut, 0.98%, the worst offending connection, which has the highest watercut in the well, will be closed. Even the entire well will be shut in for a case of extremely high watercut. This remedy has a significant effect on the reduction of total water production and watercut and the increment of total oil production. For GA scenario, the total production is 19.92 MMstb, oil total production is 90.47 MMstb and oil-water ratio is 0.22. For GA* scenario, the total production is 21.26 MMstb, oil total production is 81.63 MMstb and oil-water ratio is 0.26. The results are shown in Table 8. Therefore, it is so important that watercut should be monitored and controlled during production. There are however several limitations for this scenario. In this scenario, some costs are not considered such as costs of a workover operation to shut in the offending perforation and costs of shutting a well. In fact, by considering economic factors, a reservoir engineer should balance the costs of shutting in well/wells and the costs of oil and water production. Including these factors and finding the optimized time to shut in well/, wells can step further in future research to provide a comprehensive and practical optimization problem.
For a future study, it is suggested to use different properties and parameters such as costs of well shutting period, well perforations remediation, water treatment equipment and changes of oil-water contact. Also, it is important to investigate this problem in reservoirs that have both gas and water coning to optimize length of perforation. In this work there is no specific solution for completion design for a reservoir with expanding gas cap. The completion schedule developed in this work is for the reservoirs with strong aquifer. The optimized completion scenario in this work may not be suitable for a reservoir under water flooding process, however the basis of optimization algorithm can be used for such cases. Although production control parameters have important effects on the watercut, assuming in the reservoir with strong aquifer, optimization of perforation intervals in current and future wells can be a long-term approach that helps the reservoir to produce more oil with declining watercut. Authors showed the importance of optimization of perforation length and emphasized on the operations before starting production. In many oil fields, reservoir engineers use liquid production rates or bottom hole flowing pressure which are not solely a perfect solution to reduce water, but if perforation length optimization is coupled with production control parameters, it can be a practical solution.

Conclusions
This work presents a comparison between three optimization techniques, adaptive simulated annealing, genetic algorithm, and particle swarm optimization, to optimize the lengths of perforations in the production wells of an oil reservoir model in which there is a strong aquifer. The aim of optimization is finding the best perforation lengths of the wells in the water-drive reservoir for its life cycle to produce oil as much as possible, control watercut and water production. Based on this study, it is concluded that: • Optimization of the completion lengths should be performed prior to production. When a reservoir has a strong aquifer, determining the optimized lengths of perforation is important to avoid or postpone the additional costs of workover or water treatment. • Optimized lengths of perforation lead to increasing low-water-cut, smaller water production periods, total oil production, retarding the water breakthrough, and finally increasing the ultimate recovery. • Based on the result of optimization, GA algorithm has a better performance in comparison with ASA and PSO algorithms. However, their results are approximately close to each other. • Comparison with un-optimized perforation scenarios shows that it is possible to have cases with higher oil production than the optimized perforation scenario, however Objective Function ($)