## Abstract

In this work we present a methodology for optimal management of brownfields that is illustrated on a real field. The approach does not depend on the particular reservoir flow simulator used, although streamline-derived information is leveraged to accelerate the optimization. The method allows one to include (non-linear) constraints (e.g. a recovery factor larger than a given baseline value), which are very often challenging to address with optimization tools. We rely on derivative-free optimization coupled with the filter method for non-linear constraints, although the methodology can also be combined with approaches that utilize exact/approximate gradients. Performance in terms of wall-clock time can be improved further if distributed-computing resources are available (the method is amenable to parallel implementation). The methodology is showcased using a real field in west Siberia where net present value (NPV) is maximized subject to a constraint for the recovery factor. The optimization variables represent a discrete time series for well bottom-hole pressure over a fraction of the production time frame. An increase in NPV of 7.9% is obtained with respect to an existing baseline. The optimization methods studied include local optimization algorithms (e.g. generalized pattern search) and global search procedures (e.g. particle swarm optimization). The controls for one injection well in the real field were actually modified according to the solution determined in this work. The results obtained suggest improvement for most economic scenarios.

The importance of brownfields resurges in periods of relatively low values of oil price. The management of this kind of field is principally based on the periodical setting of well controls (this may also include shutting wells in and switching their functionality, e.g. from being a production well to an injection well). Expensive capital investments such as the drilling of new wells are not contemplated in brownfields very often, and, for this reason, these fields can be a very attractive component of a field portfolio of an oil and gas company.

The modification of well settings in brownfields can be an arduous task that cannot be performed with arbitrary frequency. Thus, it can be very important to elaborate optimized schedules of well settings that lead to efficient production of the field. However, this can be a rather difficult enterprise because of the frequently large number of control variables associated with the wells in a brownfield, the elevated computational cost of the numerical models used (which can render optimization problems challenging even for a moderate number of variables), and the presence of multiple and complicated operational constraints (finding feasible controls that satisfy all of these constraints is usually not trivial). On many occasions these schedules are determined using available rules of thumb and educated guesses, which, due to their heuristic foundations, may lack accuracy and robustness, and even fail and not provide acceptable solutions.

Therefore, in search of techniques more solid than ad hoc solutions, brownfield management has been approached by means of formal optimization and accurate simulation. In the optimization problems commonly formulated for brownfields, one aims at maximizing a performance metric such as net present value (NPV), and the optimization variables are sequences of well settings (e.g. well bottom-hole pressure (BHP) values or rates) that represent how each well is controlled as a function of time. Operational (non-linear) constraints can be given, for example, by a minimum recovery factor (that may correspond to a baseline production plan) or through pre-specified bounds for field quantities (e.g. maximum field injection rate). The cost function in these optimization problems has been often observed to be relatively smooth (unlike in well location optimization problems, where the cost function may be markedly non-linear due to significant spatial heterogeneities in the reservoir rock properties), and, in cases of moderate or high number of variables, different optimized solutions with similar cost-function values can be found (e.g. see Echeverría Ciaurri *et al*. 2011*a*; Van Essen *et al.* 2011). Consistent with that, the optimization techniques of choice in the context of brownfield management are in most cases of a local nature (in contrast to global search procedures that may be used, at the expense of additional computation cost, in those cases where the optimization cost function might present a higher degree of non-linearity than expected; global search methods are less frequently employed in production optimization than local algorithms; the reader may consult Harding *et al.* 1998; Echeverría Ciaurri *et al*. 2011*a* for examples of applications of global search methods in production optimization).

Whenever derivative information can be computed efficiently, local optimization is addressed by means of gradient-based optimization methods (e.g. see Nocedal & Wright 2006). Adjoint formulations (e.g. see Sarma *et al.* 2008; Van Essen *et al.* 2010) can be considered in simulation-based optimization for a rapid gradient computation of the cost function, provided one has complete access to and detailed knowledge of the simulator used. This is not the case in many situations in practice (e.g. when the simulation engine relies significantly on commercial software) and consequently either gradients are approximated (e.g. see Su & Oliver 2010; Dehdari & Oliver 2012; Zhao *et al.* 2013 for schemes relevant to production optimization) or not even computed. Derivative-free optimization (e.g. see Kolda *et al.* 2003; Conn *et al.* 2009; Kramer *et al.* 2011) provides a fairly efficient alternative in practice for a broad spectrum of optimization problems in oilfield operations (in this respect, the reader may consult, for example, Echeverría Ciaurri *et al*. 2011*b*). It is worthwhile noticing that global search methods proceed, in general, in a gradient-free fashion. Unlike in adjoint-based optimization techniques, the computational cost of most derivative-free methods is larger the higher the number of optimization variables the problem has. However, many gradient-free algorithms are amenable to a distributed-computing implementation, and that may yield acceptable performance in terms of wall-clock time. Another family of efficient alternatives to adjoint formulations leverages techniques that borrow concepts from streamline-based simulations (the reader may consult Batycky & Thiele 2006; Datta-Gupta *et al.* 2010; Voznyuk *et al.* 2014; Wen *et al.* 2014 for some examples of interest in brownfield management). These techniques are essentially devised to optimize well rates and, in principle, this may not be the best option in real scenarios because a considerable number of fields are controlled by means of procedures that do not provide a constant rate for relatively long periods of time (e.g. injection wells controlled through a selection of wellhead choke sizes). Well controls based on piecewise constant BHP values are much closer to existing field implementations and, since in this work we eventually apply the results of the optimization to an actual brownfield, these controls are the ones that will be considered here (however, the method proposed in this research can also be used with well rate controls).

Derivative-free optimization approaches, although applicable to a wide variety of problems, may render formulations that, if there are a relatively high number of variables, can be extremely time-consuming to solve (e.g. when models that integrate reservoir and surface facilities are considered). This might not be an attractive option in many practical scenarios where solutions are required in a short time frame (e.g. overnight or in few days). Hence, derivative-free optimization is very often combined with surrogate modelling (e.g. see Koziel *et al.* 2011, and especially Jansen & Durlofsky 2017, for applications of reduced-order models in production optimization) and with mechanisms that alleviate the computational cost associated with an elevated number of optimization variables. In this work we reduce the number of variables in each iteration of the optimization algorithm by selecting the subset of wells which contribution to a measure related, for example, to the cost function or to the constraints is the highest for all subsets with the same number of wells. It should be stressed that this selection takes into account information related to streamline-based simulation, is applicable to general well controls (not necessarily well rates) and is updated as the optimal search progresses (since the contribution of each well to the measure chosen may depend on the respective well controls). Although not included in the results presented in this paper, the use of decline models to approximate metrics based on long-term simulation runs (e.g. lifelong reservoir recovery factor) may accelerate the optimization process further. We would like to note that geological uncertainty is not considered in this research (this type of uncertainty may not be as severe for mature fields as for greenfields). In any event, the retrospective optimization framework presented in Wang *et al.* (2012) could be readily used as an add-on to the approach introduced here in order to deal with uncertainty.

This paper is structured as follows. In the next section we state the optimization problem of interest and describe the methodology proposed to solve this problem. Thereafter, in the following section, we apply this methodology to a numerical model for a field in west Siberia, and later we report on actual experiments performed on that field to validate the solution obtained with the methodology. We end the paper with some conclusions and possible future research directions.

## Problem statement and methodology

In this research we address optimization problems in which optimization variables represent a discrete time series for different well control settings. Although a wide variety of well settings can be handled through these variables, here will be BHP values for the injection wells in the field of interest. The optimization cost function, , may refer to any performance metric. In this work, net present value (NPV) for a relatively short time frame (1 and 3 years) will be considered in all cases as the cost function:
(1)where , and are the oil sale price (this price is adjusted so that the separation cost and taxes are included), lifting cost and injection cost, respectively; , and are the rates for oil production, liquid production and water injection for each well *j* and *k*th time interval of duration ; and *d* is the annual discount rate where time *t _{k}* is expressed in years. We will also include bound constraints: that is, will be defined in , where and denote lower and upper bounds, respectively, and non-linear constraints. These latter constraints are mathematically formulated by means of a vector function , in which components have to be non-positive for feasible solutions. The specific non-linear constraint in our optimization problems refers to a measure of recovery factor (RF) computed for the same time frame as the NPV (this RF measure needs to be larger than that of a baseline plan). Since the methodology presented in this research can be applied to an arbitrary number of constraints, we will keep the notation (i.e. as vector function) throughout the paper. The optimization problem solved can be formulated as follows:
(2)where indicates any (locally) optimal solution, and the precise definition of the RF non-linear constraint is adjusted as required in order to translate feasibility into being non-positive.

### Derivative-free optimization

For the reasons mentioned above, we opt for derivative-free (i.e. simulator-non-intrusive) optimization methods. Most of the experiments performed are based on generalized pattern search (GPS: e.g. see Torczon 1997; Audet & Dennis 2002; Echeverría Ciaurri *et al*. 2011*b* for more details of that method and applications to oilfield operations), which is a local optimization technique (as noticed earlier, well control optimization problems have been observed to typically present multiple optima with similar cost-function values). Additionally, particle swarm optimization (PSO: e.g. see Eberhart & Kennedy 1995; Kennedy & Eberhart 1995; Onwunalu & Durlofsky 2010 for a description of the algorithm and its potential in field development), a global search procedure, has also been tested out in one example in this work. In both algorithms, the exploration of the search space is performed by means of a collection of points (stencil for GPS and swarm for PSO) that are modified in each iteration. We reiterate that the computational cost of many derivative-free methods scales with the number of optimization variables. Therefore, the optimization of mature fields with a relatively large number of wells is on some occasions a rather time-consuming solution. GPS and PSO tend to provide fairly robust and effective implementations (e.g. see Bellout *et al.* 2012; Humphries *et al.* 2014; Isebor *et al.* 2014 in addition to the references given for these two algorithms). Derivative-free alternatives to GPS and PSO can be found, for example, in Kramer *et al.* (2011).

### Reduction in the number of variables using streamline-derived information

Here we partially alleviate the problem of having a moderate or high number of variables by reducing the number of wells optimized in each control interval. This approach is based on constructing a linear model for reservoir fluid flow during a given period of time using streamline-derived information in the form of a number of volume ratios and efficiencies defined for all injection–production well pairs (for more details, see, for example, Batycky & Thiele 2006; Wen *et al.* 2014). This model (which essentially approximates the reservoir thorough a network of pipes) allows one to express production rates as a sum of contributions associated only with injection wells (and vice versa: i.e. injection rates can be written in terms of contributions corresponding only to production wells). Thus, a performance metric such as NPV can be decomposed into a number of components linked just to injection wells. We use these constitutive components to rank wells and, in turn, to select a subset of them that can be optimized. Metrics associated with non-linear constraints (e.g. RF in our case) can also be used for the ranking (otherwise, selections based on a metric solely related to the optimization cost function may favour infeasible well configurations). Note that a different subset of wells can be obtained for each of the control intervals considered, and that the subset selection is updated in every iteration of the optimization method applied. Although this streamline-based approximation is suited to wells controlled by means of rates, the precision obtained with it in brownfield management scenarios for other types of well controls can be often satisfactory (connections between injection and production wells tend to be stable for relatively long periods of time in mature fields, provided well controls are not modified drastically). This reduction approach might be attractive in real implementations because field operators prefer production plan updates where only a small number of controls are changed (in many cases, adjusting well controls implies travelling to remote locations of difficult access). It should be added that not all production scenarios are amenable to being modelled efficiently by streamline-based simulation but the management of brownfields, among others, is one of them (e.g. see Wen *et al.* 2014).

### Filter method for non-linear constraints

The non-linear constraints in the optimization problem in equation (2) are challenging to deal with because, as is the case with the cost function, simulation is needed to evaluate and, consequently, derivative information may not be available. The filter method (e.g. see Fletcher *et al.* 2006; Nocedal & Wright 2006) has recently been used in production optimization (e.g. see Echeverría Ciaurri *et al*. 2011*a*) to handle successfully this kind of constraints. A filter is essentially a new acceptance criterion for each optimization iteration that is based on concepts borrowed from bi-objective optimization. By means of the filter method, points are accepted/rejected not only based on the values of the cost function but also on a measure of constraint violation:
(3)where the maximum is determined for all components of , which are usually scaled so that they have comparable values. Thus, the optimal search with a filter is enriched using points that violate constraints and, thanks to that additional information, optimization techniques that include a filter for non-linear constraints may outperform other approaches that are based mainly on feasible points. The filter method is really an add-on to many optimization algorithms. The output of an optimization method that incorporates a filter is a collection of solutions that represent a trade-off between cost function and constraint violation. The user can select any of these solutions according to the flexibility available: that is, in some cases, constraint satisfaction can be relaxed if the improvement in cost function is acceptably large. In this work we combine the filter with GPS and with PSO (the reader may consult Isebor *et al.* 2014 for more details on the use of PSO with the filter method to deal with computationally expensive non-linear constraints for which derivatives are not available). Alternative approaches to handle general constraints in a derivative-free manner within the context of production optimization can be found, for example, in Echeverría Ciaurri *et al*. (2011*a*).

## Real field case study

### Pilot field description

The optimization methodology has been tested on Field ‘X’ of Gazprom Neft, located in west Siberia. The formation is a layered, uplifted reservoir with the presence of structural trapping. The reservoir, which is 12 × 4.3 km and 51 m thick, is mainly composed of sandy rocks interbedded with siltstones, dense clay and carbonate rocks. The formation is divided into 9–28 permeable layers. The average thickness for the pay sand is 25 m and the initial thickness for the net pay varies from 1.1 to 25.6 m. The average porosity, permeability and initial oil saturation are 18.8%, 23.6 mD and 61%, respectively.

Oil reserves were proved for Field ‘X’ through exploration drilling in 1989. The development of the field, mainly via waterflooding, started in 2006. Production started from the (vertical) exploration wells and additional horizontal wells were progressively drilled. The average production rate for the vertical and for the horizontal wells was 723 and 3189 barrels (bbl) per day, respectively. The peak oil rate took place from 2007 to 2012 and was approximately equal to 17 000 bbl/day. The lack of abrupt injection breakthroughs and gradual increase in water-cut suggests high fluid displacement efficiency. Field rates and the number of active wells as a function of time are represented in Figures 1 and 2. The optimization period considered in this paper started on 1 May 2016, when the field was in the fourth stage of development. At that moment, there were 19 production wells (16 of them were horizontal) and 15 injection wells (one of them was horizontal), and the field oil production rate was 6440 bbl/day. Injection and production rates for all the wells in the field at the beginning of the optimization time frame are represented in Figure 3. The waterflooding optimization is addressed through the adjustment of controls for all the 15 injection wells. This optimization strategy requires no capital expenditure and looks especially promising under conditions of relatively low oil price. The simulation model used in the optimization is based on a black-oil formulation and is discretized by means of a 70 × 145 × 114 grid (the number of active grid blocks is equal to 768 000). Each of these grid blocks has a length and width equal to 100 m, and an average height of 0.3 m. The history-matching quality of the model in most regions of the reservoir was deemed satisfactory for the optimization experiments performed. However, some wells were located in areas where this quality was questionable and, as will be discussed below, control changes affecting those wells in the optimized solutions were disregarded.

### Results of optimization

In the optimization we aimed to maximize the NPV as defined earlier in equation (1) over a production time frame equal to 1 year (first numerical test) and to 3 years (the rest of the numerical tests). The optimization constraint requires that the RF (or equivalently, the volume of oil produced) during that same time frame be larger than the RF associated with an existing baseline. The optimization variables describe sequences of BHP values that represent the controls for each injection well as a function of time. We have tested two derivative-free optimization algorithms: generalized pattern search (GPS) and particle swarm optimization (PSO). In both algorithms, the non-linear constraint for RF is handled by means of procedures based on the filter method (in the optimization problem there are also bound constraints, which can be dealt with in a relatively straightforward manner). The optimization starting point for GPS is always the baseline plan for both time frames of 1 and 3 years (this plan is not included in the initial swarm of PSO in order to better check how global the search performed by this method is). The GPS algorithm terminates when the size of the compass stencil (e.g. see Echeverría Ciaurri *et al*. 2011*b*) used to explore the optimization space is smaller than a predefined value (this smallest stencil size can be selected through considerations around the resolution deemed acceptable for the control variables), while PSO, which has a swarm size equal to 100, is run multiple times due to its random nature and is always stopped after 1600 cost-function evaluations (in our experiments this number seemed to be high enough for the algorithm to converge). This swarm size being somewhat larger than typical choices in optimal well location (e.g. see Onwunalu & Durlofsky 2010; Echeverría Ciaurri *et al*. 2011*b*; Wang *et al.* 2012), where, as mentioned above, cost functions may be considerably non-linear, may enhance the global nature of PSO in our application.

First, we need to decide on the size of the time interval where BHP is constant for each well. In principle, one would like to use as many intervals as possible and exploit that flexibility in the optimization to possibly obtain higher NPV. However, a large number of intervals will be translated into optimization problems with a considerable number of variables, which can be a time-consuming solution when derivative-free methods have to be used (we reiterate than in many cases, in practice, efficient gradient-based techniques cannot be applied because the user has no access to the simulation source code). Additionally, from the point of view of the field operator, frequent modifications of the well controls may not be desirable since these changes usually imply the execution of costly processes (e.g. travel to remote locations and technically complicated procedures). Here we have approached the optimization over a time frame of 1 year using well controls that are constant over intervals of either 3 months or the entire year. Five injection wells (selected by an expert) have been considered for optimization in this test (note that all wells in the field are taken into account in the simulations). Thus, the number of optimization variables for each of the two settings just described is equal to 20 and 5. The results obtained with the GPS algorithm are shown in Figure 4. In this figure we represent the evolution of the best feasible solution for the two optimization runs based on well controls that are constant over intervals of 3 months and 1 year. The difference between the optimized NPVs is relatively small (i.e. 0.96 and 1.13%) but the computational cost associated with the use of controls for 1 year is around 10 times smaller. Hence, in the rest of numerical tests (all of them for a total production time frame of 3 years) the piecewise constant representation of BHP controls for injection wells will be based on intervals of 1 year.

Next, we compare GPS with PSO when applied to the optimization of the same five wells as in the previous case but over a production time of 3 years (divided in three intervals of 1 year during which BHP is constant for each of the five wells). Therefore, this optimization problem has 15 variables. The results obtained are illustrated in Figure 5. Again, the plots represent the evolution of the best feasible solution in each optimization process. The PSO algorithm was run four times and the average of all these runs is shown in the plot (it should be noted that all the runs yielded the same optimized NPV; this result suggests that no more runs may be needed to assess performance for this configuration of PSO in this example). Both GPS and PSO reach values very close to the converged NPV (i.e. increase of 2.26%) after around 400 function evaluations. The PSO algorithm seems to take advantage of a large swarm size (equal to 100) to explore the space initially in a more thorough fashion and improve the cost function faster in spite of not including the baseline plan as a particle in the first swarm (the stochastic component of PSO may have contributed during the earliest stages of this optimization problem to deal with the non-linear constraint in more efficient manner). The systematic search performed by GPS (supported by mathematical convergence theory) is translated into more rapid convergence in the neighbourhood of the optimal solution.

It is important to notice that the optimized BHP controls for all solutions obtained with both optimization techniques are almost identical (in Fig. 6 we show the baseline, the solution found by GPS and one of the solutions determined with PSO; we can see that the two optimized solutions differ by a small margin in the third interval of injection well INJ-13; the differences between the solutions computed in PSO are negligible). This fact may indicate that this problem has a unique (global) optimum (as indicated above, well control optimization problems with a higher number of variables may present multiple solutions but with similar cost-function values). These results suggest as well that hybridization of global and local techniques (as proposed, for example, in Isebor *et al.* 2014) may somewhat accelerate convergence.

Finally, we address the control optimization for all the 15 injection wells in Field ‘X’ during three intervals of 1 year. The attendant optimization problem has 45 variables. As explained above, in this work we have introduced a procedure to accelerate the optimal search by means of streamline-derived information and the iterative selection of distinct subsets of wells for each of the (three, in this example) control intervals and for every iteration of the optimization algorithm. Different selections can be obtained depending on the performance metric that is used to rank the wells involved in the optimization. Here we consider two different performance metrics to choose five wells: NPV and the overall operational expenditures (OPEX) associated with all the wells in Field ‘X’. The selection of five wells instead of another number did not follow any particular guideline in this work. Research on criteria to choose that number could be a useful enhancement to the methodology. It is important to notice and emphasize that, in spite of selecting only five out of the 15 wells, we may still explore the original 45-dimensional space because the well selection could be modified in each iteration of the optimization algorithm (in other words, the overall method contemplates changes in the controls for all 15 wells). We additionally compare the two mentioned selection procedures with the optimization of the controls corresponding to five wells chosen by an expert (these five wells are the same for every control interval and each iteration of GPS).

The results of the four optimization runs are shown in Figure 7. In this figure, the runs that correspond to the selection of wells by means of an expert and through NPV and OPEX as performance metric are denoted by C1, C2 and C3, respectively. As expected, the explicit optimization of the 15 wells yields the highest optimized cost-function value (i.e. 3.91% increase in NPV with respect to the baseline) and demands almost 5000 cost-function evaluations, noticeably more than the other three approaches. Although the optimization problem was formulated for a production time frame of 3 years, it is important to verify that the optimized strategy does not underperform the baseline also in the long term. To this end, we have computed its NPV and RF during the period starting on 1 May 2016 and ending when cash flow in equation (1) is equal to zero. In this period, the well controls at the end of the first 3 years (i.e. on 30 April 2019) are kept untouched until production stops. Nothing guarantees that the overall strategy will be optimal for the entire reservoir life but, since both NPV and RF were the cost function and constraint of an optimization problem defined over a relatively long period of time, we can expect that long-term performance metrics may be acceptable. The optimized strategy (followed by unchanged controls) yields around 6 additional years of profitable production and the overall increase in NPV with respect to the lifelong performance of the baseline plan is 7.89%. It is worthwhile commenting that, even though in the optimization RF is introduced as a constraint, the RF obtained for the complete reservoir life is 1.14% higher than that of the baseline (note that this is an absolute and not a relative increase, i.e. 1.14% additional recoverable reserves are produced compared to the baseline). However, we stress the fact that the degree of uncertainty in long-term predictions is so high that conclusions based on this type of forecast need to be drawn carefully.

The procedure based on NPV (OPEX) for well selection yields an increase in NPV of 3.73% (3.35%) and is around 59% (66%) faster than the direct optimization of all wells. When the five wells are chosen by an expert, the optimized solution is computed in slightly less than 1000 cost-function evaluations but the corresponding increase in NPV is only of 2.26%. We can clearly see a trade-off between accuracy in the optimization and computing cost that we can adjust according to the precision needed in the solution and the resources available (e.g. time). It is worthwhile mentioning that, as pointed out earlier when describing the filter method, if constraint satisfaction can be relaxed, slightly infeasible solutions with cost-function values higher than the optimized one (which corresponds, in general, to a feasible solution) can also be outputted. In this particular optimization problem we have observed that the optimized solutions are located in the interior of the search domain (i.e. for each solution is strictly less than zero or, equivalently, the RF is higher than that of the baseline) and infeasible points only appear in the optimization process during its first stages (note that the initial guess in the optimization is at the boundary of the search space, i.e. is equal to zero). None of the infeasible points visited in the optimization runs were regarded as being superior to the feasible solutions found.

The solution with highest optimized NPV (i.e. the solution where BHP controls for all 15 injection wells were considered as optimization variables; in practice, one may not have computing time to optimize all 15 wells and may resort to selecting a subset of them using the procedure presented in this paper) was scrutinized before proposing the attendant experiment on the real field. Although this solution is based on adjustment of 14 out of the 15 injection wells, only changes in three of these wells were initially deemed acceptable from a practical point of view (the contribution to the increase in NPV from many of these wells is relatively small and, consequently, difficult to measure experimentally; in addition, some of the wells are located in regions of the reservoir where the quality of the history matching is questionable and, in consequence, modifications of their controls may not have the effects on the real field predicted through simulation). The difference in this study between the number of wells optimized and those eventually considered practically relevant is noticeable. Further research could be directed towards reducing the gap for these two numbers (e.g. by considering new formulations of the optimization problem that promote certain types of solutions).

The three wells chosen and the respective relative rate changes that correspond to the new BHP controls are shown in Figure 8. (It is important to note that the two solutions obtained using the reduction in the number of wells via streamline-derived information also included relatively major changes for these three wells.) The most drastic recommendation is related to injection well INJ-2: the optimized solution requires shutting the well in (i.e. relative change of −100%). The other two recommendations refer to an increase in BHP for wells INJ-1 and INJ-3. The high injection rate of INJ-2 in the baseline plan might hinder effective fluid displacement from INJ-1 to the five nearby production wells indicated in Figure 8. Therefore, the shutting-in of INJ-2 combined with a modest injection rate increase of INJ-1 may distribute better fluid flow and, in turn, produce the reserves in that area in a more effective manner. Production wells around INJ-3 may also benefit from a higher rate at this injection well.

Taking all the selected changes of the well controls in total, we expect that a reduction in the overall injection rate of 23% for INJ-1, INJ-2 and INJ-3 will not be translated into a decrease in oil production (due to a decline in reservoir pressure) because other remaining reserves will be produced; additionally, the savings in injection and lifting costs will contribute further to an NPV increase. We reiterate that the well settings considered in the optimization are BHP values (the actual settings that can be modified for injection wells in the real field are choke sizes, which are easier to relate to BHP values than to rates).

### Results of field experiment

The experts responsible for direct operation and control of Field ‘X’ examined the optimized solution with the highest NPV. From the three wells identified above relevant in that solution, only one of them, INJ-2, was selected for actual implementation in the real field. In addition, this well, for which shutting-in was recommended in the optimization, would be still active although with −50% relative injection change. There are two main reasons why the complete optimized solution cannot be tested out in practice. First, there are important infrastructural constraints that were not taken into account in the numerical model used in the optimization (e.g. technological limitations for the injection redistribution in the current surface network). The use of a comprehensive model that integrates the reservoir and facilities may allow one to deal with these constraints at the expense of a considerably higher computational cost than the already time-consuming reservoir flow models. This more complex optimization process may be a possible direction for future research. In the second place, the high degree of uncertainty typically present in real-field applications is very often translated into the questioning of new plans that depart noticeably from the status quo. Hence, rather than modifying significantly the settings for a given well, it may be a reasonable strategy to introduce relatively small changes in the current controls that are representative of those associated with the drastically new settings (e.g. −50% relative injection change instead of shutting the well in): that is, perform some kind of sensitivity analysis in the real field in order to validate the solution proposed. Consistent with these limitations, the choke for INJ-2 was replaced by a smaller one, which provides, approximately, a relative 50% injection reduction. Simulation results indicated that this recommendation yields around 40% of the increase in NPV given by the complete (optimized) solution, together with satisfaction of the short-term RF constraint. The results of the corresponding field study are presented next.

In this analysis, Field ‘X’ will be divided in two regions: the area affected by the optimized policy, which will be referred to as the region of the experiment, and comprises injection wells INJ-1, INJ-2 and INJ-3 and production wells PRO-1, PRO-2, PRO-3, PRO-4 and PRO-5, as shown in Figure 8; and the rest of the reservoir. The selection of the wells that belong to the region of the experiment was made, as in the reduction stage described earlier, using streamline-derived information (i.e. the wells chosen were those with a relatively strong connection with the injection wells INJ-1, INJ-2 and INJ-3). We represent in Figure 9 the monthly operational rates for the region of the experiment before and after the optimized plan is applied. As can be seen in that figure, the optimized plan yields in the region of the experiment a significant drop in water injected together with a noticeable water-cut reduction and slight decrease in oil produced (note that this decrease in oil production, typical for a brownfield, is clearly more pronounced for the rest of the reservoir, as will be illustrated below). In order to better quantify this improvement in production performance, these results have to be compared with a reference. Since during a given time frame only one production plan can be implemented, an ideal comparison can never be made in practice. A reasonable approach would be based on using a baseline determined with the trends for the existing plan just before the optimized strategy started. However, and this is especially common in the management of mature reservoirs, field development implies continuous monitoring and corresponding adjustments of the production plan anywhere in the field. This is translated into frequent engineering activities, which may happen simultaneously and in the same region as the optimized strategy. In our case, during the time frame of the test, the operator performed a number of treatments that affected the entire field (including the region of the experiment). The effects of these activities should be somehow taken into account in the analysis and that can be a rather challenging enterprise when real data are considered. To this end we proceed as follows. First, we compare production performance before and after 1 May 2016 for the two regions into which the field is divided. This comparison, when analyzed outside of the region of the experiment, may allow us to roughly quantify the effects associated with the treatments performed by the operator and, with adequate scaling, estimate their contribution and deduct it within the region of the experiment. Although this procedure implicitly assumes that the regions have similar flow behaviour and does not include cross-effects between them, it could provide useful information as far as the impact of our methodology on performance metrics of interest in the management of Field ‘X’.

Field measured data are analyzed next according to the three different components in equation (1) for the computation of NPV: water injection, and oil and liquid production. In Figure 10, we illustrate the reduction in water injection rate for the wells in the region of the experiment and also for the rest of the reservoir (note that a different vertical axis is used for each region). The average water injection rate in the region of the experiment decreased from 7090 to 5430 bbl/day. We estimate that the reduction in the total amount of water injected during the 8 months (245 days) of the experiment was, on average, 407 000 bbl. Separately, we quantify the mean decrease in water injected for the rest of the reservoir as 172 000 bbl. If this reduction is scaled using the average baseline injection rates for both regions of the reservoir (i.e. 7090/21 400) and then is subtracted from the estimation above for the region of the experiment, we can conclude that in this region the average savings in water injected during the 8 months of the experiment could be approximately 350 000 bbl. It is apparent that these savings have a noticeable effect in the operational expenditure component of the NPV.

Figure 11 represents the reduction in fluid production rate for both regions of Field ‘X’ (the respective water-cut for each region is also plotted in the figure). If we repeat the process followed above, we now find that the average reduction in the volume of liquid produced for the region of the experiment and for the rest of the field are 19 000 and 241 000 bbl, respectively. Since the liquid production rate in the latter region is much larger than the corresponding rate in the former region (even when we consider the scaling factor, which is taken equal to 2.6 in this case), we conclude that our methodology may actually yield a relatively small increase in liquid production volume, and we estimate this increase to be around 73 700 bbl. This result, together with the fact that the water-cut in the region of the experiment decreases in somewhat clearer form than in the rest of the field, might indicate a higher relative oil production for our approach with respect to the reference. Oil production rates for the two regions of Field ‘X’ are shown in Figure 12. In order to compute differences between oil volumes we resort to a linear decline model, which is often a reasonable approximation for short time frames. We thus find that the average increase in the volume of oil produced with respect to the baseline for the region of the experiment and for the remaining part of the field are 25 300 and 24 700 bbl, respectively. If we note that the oil production rate in the second region is about 1.9 times the rate in the first region, we can conclude that our methodology may yield around 12 300 bbl of oil more than the reference.

The overall results from the experiment indicate that the application of the approach presented in this work seems to be beneficial as far as NPV is concerned. The increase in revenue and the reduction in water injection expenses outweighs in our analysis (for most logical economic scenarios) the slightly higher lifting costs. The optimized strategy appears to produce more oil than the baseline even with less water injected and at the expense of a relatively small amount of additional volume of liquid production. We would like to add that the methodology presented could be enhanced when implemented in a closed-loop optimization fashion: that is, measurements could be leveraged to update the existing flow models and with that provide recommendations based on improved production forecasts.

## Conclusions

In this paper we have presented a methodology that can be used as a decision-support system in the optimal waterflood management of brownfields. Optimal waterflooding may be an appealing solution (especially in periods of relatively low values of oil price) because it usually requires low capital investments, it generates recommendations that are in most cases reversible actions (unlike, for example, fracturing or sidetracking operations) and it allows relatively fast verification of its viability (in contrast to, for example, many treatments in enhanced oil recovery).

The approach described relies on simulation-based, derivative-free optimization techniques that are non-intrusive regarding the particular simulator(s) used (i.e. in principle, any simulator can be incorporated into the system). Gradient-based methods, when derivatives are computed efficiently, often provide rapid solutions in production optimization scenarios (the computational cost of derivative-free optimization scales, in general, with the number of control variables). However, fast computation of accurate gradients is typically a simulator-intrusive procedure that cannot be applied in many situations (e.g. when some commercial simulators are used).

Simulation-non-intrusive optimization can readily handle a wide range of problem formulations, including most types of well controls (assuming the underlying simulation incorporates these controls as input). Performance metrics such as net present value (NPV) and recovery factor (RF) can be considered as the optimization cost function and constraints (the approach can accommodate flexibility in the satisfaction of the constraints by outputting slightly infeasible solutions but with satisfactory cost-function values, as demonstrated, for example, in Echeverría Ciaurri *et al*. 2011*a*).

The number of well controls in many production optimization problems may be too high in order for derivative-free optimization to provide solutions efficiently. Rate controls can be optimized in a relatively fast manner by means of procedures that leverage streamline-based simulation. Nevertheless, these procedures are not applicable to other kinds of well controls such as, for example, BHP values that are easier to translate into real-field implementations than rates. In this research we have introduced a procedure that by using streamline-derived information selects the most important wells in the optimization according to a previously chosen metric and, in turn, accelerates the entire computing process. Although the approach has been illustrated with the optimization of BHP values, any type of well controls can be considered. It should be noted that streamline-based simulation is a modelling technique especially suited to certain production scenarios – brownfield management being one of them.

The overall approach was applied to a real brownfield located in west Siberia that currently has 15 injection wells and 19 production wells. The optimization problem relied on a history-matched model of the field and aimed at maximizing the NPV over a production time frame of 3 years subject to a constraint for the RF defined for the same period.

Multiple solutions that represent a tradeoff between precision and computational cost were computed. If only five wells are optimized and these wells are manually chosen by an expert, the increase in NPV with respect to the baseline is 2.26%. When the procedure based on streamline-derived information was used to also pick five wells, the NPV obtained was 3.73% higher than the baseline (nonetheless, the approach was computationally more expensive than resorting to an expert to reduce the number of wells). The optimized strategy with the highest NPV (which also corresponded to the most time-consuming optimization run and entailed considering controls for all 15 injection wells) yielded an increase of 3.91% with respect to the baseline (the attendant RF was also higher than that of the reference but, since this metric was introduced only as a constraint, the observed increase cannot be ascribed to the optimization process). This solution is equally satisfactory during the complete reservoir life (i.e. while cash flow is not negative), even though the optimization problem was only formulated over 3 years. The increase in lifelong NPV and RF (regarding the baseline) associated with that strategy is 7.89 and 1.14%, respectively.

The controls for one injection well were modified over about 8 months according to the optimized plan. However, other field operations than this change in the controls for one well took place concurrently (this could be a frequent situation in real scenarios because operators may not be able to cancel all previously planned actions). In order to deal with this situation, we developed a procedure to approximately quantify the impact of our strategy, in the presence of others, on the management of the field. This procedure could be useful for the operator to better assess the performance of multiple recovery mechanisms that are applied simultaneously which respective effects may be difficult to tell apart.

During the 8 months of the experiment, we estimated that the savings in water injected according to our solution were around 350 000 bbl. Our estimation for the increase in liquid and oil production was 73 700 and 12 300 bbl, respectively. Consistent with the numerical results obtained in the optimization, the extra volume of oil produced together with the reduction in injection costs clearly outweigh in most economic scenarios the lifting costs related to the additional volume of liquid produced (firmer conclusions could be drawn for longer field experiments). In any event, we can expect these promising results to extrapolate to other brownfields scenarios.

This research can be extended in several directions. In the first place, besides testing alternative optimization techniques (e.g. the use of approximate derivatives via ensemble-based and stochastic-gradient methods), the optimal search process can be accelerated if surrogates are considered at different stages of the workflow. The selection of the number of wells in the reduced set could be supported by a number of guidelines related, for example, to the computational resources available. Formulations of the optimal search problem may include measures of the plausibility of the solutions found in order to make these solutions more readily implementable on the real field. The underlying simulation-based model in the optimization may integrate reservoir and surface facilities, and in that manner capture technical constraints that are in many cases disregarded. This type of integrated model might facilitate the optimization of well controls that are more realistic than bottom-hole pressure (BHP) and rate (e.g. choke size and pump frequency). However, the computational cost of the optimization problems thus obtained may in many cases be prohibitive, and strategies to speed up computations significantly would then be essential. Uncertainty considerations may be included in the process to endow solutions with additional robustness and flexibility (different future actions can be recommended based on upcoming information and on the risk attitude of the decision-maker). Finally, cognitive computing (e.g. see Kelly & Hamm 2013) may be leveraged in the whole system to strengthen the interaction between human and computational components.

## Acknowledgements

The authors would like to thank Arseniy P. Gruzdev from IBM for his implementation of the particle swarm optimization.

## Funding

This work was supported by a funded joint research project between Gazprom Neft and IBM to develop robust technology for Gazprom Neft's mature assets in west Siberia. The theoretical outcome of this project formed a basis for a software prototype that was applied to a real oilfield.

- © 2018 The Author(s). Published by The Geological Society of London for GSL and EAGE. All rights reserved