Determination of lower and upper bounds of predicted production from history-matched models

We present a method to determine lower and upper bounds to the predicted production or any other economic objective from history-matched reservoir models. The method consists of two steps: 1) performing a traditional computer-assisted history match of a reservoir model with the objective to minimize the mismatch between predicted and observed production data through adjusting the grid block permeability values of the model. 2) performing two optimization exercises to minimize and maximize an economic objective over the remaining field life, for a fixed production strategy, by manipulating the same grid block permeabilities, however without significantly changing the mismatch obtained under step 1. This is accomplished through a hierarchical optimization procedure that limits the solution space of a secondary optimization problem to the (approximate) null space of the primary optimization problem. We applied this procedure to two different reservoir models. We performed a history match based on synthetic data, starting from a uniform prior and using a gradient-based minimization procedure. After history matching, minimization and maximization of the net present value (NPV), using a fixed control strategy, were executed as secondary optimization problems by changing the model parameters while staying close to the null space


Introduction
It is well known that assimilation of production data into reservoir models is an ill-posed problem; see, e.g.[15,20,24].This is mainly because generally the number of uncertain model parameters largely supersedes the number of measurements.Moreover, the measurements are strongly correlated because they originate from a relatively small number of sources: the wells.As a result, they contain less information about the true value of the model parameters than could be expected based solely on the number of data points.A relevant question in view of the purpose of large-scale, physics-based reservoir models is how much the long-term predictions can vary because of the ill-posedness of the assimilation problem.In other words, what may be the economic consequences of the lack of information about the reservoir in the measurements?
In most practical circumstances, this question is addressed by constructing and history-matching low-and high-case models, besides the nominal model.Alternatively, a set of model realizations can be used in a data-assimilation algorithm to obtain an entire collection of predictions, as is the case with ensemble Kalman filter (EnKF) methods; see, e.g.[1,7].However, in either way, the resulting history-matched models are heavily influenced by the prior information that went into the data-assimilation process.Hence, properly answering the question stated above requires either some (heuristic) method to translate static geological properties to flow behaviour or economic performance, or requires many forward simulation runs to obtain a proper low-or high-case prior model.These methods may be either unreliable or impractical to provide a good measure of the economic consequences of the lack of knowledge about the true field.
In this paper, a method is introduced to search for lower and upper bounds on predicted production (or any other economic objective) over the remaining life of a field, using a history-matched model.The method consists of two steps: (1) performing a traditional computer-assisted history match of a reservoir model with the objective to minimize the mismatch between predicted and observed production data through adjusting the permeability values of the model and (2) performing two optimization exercises to minimize and maximize an economic objective over the remaining field life, for a fixed production strategy, by manipulating the same grid block permeabilities, however without significantly changing the mismatch obtained under step 1.To achieve this, we make use of the fact that history matching through adjusting grid block parameters is an ill-posed problem such that many combinations of parameter values may result in (nearly) identical mismatch values.

Problem definition
The problem of determining a history-matched model that provides either a lower or an upper bound on the predicted economic performance over the life of a reservoir is essentially a multi-objective optimization problem.For a general overview of multi-objective optimization, see, e.g.[14].The first objective is to find a certain realization of model parameters that minimizes the error between the measured and simulated production data, which can be expressed through a quantitative objective function V , e.g.mean square difference.The second objective relates to finding a set of parameter values that-for a certain future production strategy-minimizes or maximizes some economic cost function J , e.g.net present value (NPV).However, the multiple objectives are not of the same importance; priority lies with obtaining a good history match, while determining a lower or upper bound on predicted economic performance serves as a secondary objective.To that end, the multiobjective optimization problem may be cast into a hierarchical optimization problem, as presented in [10] and more recently specifically for oil production optimization in [4,8,23].In this structure, optimization of a (secondary) economic cost function J is constrained by the requirement that the (primary) quantitative history-matching cost function V must remain close to its minimal value V min .This requires solving the following two (hierarchical) optimization problems, and where ū is the fixed control vector (input vector), x is the state vector (typically grid block pressures and saturations), g is a vector-valued function that represents the system equations, x 0 is a vector of the initial conditions of the reservoir, the subscript k indicates discrete time and K is the total number of time steps.The vector of inequality constraints c relates to the capacity limitations of the wells.The term ε is small value compared to V min .In order to solve the secondary optimization problem, given in Eq. 4 to Eq. 7, first, a (single) optimal solution to the primary optimization problem Eq. 1 to Eq. 3 is required to determine V min .The optimal solution to the primary problem θ * 1 can serve as a feasible initial guess for the secondary problem.Note that the second optimization problem is also solved in terms of θ , while the values of the control ū remain unchanged.The search space of the secondary problem is now constrained by the null space of the primary objective function at a value of V min , through inequality constraint Eq. 7. In other words, the redundant degrees of freedom (DOF) of the primary problem are the DOF of the secondary problem.The motivation for using the constraint Eq. 7 is actually twofold.If ε is arbitrarily small (or even equal to 0), the parameter space that remains is actually the null space within the parameter space, which can be substantial because of the generally ill-posed nature of the inverse history-matching problem.Any changes of the model parameters within that null space will have no effect on the value of the used quantitative history-match quality indicator, i.e. the objective function V .For ε > 0, the corresponding parameter space that satisfies Eq. 7 can be given the interpretation of a parameter uncertainty set, with a clear statistical interpretation, in the case of Gaussian noise disturbances on the data.The statistical uncertainty set then results from a hypothesis test based on the so-called likelihood ratio test, and is characterized by level sets of the likelihood function V ( θ).See, e.g.[18] for the case of nonlinear models and [5] for linear models.This implies that under appropriate noise conditions, we can, for every value of ε > 0, connect a probability level to the parameter uncertainty set defined by Eq. 7 and thus account for the variability of the history-matched parameters in the subsequent secondary economic optimization problem.

Methodology
In [23], the primary optimization problem was attacked using a gradient-based search algorithm.(Note that in that study the optimization variables were the inputs u, while here they are the model parameters θ .)The gradients were obtained using a system of adjoint equations which was solved backwards in time once, regardless of the number of optimization parameters.(See [11] for an overview of adjoint-based optimization in porous media flow and [12] for the specific implementation used in this study.)Subsequently, the secondary optimization problem was also attacked using a gradientbased search algorithm.However, the secondary problem was executed with the addition of projecting the search direction onto a second-order approximation of the null space with respect to the optimality constraint defined in Eq. (7).The second-order approximation was explicitly determined through a forward difference scheme using first-order information obtained with the adjoint.Unfortunately, using this approach, the number of forward and backward simulations is proportional to the number of optimization parameters.Hence, for the assimilation of production data, this method is in most cases computationally infeasible.
In [23], also an alternative method was introduced to solve the hierarchical optimization problem without explicitly calculating the null space with respect to Eq. (7).It uses an 'on-off' type weighted objective function with weighting functions 1 and 2 : where 1 and 2 are 'switching' functions of V and J that take on values of 1 and 0 ('on' and 'off') or vice versa, Here, ε is the threshold value as defined in inequality constraint Eq. (7).The gradient of W with respect to the model parameters θ for iteration n + 1 is then simply, Solving the secondary optimization problem sequentially, using W as defined in Eq. ( 8), gives improving directions for either V or J .With each iteration, the value of J either increases while the value of V decreases or the other way around, as the solution moves to and from the feasible region with respect to inequality constraint Eq. ( 7).If there exist redundant DOF with respect to the primary problem, improvement of J is possible while satisfying Eq. ( 7) and convergence of the hierarchical optimization will occur in a 'zigzag' fashion, as schematically represented in Fig. 1.
To improve convergence speed, as presented above and in [23], a small adaptation to the switching algorithm can be made.By projecting the gradient of secondary objective function J onto the first-order approximation of the null space of the primary objective function V , the resulting update of θ will keep V closer to V min .Mathematically this becomes where we use the convention that the derivative of a scalar with respect to a vector is a row vector.P | | is a matrix that projects ∂J /∂θ on ∂V /∂θ , and (∂J /∂θ)P ⊥ is the orthogonally complementary projection which ensures that the step towards the secondary objective function is taken in a direction (near-)parallel to the 'ridge' in the primary objective function.Inserting Eq. ( 11) in Eq. ( 8) gives an alternative switching search direction d for solving the hierarchical optimization problem The switching algorithm using the projected gradient d was used in the following example to illustrate the performance of the method.
Fig. 1 Schematic representation of the iterative process of solving a hierarchical optimization problem using a weighted objective function, as described by Eq. 8.The process converges towards a final solution in a 'zigzag' fashion, moving into and out of the feasible region bounded by the optimal solutions of the primary objective function.(After [23])

Primary and secondary objective functions
Two different reservoir models are used in this paper with the goal of determining lower and upper bounds on the expected economic performance over the remaining life of the field by changing the permeability field, while the model stays compliant with historic data over the history-matching period.Consequently, the primary objective function, V ( θ), is defined as the data mismatch between observations and simulated data: where θ is a vector of unknown model parameters, d is a vector of data (measurements), h is a vector-valued function that relates the model parameters to the model outputs (i.e. the simulated data) and P v is a covariance matrix of data errors.
The secondary objective function,J , is of an economic type, generally the NPV, where y wp,j is the water production rate of well j ; y o,j is the oil production rate of well j ; y wi,j is the water injection rate of well j ; r wi , r wp and r o are water injection costs, water production costs and oil revenues respectively; t k is the time interval of time step k in days; b is the discount rate for a reference time τ t ; and N inj and N prod are the number of injection and production wells.

Egg model example
In this first example, initially presented by [22], we consider a three-dimensional oil reservoir model, introduced for a different purpose in [21].The reservoir model consists of 18,553 active grid blocks, as depicted in Fig. 2, and has dimensions of 480 × 480 × 28 m.Its geological structure involves a network of fossilized meandering channels of high permeability.The average reservoir pressure is 40.0 MPa.All remaining geological and fluid properties used in this example are presented in Table 1.The reservoir model contains eight injection wells and four production wells.The near-wellbore flow is modelled using a Peaceman well model.
During the first 1.5 years of production from the reservoir, the bottomhole pressures of the producers are kept at a constant value of 39.5 MPa.During that time, the injection rates of all eight injectors are prescribed to fluctuate monthly with a uniform probability distribution around an average value of 5.52 × 10 −4 m 3 /s (300 bbl/day) and a Fig. 2 Three-dimensional oil reservoir model with eight injection and four production wells, after [21].Its geological structure involves a network of fossilized meandering channels of high permeability in a low-permeability background In this example, historic data are available over the first 1.5 years of production and lower and upper bounds on expected economic performance are determined over the remaining life of the field-from 1.5 to 6.0 yearsby changing model properties (grid block permeabilities), while the model stays compliant with historic data over the     this threshold is related to the ratio between oil revenues r o and water production costs r wp .To determine the historymatched models that provide the lower and upper bounds on the NPV for the remaining producing life, two hierarchical optimization procedures are initiated.They terminate when the feasible updates no longer result in a significant change in the NPV. Figure 3 depicts the measured production data, along with the simulated production data originating from the final lower-and upper-bound models resulting from the hierarchical optimization method.It shows that the errors between measured and simulated bottomhole pressures of the injectors and fractional flow rates of the producers are very small for both the lower-and upper-bound models.Thus, the condition that the updated models maintain a good history match is met.However, in Fig. 4, it can be observed that the permeability fields of both models are quite different.These differences have a large impact on the predicted production data given the assumed reactive production strategy, as can be observed in Fig. 5.Moreover, the change in permeability in the near-well areas around the injectors has a strong effect on the pressure response of the injectors.Finally, Fig. 6 shows the actual lower and upper bounds on the predicted NPV over time, in terms of the NPV for the entire producing reservoir life (6 years) and in terms of the incremental NPV for just the remaining (future) producing reservoir life (4.5 years).It can be observed that the upper and lower bounds of the incremental NPV are 63 % above and below their average value.

Brugge model example
In the second experiment, we use data from the Brugge benchmark workshop organized in 2009 [16,17].In the original benchmark study, the 'truth' case used to generate the data was not disclosed, and therefore, in this work, we use a new 'truth' honouring all well logs, geological descriptions and distributions of geological model parameters, porosity/permeability relations and the geological structure of the Brugge field.Figure 7 depicts the new 'true' Brugge permeability field, which is used to generate synthetic data.Blue and red bars in Fig. 7 represent injectors and producers respectively.The fluid properties and Corey exponents used in this example are given in Table 2.
The reservoir model consists of 60,048 active grid blocks and has dimensions of 3 km × 10 km × 80 m.It contains 11 injection wells located near the rim of the oil-water contact at a depth of 1678 m from the surface and 20 production wells, as depicted in Fig. 7. Wells are located in the grid block centres, and we use a standard Peaceman well inflow model.During the first 10 years of production (the history-matching period), all production wells are constrained to a minimum pressure of 4.9 MPa and a maximum liquid rate of 3.7 × 10 −3 m 3 /s and all injection wells operate at a constant water flow rate of 7.4 × 10 3 m 3 /s.Moreover, production wells are shut-in individually if the water fraction in the produced liquid is above 90 %.After the history-matching period (10 years), closed wells are reopened.Wells are drilled according to the time scheme presented in the Brugge workshop [17].

Historical data
In this example, historical data are available over the first 10 years of production and lower and upper bounds on expected economic performance are determined over the remaining life of the field-from 10 to 30 years-by changing the permeability field, while the model stays compliant with historic data over the first 10 years of production.Time-lapse seismic data as well as production data are used as historic data.Production data consist of periodic measurements of water and oil rates in the producers.Independent measurement errors are generated from Gaussian distributions with zero mean and standard deviations equal to 10 % of the original measurements.Negative production rates, after the addition of noise, are reset to zero.
Because the measurement errors are independent, the error covariance matrix is diagonal.

Multi-objective optimization settings
Using Eq. ( 13) as the primary objective function and Eq. ( 14) as the secondary objective function, two hierarchical optimization procedures are conducted to determine the history-matched models that provide the lower and upper bounds on the NPV for the remaining producing life.The procedures are terminated when the feasible updates no longer result in significant changes in objective function value.The starting point for the assisted historymatching process (primary objective function) is selected randomly out of 104 available prior models in the Brugge data set.The prior model is iteratively conditioned to historical data by adjusting the horizontal grid block permeability values.In this experiment, the water injection costs r wi , the water production costs r wp and the oil revenues r o are assumed constant at values of 5 $/bbl, −5 $/bbl and 80 $/bbl respectively.The discount rate, b, is set to 10 %.

Results: history matching based on production data
In this example, history matching is performed based on production data.We constrain the search space of the secondary problem by choosing the threshold value of Eq. ( 7) as 0.5 % of the minimum of the primary objective function.
Figure 8 depicts the historical data and the lower and upper bounds of water production in the first eight producers as an example of the typical ranges of the bounds.The historymatching and forecasting periods are separated by a dashed line.Blue and red colours represent the lower and upper bounds of oil and water production.Figure 9 depicts the injection pressures in the first four injectors.Unlike in the results for the previous example, depicted in Fig. 5, there is no jump in the pressures at the beginning of the forecasting period because they have already reached their maximum allowable values.Figure 10 depicts the historical data and the lower and upper bounds for the cumulative oil and water production of the entire field.
As can be seen in Fig. 10, the lower-bound and upperbound models produce the same history but different forecast.Moreover, Fig. 11 depicts the economic performance (NPV) of the upper-and lower-bound models over time for the entire production life, including the history and the prediction.In this experiment, the incremental NPV of the upper-bound model is 19.5 % higher than the incremental NPV of the lower-bound model.
Figure 12 shows the differences between the lower-and upper-bound permeability fields for all nine layers of the field.It can be observed in Fig. 10 that the permeability fields of both models are different, especially in the producing layers.These differences have an impact on the predicted production data while they result in the same production history, as can be observed in Fig. 8.We note that although the permeability values away from the wells are more likely to be in the null space (i.e. to have room for variation), they also have less of an effect on the output in the Fig. 8 Historical and predicted water production over 30 years of production for the first eight producers for the Brugge field example wells.Apparently, the optimization algorithm did not produce significant changes in these values because that would not have changed the resulting NPV.Computation of these results required 200 pairs of forward and backward (adjoint) simulations, where each pair took, on average, 786 s.

Effect of data type
In the previous section, we obtained lower-and upper-bound models based on production data.In order to investigate the effect of data type on the upper-and lower-bound models, two more experiments are conducted based on different data types.In the first experiment, the upper-and lower-bound models are obtained based on interpreted time-lapse seismic data (saturation maps) and production data.The saturation maps are generated by simulating the 'truth' and adding independent measurement errors by sampling from a Gaussian distribution with zero mean and standard deviations equal to 10 % of the simulated saturation values.As before, we constrain the search space of the secondary problem by choosing the threshold value of Eq. (7) as 0.5 % of the minimum of the primary objective function.The second experiment involves assimilation of both time-lapse seismic and production data while also prior information is added to the primary objective function as a regularization term.Figure 13 shows the incremental NPV differences between the lower-and upper-bound models obtained using different data types.As can be seen in Fig. 13, Fig. 9 Historical and predicted injection pressures over 30 years of production for the first four injectors for the Brugge field example the incremental NPV difference decreases by adding more information.

Effect of threshold value
In this section, we investigate the effect of the threshold value, ε, in Eq. (7).We constrain the search space of the secondary problem to different extents by choosing a range of threshold values varying between 0.15 and 1.5 % of the minimum of the primary objective function.Interpreted time-lapse seismic data (saturation maps) and production data formed the historical data, and two hierarchical multiobjective optimizations were conducted to find the lower and upper bounds for the reservoir model for different threshold values.Fig. 14 shows the incremental NPV difference between the upper-and the lower-bound models versus the threshold value ε.
Figure 15 depicts how the primary and secondary objective functions change for different values of ε.Figures 14  and 15 show that as the threshold value in Eq. ( 7) increases, the difference between the lower and upper model values of the incremental NPV increases also.However, the effect is not very large and even for the lowest threshold value (ε = 0.15 %), a difference of approximately 17 % in incremental NPV is obtained.
We note that the lower and upper bounds have been obtained by a gradient-based optimization technique which may have resulted in local rather than global optima.Lower lower bounds and higher upper bounds may therefore exist.
We also note that both the red and the blue curves can be interpreted as parts of (approximate) Pareto curves.Points on a Pareto curve are at the boundary of the feasible set of solutions in the bi-objective space, and recently, several studies have been performed to characterize full Pareto curves for bi-objective flooding optimization; see, e.g.[13].Such a curve gives the decision maker the opportunity to select between competing objectives, i.e. to achieve a large value of the secondary objective function at the price of a strong drop in the primary objective function value or a somewhat smaller secondary objective function value without losing much of the primary objective function value.

Discussion
In the current paper, we used grid block permeabilities as history-matching parameters.However, the proposed method could equally well be applied using other parameters, e.g.porosities, fault multipliers or aquifer strength.Moreover, other data types than the production data and interpreted time-lapse seismic that we used could be assimilated.We note that the use of gradients with respect to the history-matching parameters is an important ingredient in our method.This implies that we need a technique to compute those gradients.We used an adjoint method, which is computationally very efficient.However, it is, in theory, also possible to implement our method using approximate gradients obtained with, e.g. the simultaneous perturbation stochastic approximation (SPSA) technique (see [19]) or ensemble optimization (EnOpt) (see [3] for the basics of the method and [8] for an implementation in hierarchical optimization).The latter (EnOpt) approach also allows for the inclusion of uncertainty in the reservoir models; see [9].
We note that our method has theoretical links to the use of level sets to relax the primary objective function  constraint in hierarchical optimization as discussed in the last paragraph of Section 2.Moreover, we note that other weak-constrained optimization methods could be applied to solve the hierarchical optimization problem.The current implementation has shown to be robust in various applications using both adjoint-based and ensemble-based techniques [4,8,23].
The determination of lower and upper bounds of future production using different types of data, as performed in the Brugge example, can be interpreted as a means to assess the cost-effectiveness of acquiring different data types to reduce the uncertainty in the expected NPV.It is tempting to interpret this as a way to assess the value of information (VOI) of those measurements, but because we do not know the statistical properties of the forecasted NPV, we cannot draw conclusions about the change in expected value of those forecasts and therefore our method does not truly provide the VOI.(For detailed information about the concept of VOI, see [2,6])

Conclusions
In this paper, we presented a hierarchical optimization method to determine lower and upper bounds on predicted production from history-matched models.We conclude that: • The nonuniqueness of history-matched models implies that future production can only be predicted within bounds.• The nonuniqueness implies the presence of remaining degrees of freedom after history matching (i.e. after solving the primary optimization problem) which can be used to determine lower and upper bounds on future production through solving two secondary optimization problems.
• The method proposed in our paper provides a way to gain more insight in the possible economic consequences of the lack of information in historic data.These consequences can be represented by total production, ultimate recovery, (incremental) NPV or any other economic measure.• The method is not limited to historic production data.
Alternative data sources, e.g.time-lapse seismic data, can be used to determine the impact on the predicted economic performance.• Introducing more data sources, e.g.time-lapse seismic or prior information, results in smaller differences in economic performance (incremental NPV) between the lower-and upper-bound models.
appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Fig. 3 Fig. 4
Fig.3 Measured production data of the first 1.5 years of production from the (synthetic) 3D reservoir, along with the simulated production data originating from the lower− and upper−bound models

Fig. 5 Fig. 6
Fig.5 Measured production data of first 1.5 years of production from the (synthetic) 3D reservoir, along with the simulated production data for the remaining 4.5 years of production until the end of the field's life, originating from the lower-and upper-bound models

Fig. 7
Fig. 7 Permeability field with 11 injection wells and 20 production wells.The blue surface indicates the oil-water contact

Fig. 10
Fig. 10 Historical and predicted cumulative oil production (left) and water production (right) over 30 years of production

Fig. 12
Fig. 12 Difference between the lower-and upper-bound permeability fields.All permeability values are expressed as the natural logarithm of permeability in mD

Fig. 13
Fig. 13 Difference between the upper-bound and lower-bound incremental NPV values for models obtained based on different data types

Fig. 14
Fig. 14 Incremental NPV difference between the upper-bound and lower-bound models for different epsilon values

Fig. 15
Fig.15 Secondary objective function value versus its corresponding primary objective function value, both expressed as incremental NPV

Table 2
Fluid properties and Corey exponents for the Brugge field example