Open Access
Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles
Volume 76, 2021
Article Number 7
Number of page(s) 16
Published online 25 January 2021

© D. Santos Oliveira et al., published by IFP Energies nouvelles, 2021

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1 Introduction

Development and management of petroleum fields generally require estimation of many parameters (e.g., permeability, porosity, oil, water cut, etc.), so the best production scheme for each field can be applied. In reservoir simulation, physics-based models are used to obtain these parameters, but in most cases, they are difficult to be applied due to its high computational cost. An alternative to deal with this problem is to use proxy models. These models have been widely used for decision making process in recent years, and present as a powerful technique to analyze reservoir dynamics with a low computational cost and a quick response.

During the development of oil fields, one of the most used methods to maintain reservoir pressure and improve oil recovery is the waterflooding method. For this method, there are several injectors wells operating simultaneously in the field that inject water into the reservoir to increase its pressure, and the amount of water injected by a certain injector may influence the flow rate of one or more producer wells. At the surface facilities, liquid production rates are separated in, e.g., water, gas and oil rates, and can be processed to obtain valuable and useful information about the relationship between producer and injector wells. In any active field, production and injection flow rates are generally the most abundant data available, and can be used to optimize oil recovery when employed to forecast injection flow rates, determine when and where drill new wells, improve the production design and estimate inter-well connectivity.

There are several different techniques, statics and dynamics, to estimate the communication between wells. Some Static techniques are interference test (Du et al., 2002; Maizeret, 2013) and tracer test (Zhang et al., 2001). As example of Dynamic techniques, we have streamlines analyses (Datta-Gupta and King, 2007), flow diagnostics (Moyner et al., 2014) and signal processing techniques (Wang et al., 2020) and the Capacitance Resistance Model (CRM) by Yousef et al. (2005). Some traditional techniques, e.g. the interference test, require in some level reservoir characterization or fluid distribution description to run, what may limit their applicability. In this context, the CRM shows as an alternative to estimate inter-well connectivity, once only production data is needed to run the model. It also can be used to qualitatively indicate the presence of channel, sealing faults, and allocate wells (Weber, 2009). Usually its parameters are used to forecast field production. CRM is a data-driven model that can carry out field history matching and heterogeneities detections only using available field production history data. Sayarpour (2008) details the CRM development and shows an analytical solution based on the superposition in time and space for three different reservoir control volumes: the entire field (CRMT); the drainage volume of a single producer (CRMP), and the volume between an injector-producer pair (CRMPI). For example, Lins et al. (2017) applied the CRM to compute the connectivity between layers in a compartmentalized reservoir, and Mirzayev and Jensen (2019) estimate the connectivity in a low permeability reservoir using pseudo-wells.

In this study, we apply the CRMP model, widely used in recent studies (Lins et al., 2017; Tueros et al., 2018), that describes the relationship between a pair of injector and producer with only one parameter called connectivity. This model has three parameters: connectivity, that represents the rate fraction of a given injector that influences the liquid rate of a specified producer; time constant, that quantifies the degree of fluid storage between wells; and productivity index, that indicates how the pressure changes in the pore drainage volume around a producer affect its rate.

The CRMP parameters can be obtained using different strategies: Nonlinear Regression (Albertoni and Lake, 2003; Gentil, 2005; Kaviani and Valkó, 2010; Liang et al., 2007; Sayarpour et al., 2007), Nonlinear Least Squares method (Holanda et al., 2015; Lins et al., 2017; Moreno, 2013; Tueros et al., 2018; Weber et al., 2009) and Neural Networks (Liu et al., 2020; Panda and Chopra, 1998). Here, we used the nonlinear least squared method. Therefore, with CRMP parameters and injection flow rates, the total liquid rates of producer wells can be predicted, and, if well type and location are kept fixed, the injected water volume can be optimized (Holanda et al., 2015; Hong et al., 2017; Jafroodi and Zhang, 2011; Sayarpour et al., 2009; Weber et al., 2009). Theorical arguments about the physical meanings and the strategies to obtain CRM parameters can be found in (Gentil, 2005; Noetinger, 2016). For example, Shahamat et al. (2016) used the capacitance and resistance concept and developed a mathematical model for compartmentalized reservoirs and detailed analytic solutions.

Over the past decade, CRMP strategy has been used in decision making process for waterflooding operations, but as every method, it has its advantages and disadvantages. As main advantage it presents the use of injection and production data as input to estimate fluid dynamics in the reservoir, without any prior knowledge of its physical and fluid properties. However, as a disadvantage CRMP cannot separate oil and water rates from output results.

The main objective of this work is to maximize the profitability of the field, and for this, it is necessary to know oil production, and water injection and production of wells. However, the CRMP cannot separate oil and water produced volumes. Therefore, fractional flow models can be used for this task. The most used models in the literature are Koval (Cao, 2014) and Gentil (Gentil, 2005).

The Koval model has been used for cases where the wells have low water production, immature fields (Cao, 2014), but, Hong et al. (2017) proposed a modified model for mature fields. Gentil (2005) proposed a model for mature fields with Wcut above 50%. In this work, we propose to combine the fractional flow models, Kogen, to address the weaknesses of the Koval and Gentil models using Wcut as a transition parameter between both models. The parameters of Koval and Gentil models can be obtained by NLS.

The SQP algorithm is used to determine Kogen parameters, and well control trajectories in the waterflooding problem (Alim, 2013; Dehdari and Oliver, 2012; Pinto et al., 2020; Tueros et al., 2020). SQP is a gradient-based algorithm, so derivatives are required. Computing derivatives by finite differences in large scale problem becomes computationally prohibitive. A viable alternative is to use the ensemble-based method proposed by Chen et al. (2009). This method has been widely used in reservoir engineering by (Chen and Oliver, 2012; Chen and Reynolds, 2016; Do and Reynolds, 2013; Tueros et al., 2015; Zhang et al., 2018). Tueros et al. (2018) proposed refinement techniques on the ensemble-based method. The main advantages of the ensemble-based method are reduced number of simulations in the optimization process and readily coupling with reservoir simulators.

As mentioned earlier, we use the combined Kogen model to obtain oil and water cumulative productions and the ensemble-based method for calculating derivatives in the optimization process. Both strategies are applied to two reservoir models existing in the literature. We observed that the proposed strategy reduces the computational cost and the number of iterations and function evaluations during the optimization process.

2 Optimization problem formulation

In this section we present the waterflooding problem, optimization algorithms used to solve the problem, approximate computation of gradient and smoothing strategies to avoid wild changes of controls in the optimized trajectories.

2.1 Waterflooding problem

Here, the waterflooding problem with capacity operations may be formulated as an optimization problem with linear constraints where design variables are the controls of injector wells:


In equation (1), the NPV is the objective function, d is the discount rate and τt is the reference time period at control cycle t. qt = the vector of rates at control cycle t and qw,t = injection rate of the well w at control cycle t. I = set of indices to the NI injectors, the life-cycle period subdivided into Nt control cycles and Qinj,max = the maximum allowed total injection rate of the field. Superscripts l and u represent the lower and upper bounds of design variables.

The cash flow at the control cycle t represents the oil revenue minus the cost of the water injection and production, given by:


where P = set of indices to the NP producers, ∆τt = time interval of the tth control cycle; and are average oil and water rates at the wth production well at tth control cycle, respectively;  = average water rate at the wth injector well at tth control cycle; ro = oil price, cwp and cwi = cost of water producer and injected, respectively. Indices op, wp and wi denote the oil, water production and water injection phases, respectively. On the other hands, in the problem, equation (1), we can also impose nonlinear constraints as maximum field liquid production.

In this study, the control variables are the injection rates in each control cycle which are normalized by the respective maximum injection rate of the well, qw,max:


With normalization of controls in equation (3), we consider the following alternative formulation:


where x = column vector that contains all injection well controls;  = normalization factor of the function f(x) normally taken as NPV(x0); x0 = initial point; if xw,t = 0 means the tth control has zero flow and the well is closed at cycle t and the value 1 is the maximum allowed rate for injector wells.

2.2 Optimization algorithms

In this study, we use two optimization algorithms to solve our problems. Among the gradient-based optimization algorithms for solving large-scale problem, one of the most efficient both in number of iterations and function evaluation required in the optimization process is the SQP algorithm (Dehdari et al., 2012; Liu et al., 2018). A detailed description of the SQP can be found in Biegler (2010). In addition, the SQP solves a subproblem, in each iteration to approximate the original problem. Another method that used is the Nonlinear Least Squares that is detailed below.

2.2.1 Nonlinear Least Squared method

In this work, we use the NLS method to solve the uncoupled problem of Section 4.1. In this method, the objective function is defined as:


where Ti(x) = is the difference between the numerical reservoir simulator and the approximate model responses. For example, the real value may be the liquid rate of a given producer well obtained by the simulator and the approximate value is obtained using the CRMP model. By minimizing this function, we find values for the parameters that best correlate the approximate model and the real response.

The gradient vector and Hessian matrix of the objective function in equation (6) are given by:


Assuming that Ti(x) tends to zero near the solution, then the first term of matrix of second derivatives (Eq. (6)) will tend to zero, and the Hessian matrix of the second derivatives f(x) can be approximated using only the first order derivatives of .

The NLS is generally based on the Gauss-Newton approximation, which displays a quadratic convergence ratio near the solution, for those cases where the Hessian approximation is accurate, that is, when the residuals tend to zero in solution.

The advantage of using this method is that due to its mathematical structure it allows obtaining the Hessian matrix using only first-order derivatives. Moreover, by exploring this problem structure, the second-order convergence characteristic of Newton’s algorithm can be obtained using only first-order derivatives.

2.3 Computation of the gradient by ensemble-based method

The SQP algorithm required the gradient vector in the optimization process. Here, we approximate the gradient using the ensemble-based method by Chen and Oliver (2012). The gradient vector computed by ensemble-based method is a nonlocal estimate considering the subspace spanned by the adopted control realizations that has a strong dependence on ensemble size (Chen and Oliver, 2012). The advantages are the reduced number of simulations in the optimization process, and the use of the simulator as a black-box and being therefore independent of the specific adopted software. It is necessary to generate different control realizations that must be simulated to compute the NPV. In this work, we use the steps proposed by Chen et al. (2009), Dehdari and Oliver (2012), Su and Oliver (2010), Tueros et al. (2018) to generate the realizations of time-correlated controls:

  1. Use Spherical covariance function (Do and Reynolds, 2013; Tueros et al., 2015) to generate the prior covariance matrix, C, for control variables, avoiding abrupt changes with time.

  2. Compute the Cholesky decomposition of the covariance matrix, LLT = C.

  3. Inform the initial values of the controls to be perturbed, (this is only informed on the first iteration, for the following iterations use the optimal point of the previous iteration), where Nx = total number of control variables.

  4. Adopt a standard deviation to produce perturbations of the control variables δ. The magnitude of the perturbations set by δ needs to be balanced not to be too small or too large, which is problem dependent (Su and Oliver, 2010).

  5. Create a vector of independent random variables with zero mean and unit variance, Z.

  6. A realization of the vector of control variables, , is given by:


The second term in equation (7) generates smooth control realizations for gradient computation. This does not guarantee smoothness of controls during the optimization process (Tueros et al., 2020). In this study, the realizations of control variables for ensemble-based satisfy all constraints (e.g., wells and capacity allowed).

2.4 Computation of the gradient using sensitivity matrix

The following steps are used to calculate the sensitivity matrix and to obtain the ensemble-based approximate gradient (Tueros et al., 2020):

  1. Generate a matrix containing the set of all feasible control realizations using equation (7), where Nr = ensemble size. In this work, we used ensemble size with 10, 20 and 30 realizations.

  2. Compute the ensemble-based covariance matrix (Chen and Oliver, 2012).

  3. Calculate the ensemble-based cross-covariance matrix between controls and well NPV’s.

  4. Compute the sensitivity matrix as the product of the pseudo-inverse of the covariance matrix by the cross-covariance matrix.

  5. The approximate gradient is the sum of the columns of the sensitivity matrix.

The accuracy of the sensitivity matrix is strongly dependent on ensemble size. However, using large sizes increases computational cost during the optimization process.

2.5 Smoothing factor

In this study, we use the smoothing technique proposed by Tueros et al. (2020) devised to avoid abrupt jumps in control trajectories during the optimization process. The strategy is based on adding a new term to the objective function of equation (4). The optimization problem with the smoothing term can be expressed as:


where Nv = total number of controls by wells; β = smoothing factor. We obtain the original optimization problem, if this factor is zero. In this study, we used β = 10−3 for all examples as recommended by Tueros et al. (2020). The gradient of the smoothing term can be easily computed analytically.

2.6 Illustrative example

To develop an understanding of the strategies proposed, we present a synthetic model based on the Brush Canyon Outcrop (BCO; Oliveira, 2006), with an internal sealing fault (Fig. 1). The cells that model fault have permeability, porosity, and transmissibility set to zero. The model has six layers, with 55 × 53 active cells. The horizontal permeability varies from 157 to 2592 mD (Fig. 1b). The vertical permeability is 30% of the horizontal permeability, porosity varies between 16% and 28%, and viscosity is 0.11 cP. The reservoir has twelve wells, seven producers and five injectors, and wells are completed vertically (Fig. 1a). Producer wells operate with constant BHP (11 000 kPa). Injector rate has an operational range between 0 and 2500 m/day. The lifetime of the reservoir is 24 years, 4 for production history and 20 for the optimization process. The wells opening during production history is shown in Table 1.

thumbnail Fig. 1

(a) Wells locations for BCO-Fault model; (b) permeability field.

Table 1

Production control history for 4 years (BCO-Fault model).

3 Capacitance and resistance model

This work uses the Capacitance and Resistance Model based on Producer wells (CRMP) proposed by Sayarpour et al. (2007), which is derived from the material balance equation (oil and water) and does not depend on well locations to quantify the communication (connectivity) between a pair of injector and producer wells. The CRMP has three parameters that collectively produce a good representation of the reservoir flow dynamics. Its great advantage is not to require any geological or fluid reservoir characterization to construct the model, using only the available data from the field production history.

3.1 CRMP

In waterflooding operations, there are multiple producers and injectors operating simultaneously and, consequently, a given producer can have its production rate influenced by more than one injector at a certain period. CRMP considers injector rates and BHP of producers as inputs signals that perturb the reservoir dynamic, generating the total liquid production rates as outputs (Fig. 2).

thumbnail Fig. 2

Reservoir processing system.

The derivative continuity equation for a producer well j is expressed as shown by equation (9) and, it is developed from the material balance equation. Its parameters and derivation can be found in Sayarpour et al. (2007):


where τj = time constant for a producer j; qj(t) = total fluid production of producer j at time t; λij = connectivity for a pair injector/producer well; Ii(t) = rate of injector i at time t; Jj = productivity index of producer j; Pwf,j = bottom hole pressure of producer j; ct = total compressibility; Vp = drainage pore volume of a producer j (connected porous volume).

Integrating equation (9) over a discrete time period, ∆t, we obtain the total liquid rate of a producer j, at reservoir conditions, at a certain time k, which can be written as (Holanda et al., 2015; Lins et al., 2017; Sayarpour et al., 2007; Weber, 2009):


where, the first term in the right-hand side corresponds to the initial response of the production rate associated with primary depletion, the second is associated with the injection contribution and the last term with the influences of BHP changes in the producers. Figure 3 shows a schematic representation of the effect of injector wells on the liquid rate of a given producer in waterflooding operations. We can also observe that for a certain field, there will be multiple producers and injectors operating simultaneously and, subsequently, more than one injector well can influence the production of a given producer well at a certain period.

thumbnail Fig. 3

Influences of injector wells on producer well (Oliveira et al., 2020).

Gentil (2005) and Yousef et al. (2005) discuss the physical meaning of the connectivity and time constant parameters, respectively, and their insights were used for a better understanding of the CRMP model assumptions.

To determine the CRMP modal parameters λij, τj, Jj, i = 1…NI; j = 1 … Np, we solve the following constrained nonlinear least squares problem (Weber, 2009):


where, = observed liquid rate of the producer j at time step k and = liquid production rate calculated by the CRMP for the producer j at time step k. The observed data is obtained from the actual field production history, but here we use a commercial reservoir simulator CMG-IMEX (2018) to generate them. It is important to highlight that the number of unknown parameters depends if BHP is considering constant or variable. In this work, we consider constant minimum BHP in all producer wells.

The parameters of the CRMP problem are obtained by solving two optimization problems: an uncoupled problem, where only limit constraints are imposed, and another that considers all constraints, coupled. The uncoupled problem response serves as a warm start point for the other problem and has the advantage of reducing the number of iterations during the optimization process in the coupled problem.

The parameters are computed using one realization of controls, knowing that the approximation depends mainly on the available data (Albertoni and Lake, 2003). The following steps are followed for the parameter computation (Lins et al., 2017; Tueros et al., 2018; Weber, 2009):

  1. Set random BHP values for producers and injection rates on the simulator. These values are set to change every 30 days during the simulation period (Holanda et al., 2015; Tueros et al., 2018).

  2. Use production data after all producer’s present go through water breakthrough.

  3. Using NLS method to solve the uncoupled problem of equation (12) for each producer well, j = 1 … Np, to warm start the main coupled problem.


  1. For each injector well, i = 1 … NI, impose feasibility with respect to main problem unit connectivities sum constraint (Lins et al., 2017; Tueros et al., 2018).

  2. Using the computed starting point, solve coupled problem of equation (11) applying an SQP algorithm.

If a pair of injector/producer wells presents a high connectivity value, it evidences the presence of a channel or a fracture between them. However, a low connectivity value indicates the existence of a low permeability zone or barrier in the interwell region.

Table 2 presents the parameter results for BCO model. Wells P-1 and P-4 present high connectivity values with respect to I-1 and I-4, which are the injector wells at the same side of the fault. Although, their connectivities related to the remaining injectors are smaller with exception for the pair P-1/I-2 with a high interwell permeability zone that increases its connectivity value.

Table 2

BCO-Fault connectivities.

Regarding the time constant, P-4 presents the highest value which is explained by its location in the reservoir. This well is surrounded by low permeability regions and its closest injector is I-4. By contrast, P-2 shows the lowest value for the time constant indicating that the pressure changes in its surroundings are felt in shorter amount of time, as a result of the proximity to I-2 and I-3.

For the other producer wells, the same analysis was made to validate the results. Also, a history match and rate prediction were conducted to check the CRMP parameters quality and reliability. The time window for history matching was chosen such that water breakthrough had already occurred at all producers to reduce effects of nonlinearities in the CRM parameters (Holanda et al., 2015). For this case, the chosen time window was between 0 and 2200 days. The curves for history match and rate prediction for the wells P-1 and P-5 are shown in Figure 4. From the results, one can notice that CRMP model captures the reservoir dynamics well using only production data.

thumbnail Fig. 4

Matching and prediction of total liquid rate for well P-1 and P-4, BCO-Fault model.

4 Fractional flow models

In the previous section, the applied CRMP model computes total liquid rate of producer wells as a single phase. However, in most of the field operations and waterflooding optimization problems is important to separate oil and water rates.

From the liquid rate, qj,k, it can be obtained water production rate, , and oil rate, , of producer j at time k as:


where Wcutj,k = water cut of producer j at time k. In equation (13) it is necessary to calculate the Wcut of producer wells.

Therefore, we proposed a combined fractional flow model (Koval + Gentil) to adjust Wcut curves of producer wells. Koval, Gentil and Kogen fractional flow models are presented in the next sections.

4.1 Koval model

Koval model was introduced by Sayarpour (2008) and it is derivative from the Koval Theory (1963), which is used to capture the effects of water fingering generated by reservoir heterogeneity and unfavorable mobility ratio between displacing and displaced fluid. Cao et al. (2015) combine the CRMP and Koval models to separate oil and water production. This strategy shows a good forecast when the reservoir is immature.

From the Buckley–Leverett equation, modifying the viscosities ratios and after some rewriting, we obtain the following equation to determine the Wcut (Cao et al., 2015; Hong et al., 2017):


where kvalj = Koval factor of producer j; tD,j = the dimensionless time for a producer j; CIj,k = the cumulative water injected contribution from all injectors to producer j during a certain period ∆t and, Vp,j = the drainage pore volume of producer j. We assume that the water injection rate of an injector is constant over ∆t and equal to Ij,k, respectively. Using the connectivity parameters, λij, from CRMP, the CIj,k is defined as:


The dimensionless time, tD,j, may be interpreted as the volume of cumulative water injected in terms of pore volumes, implying that we must identify the contributions from all injectors to the production of a certain producer j in a reservoir, which can be determined using the CRMP model as shown in equation (15) (Cao, 2014).

The modal parameters kvalj and Vp,j, can be found by solving the following constrained nonlinear least squared problem:


where, Vp, field = the total pore volume of the field. In this work, the kvall = 1. It is important to highlight that if Vp, field is not available, the constraint can be neglected. In this work the problem is solved with limit constraint only unlike procedures proposed by Cao et al. (2015) and Hong et al. (2017). Therefore, parameters can be obtained independently for each well.

Figure 5 compares Wcut values obtained by the simulator and Koval model for wells P-1, P-4 and P-5. We can see a good forecast in the early years and the remaining time shows loss of quality. Therefore, using this Wcut provides a strong inaccuracy in the calculation of oil rate. Similar results were also observed in the works of Cao (2014) and Hong et al. (2017).

thumbnail Fig. 5

Water cut obtained using koval model.

4.2 Gentil model

In Section 4.1 it was discussed a fractional model that performs better when Wcut is small and the reservoir did not achieve significant water breakthrough. However, for large multi-wells reservoir, the time that breakthrough is achieved for each producer well will vary significantly, with some wells producing high amount of water whereas others do not. This diverse system highlights the necessity of a data driven model applicable when the Water-Oil Ratio (WOR) is greater than 0.5, on other words, when the producer wells have high Wcut value.

The Wcutj,k can be obtained from water-oil ratio as:


where WORj,k = Water-Oil Ratio of producer j at time k. Gentil (2005) proposed an empirical power law to fit the historical data by assuming a nonlinear relationship between water-oil ratio and the cumulative water injected, CIj,k (Eq. (15)) and express as Liang et al. (2007):


where, aj and bj are the fractional flow model parameters for each producer j. This expression can be transformed into a linear form given by Gentil (2005) and express as:


Equation (19), has no constraints that couple the Gentil model parameters for all producer wells, so their fit can be determined independently to history values of the water/oil rate natural log, LWORj,k.

Equation (19) is valid only for later time data (e.g. the water cut must be greater than 0.5 (Gentil, 2005). This happens because applying the natural logarithm to equation (18), the linear relationship between water oil ratio and “CI” is only observed when Wcut is large. This model is usually fitted to a smaller time window than the CRMP model.

The decision variables in the Gentil model can be obtained by solving the following nonlinear least squared problem:


where = observed natural log of the water-oil rate for producer j at time step k. = the log of the calculate water/oil rate for producer j at time step k, equation (19), and koil = the initial value of the time window used to calculate the parameters. A similar procedure was used by Weber (2009).

Figure 6 shows the results obtained for wells P-1, P-4 and P-5. It can be noticed that even though Gentil approach is applied for mature waterflood reservoirs, if the Wcut change with time is not so abrupt, the model is capable to capture the Wcut behavior for values lower than 0.7. However, most results confirm that analysis data window should be in late horizon time, where the Wcut is greater than 0.7 in which case numerical errors are consistently smaller.

thumbnail Fig. 6

Water cut obtained using Gentil model.

4.3 Kogen model

In Figure 7a we compare the water cut obtained by the simulator with Koval and Gentil models. It can be observed that Koval shows good fit for initial times and Gentil for later times. To overcome the limitations of Koval and Gentil, we propose the combined Kogen model based on the previously mentioned models. Kogen fits well to immature and mature fields, using only available production history data and CRMP parameters. Kogen main parameter is based on the difference in water cuts of the Koval and Gentil models at the moment of transition between the two (Fig. 7b), hereafter referred to as ts.

thumbnail Fig. 7

(a) Koval and Gentil model approaches; (b) Kogen model scheme.

The combined Kogen model can be formulated as a constrained nonlinear curve fitting, expressed as:


In equation (21) the objective function, MisFitj represents a measure of the difference between calculated and observed Wcut or NPV. The nonlinear constraint is the difference in Wcut values of the Koval, Wcutj,sk, and Gentil, Wcutj,sg, models at the time of transition, ts. In Equation (21), α is the jump allowed between Koval and Gentil models. In this work we used α = 20%. The other parameters are associated with the Koval and Gentil models. The problem is solved using the SQP algorithm and its gradient is obtained by finite differences.

As commented previously MisFitj function can be defined by either of the following two expressions:



where, = net present value of the producer j observed by simulator in the step time k; = net present value of the producer j calculated in the step k; = water cut of the producer j observed by simulator in the step time k. = water cut of the producer j calculated in the step k.

The performance of both expressions for misfit function are compared in Figures 8 and 9. It can be observed the Wcut-based function performs slightly better than the NPV-based.

thumbnail Fig. 8

Water cut obtained using simulator and Kogen model, function NPV, for producers P-1 (a), P-4 (b), and P-5 (c).

thumbnail Fig. 9

Water cut obtained using simulator and Kogen model, function Wcut, for producers P-1 (a), P-4 (b), and P-5 (c).

Results show that the proposed strategy is able to better approximate water cut with the simulator result, when compared to the individual Koval and Gentil models. Even though Kogen has more variables than the other models, it is fast and requires only minimal information.

Furthermore, in order to be able to evaluate waterflood performance we develop a combined fractional flow model that together with CRMP technique generates an efficient proxy model, which produces promising starting points for robust optimal management problems.

5 Examples

Two synthetic reservoirs, BCO and Brugge models, are used to assess the proposed strategies, the proxy model CRMP + Fractional flow models, Koval, Gentil and Kogen, for function evaluation during the optimization process. The optimization algorithm of choice is the SQP available in MATLAB Optimization Toolbox (2014). The gradient is computed by an ensemble-based method. Commercial simulator IMEX is used to generate high fidelity field production data.

The results statistics of the objective function are presented in box-plots. The limits of the box-plot are 25% and 75% percentile, the line (red) in the box-plot is the median, the extreme lines are the maximum and minimum values, and the crosses are the outliers. The economic parameters used to calculate the NPV for the all examples are oil revenue, 80 $/bbl, water injected and produced price, 5$/bbl, and annual discount rate of 10% per year.

5.1 Example 1. BCO model

The BCO model described in Section 2.6 is used in three different studies using an ensemble-based method for the optimization process. The first two use the same random starting controls. First, we analyze the most appropriate ensemble size for computing the gradient vector by comparing optimization results, using simulator to evaluate production response. Second, we make a comparative study of the optimization results using the fractional flow models against the previous best results. Finally, we use the best solution of the fractional flow models as a warm start point for a new high-fidelity optimization process in order to assess the quality of proxy-based results.

The objective is to maximize NPV for a period of 24 years, four for history and 20 for optimization. Well controls are changed every six months, i.e., there are 40 control cycles. The total number of variables is 200. The rate of injector wells can vary from 0 to 2500 m3/day, while the maximum allowed field injection rate is 7000 m3/day.

5.1.1 BCO by ensemble-based method

In this exercise, we analyze the size of the ensemble for calculating the gradient vector. The initial NPV for this case is 3.991 × 109. The ensemble-based method is of a stochastic nature, so to have efficient comparisons, 10 optimization process runs were performed.

Figure 10 compares optimization results obtained with different ensemble sizes. The first plot (left) was obtained with an ensemble of Nr = 10 and the others with 20 and 30 respectively. The median obtained with Nr = 30 was 5.45 × 109 and shows a small variability in the results. The best response was 5.507 × 109, requiring 1400 calls from the simulator during the optimization process. The gain compared to the medians of the plots Nr = 10 and Nr = 20 was 5% and 2%, respectively. From the results presented in Figure 10, we observe that when we increase the size of the ensemble, from 20 to 30, the gain is small. For these reasons, we adopt Nr = 30 for the other examples in this work.

thumbnail Fig. 10

Results using different ensemble sizes, Nr = 10, Nr = 20 and Nr = 30.

In Figure 11 we show the optimal trajectories of the wells for the best result obtained in Figure 10. The smoothing strategy helped to reduce abrupt changes in the trajectory of the controls. The injectors I-4 and I-5 are kept at maximum rate. Injectors I-2 and I-3 remain closed for the first few years.

thumbnail Fig.11

Optimal control for injector wells using Nr = 30, best solution.

5.1.2 Comparison of different fractional flow models

With the injection rates and CRMP parameters we predict the liquid rates for the production wells. To calculate the NPV, it is necessary to determine the oil and water rates. In this example, we analyze two scenarios: first, we compare the performance of the three fractional flow models, Koval, Gentil and Kogen, to separate the oil and water rates during the optimization process. Finally, the best solution is compared with the simulator production results. The initial NPV obtained with the Koval, Gentil and Kogen models are 3.351 × 109; 3.873 × 109 and 3.941×109, respectively (3.991 × 109, using the simulator).

Figure 12a shows the optimization results using the fractional models (Koval, Gentil and Kogen; from left to right) in three box plots. The median for Kogen is 4.975 × 109, showing a gain of 20.2% and 2.65% with respect to the median values of Koval and Gentil. The best response for each of the models in Figure 12a is 4.188 × 109; 4.878 × 109 and 5.015 × 109.

thumbnail Fig. 12

Comparison of results: Koval, Gentil and Kogen (a) and iteration history (b).

Figure 12b shows the iteration histories for the best result of each fractional model. The best NPV by Kogen was obtained in 8 iterations, while Koval and Gentil in 7 and 13 iterations, respectively. Finally, Kogen proved to be the most efficient.

Another case studied was to compute optimum NPV using Kogen and replacing the ensemble-based method with the Finite Difference (FD) method for calculating the gradient. The final NPV was 5.388 × 109 in 30 iterations. The FD solution improved the NPV but at larger number of Kogen calls.

In Table 3 we compare the response of the fractional flow models with the results obtained directly by the simulator. In all cases we observed that there is a reduction in NPV.

Table 3

Comparison of NPV results (× 109 $): Koval, Gentil and Kogen.

5.1.3 Ensemble-based optimization with a warm start

In this case we use the best solution for fractional flow models as a starting point for high-fidelity optimization (ensemble-based method using simulator response). The comparison of the results is shown in the four box plots of Figure 13. For example, the curve for Koval represents the results after high-fidelity optimization with the best result obtained by Koval model. We find that the results are consistent with those obtained in Section 5.1.1. The main advantage of the strategy is the reduction in the number of simulator calls during the optimization process. For example, the best responses using Koval, Gentil and Kogen as warm starts were obtained in 12, 10 and 9 iterations with 540, 420 and 420 simulator calls, respectively, which are significantly less than the 1400 simulations used to the best response of Figure 10. In the following examples, we use only the Kogen model.

thumbnail Fig. 13

Results of high-fidelity optimization using the best results obtained by Koval, Gentil and Kogen models.

thumbnail Fig. 14

Well locations of Brugge model.

5.2 Example 2. Brugge model

The example is based on the Brugge model with single completion per well. A detail description of the model is documented in Peters et al. (2009). We used the realization #28 (Fig. 14). The objective is to maximize NPV for a period of 20 years, using the preceding 10 years for production history. All wells operate from the very first day. An explanation about the connectivity parameters computation of the Brugge model can be found in Oliveira et al. (2020).

All 400 design control variables are the rates of the 10 injectors well using six months control steps. The rate of injector wells can vary from 0 to 4000 bbl/day, while the maximum field total injection capacity is 30 000 bbl/day. The producer wells operate with constant BHP of 1000 psi.

Below we present two studies: first, in order to establish the appropriate ensemble size, we compare high fidelity optimization results varying the number of ensemble members for computing the gradient vector. Second, we use the ensemble-based method with Kogen model to maximize NPV and the best response is used as a warm start point for high-fidelity optimization. The initial well controls were randomly selected.

5.2.1 Brugge model by ensemble-based optimization

In this exercise, the initial NPV is 1.12 × 109. The results of high-fidelity optimization are shown in two box plots of Figure 15. The leftmost plot is obtained for ensemble of Nr = 20 and the other with Nr = 30. The median value with Nr = 30 is 1.54 × 109 and shows small variability in responses. The best NPV is 1.555 × 109, where 1330 simulator calls were needed in 30 iterations. The gain obtained with Nr = 30 is 3% with respect to the median value for Nr = 20. Henceforth we adopt Nr = 30. Figure 16 shows optimal control trajectory for the best result. We observed that optimal controls do not show abrupt jumps.

thumbnail Fig. 15

Results of high-fidelity optimization varying ensemble size.

thumbnail Fig. 16

Optimal control trajectories for all injector wells. (a) wells 1–4; (b) wells 5–7 and (c) wells 8–10.

5.2.2 Ensemble-based optimization with Kogen model

Initial NPV using is 1.01 × 109 which is close to 1.12 × 109 using the simulator responses. The leftmost box plot of Figure 17 refers to results obtained using the same random initial controls of the previous section. The other plot uses the controls of the best solution as starting point. The median value of the first plot is 1.30 × 109 and maximum response is 1.348 × 109 requiring 15 iterations. Using simulator results we obtain an NPV of 1.328 × 109. High-fidelity simulator results are used for the second plot. The median NPV is 1.537 × 109, while the best result 1.543 × 109 was obtained in 10 iterations with 480 simulator calls.

thumbnail Fig. 17

Result using Kogen and Ensemble-based.

6 Discussion and conclusion

Fractional flow models present in the literature are used specifically for immature or mature fields. In this work, we present a general combined model, Koval + Gentil, called Kogen, to separate oil and water from total liquid production of a producer well, using only the production data and CRMP parameters. The performance of the proposed model was compared with the fractional flow models, Koval and Gentil, for optimization processes. The main advantage of Kogen model is its capacity to forecast the reservoir production during the whole concession period.

Two strategies were presented to determine Kogen model parameters: a misfit function based on water cut and another based on NPV. The results show that the first presented a slightly better adjustment. The adopted maximum water cut difference at the switch time between Koval and Gentil models was 20%.

The combination of CRMP and Kogen models, proxy CRMP-KOGEN, works as a tool for forecasting oil and water production, which is essential for optimizing the field production. The results were used as a warm start point for high-fidelity optimization, which uses the simulator, in order to decrease the number of functions evaluations.

Two examples from the literature, BCO and Brugge, were used to evaluate the advantages of the proposed strategy. The algorithm selected to solve the optimization problem was the SQP, and the ensemble-based method was used to calculate the gradient vector.

A parametric study was developed to define the most appropriate size of ensemble for the high-fidelity optimization process. It was observed that with Nr = 30 the gains obtained were 2% and 3% when compared with results for Nr = 20, in the BCO and Brugge models, respectively. Therefore, it is recommended to use 30 realizations per ensemble for future examples. With this value, water injection optimizations were performed using proxy models with CRMP coupled to the three fractional flow models comparing final NPV values. Proxy-Kogen model showed gains of 20.2% and 2.65% with respect to the proxy models with Koval and Gentil fractional flow models, respectively. When the optimal controls of proxy models were run in the simulator it was observed a reduction of about 3.22% in NPV values.

When the proxy-Kogen model result was used as a warm start for the high-fidelity optimization process, we observed gains of 15% in the final response. This strategy also allowed to reduce the number of simulations by approximately 65%, during the high-fidelity optimization process.


The authors acknowledge the financial support for this research by PRH-47 Human Resources Program, PETROBRAS, CAPES and Energi Simulation.


  • Albertoni A., Lake L. (2003) Inferring interwell connectivity only from well-rate fluctuations in waterfloods, SPE Reserv. Eval. Eng 6, 1, 6–16. [CrossRef] [Google Scholar]
  • Alim M. (2013) Constraint handling in life-cycle optimization using ensemble gradients, MSc Thesis Report, Delft University of Technology, The Netherlands. [Google Scholar]
  • Biegler L.T. (2010) Nonlinear programming: concepts, algorithms, and applications to chemical processes, SIAM – Society for Industrial and Applied Mathematics, Philadelphia. [CrossRef] [Google Scholar]
  • Cao F. (2014) Development of a two-phase flow coupled capacitance resistance model, PhD Dissertation, The University of Texas at Austin, USA. [Google Scholar]
  • Cao F., Luo H., Lake L.W. (2015) Oil-rate forecast by inferring fractional-flow models from field data with Koval method combined with the capacitance/resistance model, SPE Reserv. Eval. Eng. 18, 4, 534–553. [CrossRef] [Google Scholar]
  • Chen Y., Oliver D.S., Zhang D. (2009) Efficient ensemble-based closed-loop production optimization, SPE J 14, 4, 634–645. [CrossRef] [Google Scholar]
  • Chen Y., Oliver D.S. (2012) Localization of ensemble-based control-setting updates for production optimization, Society of Petroleum Engineers, [Google Scholar]
  • Chen B., Reynolds A.C. (2016) Ensemble-based optimization of the water alternating-gas-injection process, SPE J. 21, 3, 786–798. [CrossRef] [Google Scholar]
  • CMG (2018) IMEX: User’s Guide, Computer Modeling Group LTD, Calgary-Canada. [Google Scholar]
  • Datta-Gupta A., King M.J. (2007) Streamline Simulation: Theory and Practice, Textbook Series #11, Society of Petroleum Engineers, Richardson, TX, ISBN 978-1-55563-111-6. [Google Scholar]
  • Dehdari V., Oliver D.S. (2012) Sequential quadratic programming for solving constrained production optimization – case study from Brugge field, SPE J. 17, 874–884. [CrossRef] [Google Scholar]
  • Dehdari V., Oliver D.S., Deutsch C. (2012) V.: Comparison of optimization algorithms for reservoir management with constraints – A case study, J. Pet. Sci. Eng. 100, 41–49. [CrossRef] [Google Scholar]
  • Do S.T., Reynolds A.C. (2013) Theoretical connections between optimization algorithms based on an approximate gradient, Comput. Geosci. 17, 959. [CrossRef] [Google Scholar]
  • Du W.B., Zhang J.G., Hao S.Z. (2002) The application of pulse well tests for the well group of Forerunner Test of tertiary oil recovery in Yumen oilfield, Well Testing 11, 4, 32–33. [Google Scholar]
  • Gentil P.H. (2005) The use of multilinear regression models in patterned waterfloods: Physical meaning of the regression coefficient, Master’s thesis, The University of Texas at Austin, USA. [Google Scholar]
  • Holanda R.W., Gildin E., Jensen J.L. (2015) Improved waterflood analysis using the capacitance-resistance model within a control systems framework, Society of Petroleum Engineers. [Google Scholar]
  • Hong A.J., Bratvold R.B., Nævdal G. (2017) Robust production optimization with capacitance-resistance model as precursor, Comput. Geosci. 21, 1423–1442. [CrossRef] [Google Scholar]
  • Jafroodi N., Zhang D. (2011) New method for reservoir characterization and optimization using CRM–EnOpt approach, J. Pet. Sci. Eng. 77, 2, 155–171. [CrossRef] [Google Scholar]
  • Kaviani D., Valkó PP (2010) Inferring interwell connectivity using multiwell productivity index (MPI), J Pet Sci Eng. 73, 1, 48–58. [CrossRef] [Google Scholar]
  • Koval E.J. (1963) A method for predicting the performance of unstable miscible displacement in heterogeneous media, SPE J 3, 2. [Google Scholar]
  • Liang X., Weber D., Edgar T.F., Lake L.W., Sayarpour M., Al-Yousef A. (2007) Optimization of Oil Production Based on A Capacitance Model of Production and Injection Rates, Society of Petroleum Engineers, Richardson, TX, USA. [Google Scholar]
  • Lins H.K.F., Horowitz B., Tueros J.A.R. (2017) Numerical experience using capacitance resistance multilayered models, CILAMCE 2017, in: Ibero-Latin American Congress in Computational Methods in Engineering. Florianopolis, Brazil. (in Portuguese). [Google Scholar]
  • Liu Z., Forouzanfar F., Zhao Y. (2018) Comparison of SQP and AL algorithms for deterministic constrained production optimization of hydrocarbon reservoirs, J. Pet. Sci. Eng. 171, 542–557. [CrossRef] [Google Scholar]
  • Liu W., Liu W.D., Gu J. (2020) A Machine Learning Method to infer inter-well connectivity using bottom-hole pressure data, J. Energy Resour. Technol. 142, 103007. [CrossRef] [Google Scholar]
  • Maizeret P.D. (2013) Best practice to design and interpret interference tests based on features of the line-sources solution: Theory and Applications, Society of Petroleum Engineers. [Google Scholar]
  • MATLAB’s Optimization Toolbox User’s Guide, R2014. [Google Scholar]
  • Mirzayev M., Jensen J.L. (2019) Interwell connectivity analysis in a low-permeability formation using a modified Capacitance Model with application to the East Pembina Field, Cardium Formation, Alberta, Oil Gas Sci. Technol. - Rev. IFP Energies Nouvelles 74, 26. [CrossRef] [Google Scholar]
  • Moreno G.A. (2013) Multilayer capacitance–resistance model with dynamic connectivities, J. Pet. Sci. Eng. 109, 298–307. [CrossRef] [Google Scholar]
  • Moyner O., Krogstad S., Lie K.A. (2014) The Application of Flow Diagnostics for Reservoir Management, Internal Report SINTEF 1, 1–30. [Google Scholar]
  • Noetinger B. (2016) About the determination of quasi steady state storativity and connectivity matrix of wells in 3D heterogeneous formations, Math. Geosci. 48, 6, 641–662. [Google Scholar]
  • Oliveira D.F.B. (2006) Production optimization techniques for petroleum reservoirs: Derivate free approaches to dynamic rate allocation for injection and production, Master Thesis, Civil Engineering Department, UFPE, Recife, Brazil (in Portuguese). [Google Scholar]
  • Oliveira D.S., Horowitz B., Tueros J.A.R. (2020) Kogen-Combined Koval/Gentil Fractional Flow Model, in: Proceeding of the 17th European Conference on the Mathematical Oil Recovery, 14–17 September, 2020. [Google Scholar]
  • Panda M., Chopra A. (1998) An integrated approach to estimate well interactions, in: SPE India oil and gas conference and exhibition, Society of Petroleum Engineers. [Google Scholar]
  • Peters E., Arts R.J., Brouwer G.K., Geel C.R. (2009) Results of the Brugge benchmark study for flooding optimisation and history matching, Paper SPE 119094 presented at the SPE Reservoir Simulation Symposium, The Woodlands, Texas, USA, 2–4 February. [Google Scholar]
  • Pinto J.W.O., Tueros J.A.R., Horowitz B., Afonso da Silva S.M.B., Willmersdorf R.B., Oliveira D.F.B. (2020) Gradient-free strategies to robust well control optimization, Comput. Geosci. 24, 1959–1978. [CrossRef] [Google Scholar]
  • Sayarpour M., Zuluaga E., Kabir C.S., Lake L.W. (2007) The use of capacitance–resistive models for rapid estimation of waterflood performance, SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers, [Google Scholar]
  • Sayarpour M. (2008) Development and application of capacitance-resistive models to water/CO2 floods, Ph.D. Dissertation, The University of Texas at Austin, Austin, Texas. [Google Scholar]
  • Sayarpour M., Zuluaga E., Kabir C.S., Lake L.W. (2009) The use of capacitance-resistive models for rapid estimation of waterflood performance and optimization, J. Petrol. Sci. Eng. 69, 3–4, 227–238. [CrossRef] [Google Scholar]
  • Shahamat M.S., Hamdi H., Mattar L., Aguilera R. (2016) A novel method for performance analysis of compartmentalized reservoirs, Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles 71, 38. [CrossRef] [Google Scholar]
  • Su H.-J., Oliver D.S. (2010) Smart Well Production Optimization Using An Ensemble-Based Method, SPE Res. Eval. Eng. 13, 6, 884–892. [CrossRef] [Google Scholar]
  • Tueros J.A.R., Horowitz B., Willmersdorf R. (2015) Numerical experience with an ensemble-based method for constrained waterfloonding optimization, in: Cilamce 2015, Ibero-Latin American Congress in Computational Methods in Engineering, Rio Janeiro, Brazil. [Google Scholar]
  • Tueros J.A.R., Horowitz B., Willmersdorf R.B., Oliveira D.F.B. (2018) Non-distance-based localization techniques for ensemble-based waterflooding optimization, J. Pet. Sci. Eng. 170, 440–452. [CrossRef] [Google Scholar]
  • Tueros J.A.R., Horowitz B., Willmersdorf R.B., Oliveira D.F.B. (2020) Refined ensemble-based waterflooding optimization subject to field-wide constraints, Comput. Geosci. 24, 871–887. [CrossRef] [Google Scholar]
  • Wang Y., Kabir C.S., Reza Z. (2020) Deciphering well connectivity with diagnostic signal processing techniques, J. Pet. Sci. Eng. 185, 106610. [CrossRef] [Google Scholar]
  • Weber D. (2009) The use of capacitance-resistance models to optimize injection allocation and well location in water floods, Ph.D Dissertation, The University of Texas at Austin, Austin, Texas. [Google Scholar]
  • Weber D., Edgar T.F., Lake L.W., Lasdon L.S., Kawas S., Sayarpour M. (2009) Improvements in capacitance–resistive modeling and optimization of large-scale reservoirs, in: SPE Western Regional Meeting, Society of Petroleum Engineers. [Google Scholar]
  • Yousef A.A., Gentil P.H., Jensen J.L., Lake L.W. (2005) A capacitance model to infer interwell connectivity from production and injection rate fluctuations, in: SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers. [Google Scholar]
  • Zhang Y., Jiang R.Z., Zheng X.Q. (2001) Analysis technology of inter-well tracer, J. China Univ. Pet. (Ed. Nat. Sci.) 25, 2, 76–78. [Google Scholar]
  • Zhang Y.T., Stordal A.S., Lorentzen R.J., Chang Y. (2018) A Novel Ensemble-Based Conjugate Gradient Method for Reservoir Management, Society of Petroleum Engineers. [Google Scholar]

All Tables

Table 1

Production control history for 4 years (BCO-Fault model).

Table 2

BCO-Fault connectivities.

Table 3

Comparison of NPV results (× 109 $): Koval, Gentil and Kogen.

All Figures

thumbnail Fig. 1

(a) Wells locations for BCO-Fault model; (b) permeability field.

In the text
thumbnail Fig. 2

Reservoir processing system.

In the text
thumbnail Fig. 3

Influences of injector wells on producer well (Oliveira et al., 2020).

In the text
thumbnail Fig. 4

Matching and prediction of total liquid rate for well P-1 and P-4, BCO-Fault model.

In the text
thumbnail Fig. 5

Water cut obtained using koval model.

In the text
thumbnail Fig. 6

Water cut obtained using Gentil model.

In the text
thumbnail Fig. 7

(a) Koval and Gentil model approaches; (b) Kogen model scheme.

In the text
thumbnail Fig. 8

Water cut obtained using simulator and Kogen model, function NPV, for producers P-1 (a), P-4 (b), and P-5 (c).

In the text
thumbnail Fig. 9

Water cut obtained using simulator and Kogen model, function Wcut, for producers P-1 (a), P-4 (b), and P-5 (c).

In the text
thumbnail Fig. 10

Results using different ensemble sizes, Nr = 10, Nr = 20 and Nr = 30.

In the text
thumbnail Fig.11

Optimal control for injector wells using Nr = 30, best solution.

In the text
thumbnail Fig. 12

Comparison of results: Koval, Gentil and Kogen (a) and iteration history (b).

In the text
thumbnail Fig. 13

Results of high-fidelity optimization using the best results obtained by Koval, Gentil and Kogen models.

In the text
thumbnail Fig. 14

Well locations of Brugge model.

In the text
thumbnail Fig. 15

Results of high-fidelity optimization varying ensemble size.

In the text
thumbnail Fig. 16

Optimal control trajectories for all injector wells. (a) wells 1–4; (b) wells 5–7 and (c) wells 8–10.

In the text
thumbnail Fig. 17

Result using Kogen and Ensemble-based.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.