Reservoir Simulator Runtime Enhancement Based on a Posteriori Error Estimation Techniques
Amélioration des performances des simulateurs de reservoir basée sur des techniques d’estimateur d’erreur a posteriori
* Corresponding author
Accepted: 25 May 2016 Revised: 4 May 2016
In this work, we show how the a posteriori error estimation techniques proposed in [Di Pietro et al. (2014) Computers & Mathematics with Applications 68, 2331-2347] can be efficiently employed to improve the performance of a compositional reservoir simulator dedicated to Enhanced Oil Recovery (EOR) processes. This a posteriori error estimate allows to propose an adaptive mesh refinement algorithm leading to significant gain in terms of the number of cells in mesh compared to a fine mesh resolution, and to formulate criteria for stopping the iterative algebraic solver and the iterative linearization solver without any loss of precision. The emphasis of this paper is on the computational cost of the error estimators. We introduce an efficient computation using a practical simplified formula that can be easily implemented in a reservoir simulation code. Numerical results for a real-life reservoir engineering example in three dimensions show that we obtain a significant gain in CPU times without affecting the accuracy of the oil production forecast.
Dans ce travail, nous montrons comment les techniques d’estimateur d’erreur a posteriori proposées dans la référence [Di Pietro et al. (2014) Computers & Mathematics with Applications 68, 2331-2347] peuvent être efficacement utilisées pour améliorer les performances d’un simulateur compositionnel de réservoir dédié aux procédés de Récupération Assistée des Hydrocarbures (RAH). Cet estimateur d’erreur a posteriori permet de proposer un algorithme de raffinement de maillage adaptatif conduisant à un gain significatif en termes de nombre de cellules du maillage par rapport à une résolution avec un maillage fin, et de formuler des critères d’arrêt pour les solveurs itératifs algébriques et non linéaires sans perte de précision. Cet article se concentre sur le coût de calcul des estimateurs d’erreur. Nous introduisons un calcul efficace en utilisant une formulation simplifiée qui peut facilement être mise en oeuvre dans un code de simulation de réservoir. Des résultats numériques pour un vrai cas d’étude 3D montrent que nous obtenons un gain significatif de temps CPU, et ce sans affecter la précision des prévisions de production du champ.
© J.-M. Gratien et al., published by IFP Energies nouvelles, 2016
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Reservoir simulation models are nowadays a key element for oil and gas companies in the development of fields. Reservoir simulation is used extensively to design advanced EOR processes like miscible displacement by CO2 or nitrogen, chemical flooding (polymer, alkaline, surfactant, etc.) or thermal processes with in situ combustion or water steam injection, see Lake . Such models with very complex physics are usually time-consuming and expensive as they lead to solve strongly non linear numerical models. To preserve the predictivity and the accuracy of the models with reasonable running simulation time, it is therefore crucial to improve the runtime performance of the simulators.
A posteriori estimate aims at controlling the error during the simulation and enables to develop adaptive algorithms leading to a significant computational saving. First rigorous approachs can be found in [2-4] suitable for unsteady nonlinear models, and [5-7] for degenerated problems. Recently, abstract frameworks for a posteriori estimates appeared, for two-phase flow [8, 9], multiphase compositional model , and the thermal multiphase compositional global model . In these frameworks the dual norm of the residual augmented by a nonconformity evaluation term is controlled by fully computable estimators decomposed into space, time, linearization, and algebraic error components. This allows to formulate criteria for stopping the iterative algebraic solver and the iterative linearization solver when the corresponding error components do not affect significantly the overall error. Moreover, the spatial and temporal error components can be balanced by adaptations of both time step and space mesh. The aforementioned approach is promising. Significant gains in terms of the total number of algebraic solver iterations and important reductions in terms of the number of cells can be obtained in comparison with standard options of reservoir simulators. However, no straightforward comparison of the computational cost has been measured on reservoir engineering practices.
In this paper we propose an engineering approach of the a posteriori error estimates, provided by recent works [10, 11], in order to establish an efficient computation of the estimators in such a way that we can ensure important computational savings for realistic complex models. Our approach is based on providing a low-cost evaluation formula by avoiding flux reconstructions. As mentioned before we rely here on the general ideas of [10, 11], we apply a posteriori error estimators allowing to distinguish the different components of the global error, and we use them to formulate stopping criteria for the iterative algebraic and nonlinear solvers together with refinement/coarsening criteria for both the time step and the space mesh. The present work brings two major novelties: first, we propose a well-implemented evaluation of the estimators using a practical simplified formula making them usable in engineering practices; second, we apply the stopping criteria for the iterative algebraic solver in the context of a three-dimensional realistic reservoir simulation study.
The paper is organized as follows. In Section 1 we describe a model problem, the two-phase flow model, and we discuss its numerical resolution. In Section 2 we state the corresponding a posteriori error estimate. Then, we introduce in Section 3 a practical simplified formula for the evaluation of the estimators. In Section 4 we discuss the use of a posteriori error estimation tools to enhance the runtime reservoir simulator performance. Finally, in Section 5 we validate our approach on a real-life reservoir engineering example in three dimensions.
We introduce in this section a two-phase flow model, its finite volume discretization, the linearization of the resulting discrete system and the algebraic resolution. In [10, 11] a posteriori error estimates are driven for the thermal multiphase compositional general model. However, in this work we focus on the two-phase flow since it is the model of our chosen study case (SPE10 case, see Sect. 5), and for the clarity of the simplified evaluation process.
We consider an immiscible two-phase model where the phases are water and oil. Each of them is composed of only one component, water component and hydrocarbon component, respectively. As we have an immiscible model we will use the same index to represent the phase and its component. Consequently, we use the indices w,o to represent the water and the oil, respectively.
where is the capillary pressure. Consequently, the unknowns of the model are the water pressure and the saturations , with . These unknowns are collected in the vector . We denote by the porosity of the medium and by the permeability tensor. For a phase is the relative permeability, is the dynamic viscosity and is the density. The system of governing equations is given by(2)
where denotes the gravity vector acting in the negative z direction and g its Euclidian norm.
To discretize our model we choose a popular scheme in the oil industry, the implicit in time cell centered finite volume scheme, that is fast and simple to implement and well-known by its stability.
In order to define this numerical scheme we start by giving some notations for the space-time mesh.
Let denote a sequence of positive real numbers corresponding to the discrete time steps such that . We consider the discrete times such that and, for ; then we define the time intervals . Let denote a family of meshes of the space domain superadmissible in the sense of Eymard et al. [12, Definition 3.1]. In what follows, for all , we denote by its -dimensional Lebesgue measure and by its diameter. For all ; the distance between the cell center and the face center , where denotes the faces of an element not lying on . Finally, we set the vertices of , and the elements from sharing the vertex .
The cell centered finite volume scheme is a two-point flux approximation scheme, where the unknowns of the model are discretized using one value per cell: For all we let , with . System (2) is then discretized as follows: For all , and , we require(4)with denoting the upstream mobility and the two-point finite volume approximation of the normal component of the average phase velocity on is given by(5)where and is an approximation of the mass density of the phase on the face . Boundary fluxes are set to zero to account for the homogeneous natural boundary condition.
At this stage, we need to solve, at each time step, the system of nonlinear algebraic equations resulting from the discretization method of Section 1.4. First, for we apply the Newton linearization algorithm generating, for fixed, a sequence with solution to the following system of linear algebraic equations: For all and all ,(6)with, , we have
Second, for , and a given Newton iteration , to approximate the solution of System (6), we use an iterative algebraic solver generating, for fixed, a sequence solving Equation (6) up to the residuals, given for all and all by(7)with linearized phase fluxes(8)
The spirit of the a posteriori error estimates is to monitor the computational error. In [10, 11] we derive fully computable a posteriori error estimates of the dual norm of the residual augmented by a nonconformity evaluation term for the thermal multiphase compositional general model. Moreover, in  the authors bound an energy-type norm of immiscible incompressible two-phase flow.
In the previously mentioned articles, similar types of estimators (more simple for our model) are obtained, whether we choose to bound the energy-type norm or the dual norm of residual. Therefore, we will use the notation for the error norm between the known numerical approximation and the unknown exact solution, without giving more precision. The local-in-time error norm is denoted then by at time step , Newton iteration , and linear solver iteration .
As mention in the introduction, our main goal is to propose a simple evaluation of the estimators proposed in  in order to be usable for engineering practices. Accordingly, we will directly apply these a posteriori error estimates on our model in order to obtain an a posteriori error estimate distinguishing the space, time, linearization, and algebraic errors. To this end, first we have to recall the required reconstructions necessary to get the upper bound, see  for more details.
We recall here the necessary reconstructions for the theory of the a posteriori error estimates . For these reconstructions we need the Raviart–Thomas–Nédélec (RTN) finite-dimensional subspaces of . Since the meshes , used in the numerical experiments of Section 5 below, consist of rectangular parallelepipeds then we define the RTN subspaces bysee Brezzi and Fortin  for details. We recollect all the formulas since they are necessary for the evaluation process of Section 3.
For all we define such that, for all and all ,(9)with defined by (5) and on . Then, for each , we introduce the piecewise quadratic phase pressure such that, for all ,(10)and The space–time function is then continuous and piecewise affine in time, given by at the discrete times . We need also phase pressure reconstructions, , such that for all . These reconstructions are typically piecewise polynomial continuous in space and piecewise affine continuous in time.
with defined by (5), while on . We also define a linearization error flux reconstruction such that, for all and for all ,(11b)with defined by (8). Similarly, an algebraic error flux reconstruction is defined such that, for all and for all ,(11c)with defined by (7). To complete both (11b) and (11c), we set respectively and on . The total flux reconstruction is then given by(11d)
Following , under the assumption of the existence and the uniqueness of a weak solution of the considered problem and the reconstructions of phase pressures and fluxes evoked in Section 2.1.2, the local-in-time error measure can be bounded as follows(12)with
In order to implement the different estimators in a reservoir simulation code, we propose here to simplify them. Our simplification is based on the following observations:
our estimators are expressed in terms of norms of various discrete quantities. Therefore we only need to compute these norms and not necessarily the quantities themselves;
all these norms are integrals. Therefore they can be computed/approximated on each adequate quadrature formulas by mesh element;
we thus only need to know the values of the various quantities at the quadrature points;
the implementation of RTN spaces and the physical flux reconstructions can be consequently avoided.
Let us now describe our simplification process in all details.
The a posteriori estimates of Section 2.2 use the reconstruction of the different fluxes and in the space (Sect. 2.1.2). Furthermore, the pressure reconstructions also involve spaces (Sect. 2.1.1) in a way that for all ,
Hence, a priori, implementing these estimators requires operations with RTN spaces as seen in the previous chapter. To the best of our knowledge, so far RTN spaces are not implemented in petroleum industrial codes. Therefore we search for a simplification of these estimators that avoids the use of RTN spaces.
Let a mesh cell be given and let . Then can be expressed by(18)where are the degrees of freedom on the face , i.e., the face fluxes , and are the lowest-order Raviart–Thomas–Nédélec basis functions, given by (for parallelepiped meshes such as the ones used in the test cases of Sect. 5 below)where is the barycenter of the face σ′ opposite to the face , and(19)with,
Here , and . We can also rewrite this formula with vector symbols by(20)with , and as given in Table 1.
To apply the formula (20) to the integration of a function over a cell of our parallelepiped mesh, we consider the transformation such that,(21a) (21b) (21c)
As a conclusion, to compute an norm of reconstructed flux functions in the space , we just need to obtain the degrees of freedom , represented by the face normal fluxes. These can be obtained directly without any flux reconstruction. In what follows we use these simplifications to evaluate the estimators for the immiscible two-phase model.
We begin with the spatial estimators. Recall the reconstructions of the conservative fluxes and given by Equation (11d). For , one has
Owing to the reconstruction of in Equation (11a), in (11b), and in (11c), we satisfy the conservation propertyup to a neglected maladjustment from the practical construction of (Corollary 3.4). Then, the residual estimators and given by Equation (14) will be neglected. Note that, however, if we construct the algebraic flux following Equation (11c) then the residual estimators and are zero.
The non conformity estimators are also evaluated in the same way. These estimators involve the difference of two terms: and In order to apply formula (22), these two terms must be expressed in the space. This is the case of the first term (by construction, see Eq. 10), but not the case of the second . Therefore, we lift this term into the space by preserving the corresponding normal fluxes over the faces: Let be the piecewise linear nodal basis functions of on the three-dimensional element . Then for all the continuous pressure can be expressed as
To alleviate the industrial implementation, we then approximate the value of the pressure in the nodes by(25)where stands for the cardinality of the set ; it is equal to 8 for interior vertices in the three-dimensional case. But the values of are preferable since they are more precise. The normal fluxes of over the faces of can then be expressed explicitly by: for all (26)
This leads finally to a simple evaluation of the local spatial estimators for all . Note that a more precise formula could be obtained by replacing the piecewise linear nodal basis functions of by the piecewise quadratic basis functions of .
Summarizing the above developments, we have:
Finally to evaluate the global spatial estimator we approximate the time integral by a one-dimensional integration formula.
Now for the temporal estimators, recall that for all , the reconstructions of the phase pressures , are such that are in the RTN space. Thus, the evaluation of the local temporal estimators can be again done using formula (22).
Corollary. 3.2 (simple formula to evaluate the temporal estimators). By applying formula (22) to Equation (16) and using the reconstructions formulations (11a)–(11c), we can approximate the temporal estimators as (29) where and are given by (21) and (19),
For the linearization estimators, we remark that the fluxes which compose the estimators are the reconstructed fluxes in the RTN space: , given by Equation (11b). Thus evaluating of these local estimators is straightforward using formula (22):
Corollary 3.3 (simple formula to evaluate the linearization estimators). By applying formula Equation (22) in Equation (17) and using the reconstructions formulations (11b), we can approximate the local linearization estimators as (30) where and are given by (21) and (19), respectively.
Finally, for the algebraic estimators, the evaluation of the local norms is also straightforward using formula (22), as the contributions given by (11c), involved in these norms, are in the RTN space.
Following [10, Remark 3.3], we arrive at the following approximation formula:
Corollary 3.4 (computing practically the algebraic error). Following [10, Remark 3.3], to compute approximately the algebraic error, we perform j additional iterations of the algebraic solver from the stage (7) where j is a user-defined fixed number. Then, the local algebraic error estimators can be evaluated using formula (22) as follow:(31a) where and are given by (21) and (19), respectively.
As a conclusion, the evaluation of the different estimators can be carried out while avoiding the physical reconstructions of fluxes in the space. The key is the use of a quadrature formula for computing the norms, knowing the normal fluxes over all faces of any cell . Using this technique greatly simplifies the implementation of the estimators (in particular into industrial legacy codes) and yields an important computational saving compared to the previous technique where we need to build the RTN flux reconstructions.
The presented a posteriori error estimate framework has been implemented in a parallel reservoir prototype simulator , a thermal multi-purpose simulator with the capability to simulate various EOR processes like water injection, steam injection, etc. It is part of the next generation IFP Energies Nouvelle research simulators, based on Arcane framework . To improve the runtime performance of our simulator, we elaborated two strategies: the first consists in improving our linear solver by reducing the number of iterations using adaptative linear stopping criteria based on the a posteriori error estimate framework; the second strategy consists in reducing the global number of degrees of freedom by reducing the number of grid blocks of the mesh combining the Arcane adaptive mesh refinement feature to an advanced local space a posteriori error estimator.
We have implemented an adaptative linear solver stopping criteria which allows by a call back mechanism to evaluate regularly the a posteriori linear error estimate described in the previous section. This mechanism allows to stop iterative linear solvers once the relative linear stopping criteria has been reached, or when the linear error estimation is lower than the non linear error estimation by a given factor. We can, in that way, reduce the number of linear solver steps during the steps of the Newton iterative non linear solver.
In our simulator, we developed an adaptative mesh refinement mechanism that allows to refine or coarsen cells based on a space error criteria. This mechanism is aimed at adapting the grid block size to optimize the total number of degrees of freedom for a given global discrete error criteria. The idea is to have coarse grid blocks where the local error criteria is low, and fine grid block where the local error criteria is high. As the performance of the simulator is directly linked to the total number of grid blocks, minimizing this number has a direct impact on the global simulator runtime performance. We implemented a space error criteria based on the a posteriori error estimate tools described in the previous section. The idea is to reduce the standard deviation of a computed error field, used to choose the cells to refine and the cells to coarsen. The adaptative mesh refinement mechanism is then applied to minimize the required number of degrees of freedom for a given error criteria. During the coarsening and refinement phases, we use some upscaling and interpolating mechanisms to ensure the mass balance. Geological data (rock model properties like porosity, permeability, etc.) are stored at the finest geological level as data on a cloud of space points. These data are used, during mesh adaptation by standard upscaling algorithms to compute the values of the rock model properties on new created cells .
We validate our approach on the tenth SPE comparative solution project model . It is an incompressible water-oil, black-oil, two-phase flow problem. It is built on a Cartesian regular geometry with no top structure or faults and simulates a 3D waterflood in a geostatistical model. In this five-spot production model one water injection well is located at the center of the reservoir and four production wells at the four corners of the reservoir grid. Producers are controlled by a bottom hole pressure target. We present some results obtained on the fine grids comparing the runtime performance between the standard options of the simulator and the adaptative linear solver stopping criteria presented in Section 4.1. Others runs experiment the adaptative mesh refinement feature with the space criteria presented in Section 4.2 and results are compared against fine and coarse grid simulations. The coarse grid is obtained from the fine reference grid after two levels of coarsening applied on all cells. To validate the simulation results, we compared production curves for different wells to ensure that the results are still valid from a reservoir engineering point of view.
We have run our test cases using 16 cores of a dual socket linux workstation with 2 Sandy Bridge 10-cores Intel® Xeon® E5-2680 CPU with a tact rate of 2.8 GHz, 25600 KB cache size.
To have a better control of the linear solver, we applied here a stopping criteria similar to that of the adaptive algorithm proposed in  with the main difference being that we avoid the physical flux reconstructions and the evaluation of the estimators is based on the simplification proposed in Section 3. We consider here a simulation of 200 days. In Figure 1, we compare the oil production rate from both the resolution with the standard options and then with the stopping criteria. In the left part of Figure 1 we depict cumulative oil production for the first two production wells. Cumulative oil production for the third and fourth production wells are then presented in the right part of Figure 1. We observe that the stopping criteria does not have any effect on the accuracy of the well productions.
Cumulative oil production of wells during the simulation.
Regarding the numerical behavior of the previous two runs we have had the following results:
with standard options the number of resolutions related to the time steps is 67, the total number of non linear steps is 272, the total number of linear solver steps is 82 053, and the global CPU time of this simulation is 4 653 seconds;
with the stopping criteria the number of resolutions related to the time steps is 67, the total number of non linear steps is 272, the total number of linear solver steps is 56 732, and the global CPU time of this simulation is 2 754 seconds.
As it is obvious, significant saving is obtained in terms of the cumulated number of linear solver iterations. We further illustrate in Figure 2 a comparison in terms of the consumed CPU time as well as the cumulated CPU time at several time steps of the simulation. The results show an over-all gain in CPU time, with a speed-up factor reaching the value 2.
CPU time (left) and total CPU time (right) of the simulation.
We evaluate the performance of our simulator using the AMR feature combined with the a posteriori space error criteria. We compare in Figure 3 the oil production rate of the producer wells resulting from the resolution on the fine mesh, the coarse mesh with and without the AMR option. The fine grid results are supposed to be the reference. We observe that for all wells except well1 the AMR simulation behaves identically to the fine grid and much better than the simulation on the coarse grid.
Cumulative oil production of wells during the simulation.
In Figure 4 we compare the CPU time and the number of active cells regarding the simulation time. We observe a gain of performance for each time step for the same accuracy than the resolution on the fine mesh, due to a reduce number of required degrees of freedoms intermediate between the number of cells of the fine and the coarse mesh. We can notice some peaks of CPU time. This is due to the fact that the time step management is not optimal with the AMR feature, there are more failed time steps than in the reference run. We need more investigation to better understand this difference of behaviour in our simulator time step management.
CPU time (left) and active cell number (right) of the simulation.
With our new a posteriori error estimation framework, we have improved the performance of our compositional reservoir simulator dedicated to EOR processes. We have first optimized our non linear solver by implementing an adaptative linear stopping criteria to avoid to over-solve the linear systems in the first non linear steps. We have implemented secondly a space criteria to better manage the use of the AMR feature by optimizing the number of grid blocks required for a given accuracy of the well production results. Our approach has been validated on the SPE10 benchmark, a 3D study case representative of real-live reservoir study cases. Indeed, we have obtained significant gain in CPU times without affecting the accuracy of the oil production forecast. However, we still need to study the effect of the AMR feature on our dynamic time step management that we may improve in future works, implementing new criteria based on time error estimation.
- Lake L.W. (1989) Enhanced Oil Recovery. Prentice Hall Inc., Old Tappan, NJ.
- Eriksson K., Johnson C. (1995) Adaptive finite element methods for parabolic problems, IV. Nonlinear problems, SIAM J. Numer. Anal. 32, 6, 1729–1749.
- Verfürth R. (1998) A posteriori error estimates for nonlinear problems: Lr(0, T;W1,ρ(Ω))-error estimates for finite element discretizations of parabolic equations, Numer. Methods Partial Differential Equations 14, 4, 487–518. [CrossRef] [MathSciNet]
- Verfürth R. (1998) A posteriori error estimates for nonlinear problems. Lr(0, T; Lρ(Ω))-error estimates for finite element discretizations of parabolic equations, Math. Comp. 67, 224, 1335–1360. [CrossRef] [MathSciNet]
- Nochetto R.H., Schmidt A., Verdi C. (2000) A posteriori error estimation and adaptivity for degenerate parabolic problems, Math. Comp. 69, 229, 1–24. [CrossRef] [MathSciNet]
- Ohlberger M. (2001) A posteriori error estimate for finite volume approximations to singularly perturbed nonlinear convection–diffusion equations, Numer. Math. 87, 4, 737–761. [CrossRef] [MathSciNet]
- Di Pietro D.A., Vohralík M., Yousef S. (2015) Adaptive regularization, linearization, and discretization and a posteriori error control for the two-phase stefan problem, Math. Comput. 84, 291.
- Cancès C., Pop I.S., Vohralík M. (2013) An a posteriori error estimate for vertex-centered finite volume discretizations of immiscible incompressible two-phase flow, Math. Comp. doi: 10.1090/S0025-5718-2013-02723-8.
- Vohralík M., Wheeler M.F. (2013) A posteriori error estimates, stopping criteria, and adaptivity for two-phase flows, Computational Geosciences 17, 5, 789–812. [CrossRef] [MathSciNet]
- Di Pietro D.A., Flauraud E., Vohralík M., Yousef S. (2014) A posteriori error estimates, stopping criteria, and adaptivity for multiphase compositional Darcy flows in porous media, J. Comput. Phys. 276, 163–187. [CrossRef]
- Di Pietro D.A., Vohralík M., Yousef S. (2014) An a posteriori-based, fully adaptive algorithm with adaptive stopping criteria and mesh refinement for thermal multiphase compositional flows in porous media, Computers & Mathematics with Applications 68, 12, Part B, 2331–2347. [CrossRef] [MathSciNet]
- Eymard R., Gallouët T., Herbin R. (2000) The finite volume method, Vol. 7, Handbook of Numerical Analysis, P.G. Ciarlet and J.-L. Lions (eds), North Holland.
- Brezzi F., Fortin M. (1991) Mixed and hybrid finite element methods, Vol. 15 of Springer Series in Computational Mathematics, Springer-Verlag, New York, ISBN 0-387-97582-9. [CrossRef]
- Mesri Y., Gratien J.-M., Ricois O.M., Gayno R., et al. (2013) Parallel adaptive mesh refinement for capturing front displacements: Application to thermal eor processes, in SPE Reservoir Characterization and Simulation Conference and Exhibition, Society of Petroleum Engineers. doi: 10.2118/166058-MS.
- Grospellier G., Lelandais B. (2009) The arcane development framework, in Proceedings of the 8th workshop on Parallel/High-Performance Object-Oriented Scientific Computing, POOSC ’09, pp. 4:1–4:11, New York, NY, USA.
- Mesri Y., Ricois O. (2015) Construction process for an improved meshing for the simulation of a reservoir in an underground formation, URL http://brevets-patents.ic.gc.ca/opic-cipo/cpd/fra/brevet/2886110/
- Christie M.A., Blunt M.J. (2001) Tenth SPE Comparative Solution Project:A Comparison of Upscaling Techniques. Reservoir Simulation Symposium.
Cite this article as: J.-M. Gratien, O. Ricois and S. Yousef (2016). Reservoir Simulator Runtime Enhancement Based on a Posteriori Error Estimation Techniques, Oil Gas Sci. Technol 71, 59.