ECMOR XI - 11th European Conference on the Mathematics of Oil Recovery

Although many advanced methods have been devised for the hyperbolic transport problems that arise in reservoir simulation, the most widely used method in commercial simulators is still the implicit upwind scheme. The robustness, simplicity and stability of this scheme makes it preferable to more sophisticated schemes for real reservoir models with large variations in flow speed and porosity. However, the efficiency of the implicit upwind scheme depends on the ability to solve large systems of nonlinear equations effectively. To advance the solution one time step, a large nonlinear system needs to be solved. This is a highly non-trivial matter and convergence is not guaranteed for any starting guess. This effectively imposes limitations on the practical magnitude of time steps as well as on the number of grid blocks that can be handled. In this paper, we present an idea that allows the implicit upwind scheme to become a highly efficient. Under mild assumptions, it is possible to compute a reordering of the equations that renders the system of nonlinear equations (block) lower triangular. Thus, the nonlinear system may be solved one (or a few) equations at a time, increasing the efficiency of the implicit upwind scheme by orders of magnitude. Similar ideas can also be used for high-order discontinuous Galerkin discretizations. To demonstrate the power of these ideas we show results and timings for incompressible and weakly compressible transport in real reservoir models.

In the recent years of high oil prices many companies turned their attention to the advanced IOR/EOR methods including cyclic water flooding, miscible gas and steam injection. This paper summarizes experience accumulated in International Research Institute of Stavanger over several years of using in-house analytical pre-screening tool for evaluation of possible IOR/EOR strategies. An analytical models of cyclic water, miscible gas and steam injection is presented. Strong and week points of these methods are analyzed. The results of field case studies both on the Norwegian Continental Shelf and outside Norway are presented and discussed. We conclude that the analytical models may be successfully used for prescreening of the IOR methods. The analytical methods can, in most cases, correctly handle the behavior of the reservoir system as a function of IOR process critical parameters. However final decision to apply a certain IOR procedure and its optimization should be done utilizing full-scale reservoir models. The analytical methods can significantly decrease the decision making time by either: a) decreasing a number of computationally expensive full-scale simulations for mature developed reservoirs; or b) quickly screening new prospect where reservoir knowledge is limited.

The method of characteristics, or fractional-flow theory, is extremely useful in understanding complex IOR processes and in calibrating simulators. One limitation has been its restriction to Newtonian rheology except in rectilinear flow. Its inability to deal with non-Newtonian rheology in polymer and foam IOR has been a serious limitation. We extend fractional-flow methods for two-phase flow to non-Newtonian fluids in one-dimensional cylindrical flow, where rheology changes with radial position r. The fractional-flow curve is then a function of r. At each position r the "injection" condition is the result of the displacement upstream; one can plot the movement of characteristics and shocks downstream as r increases. We show examples for IOR with non-Newtonian polymer and foam. For continuous injection of foam, one can map the 1D displacement in time and space exactly for the non-Newtonian foam bank, and approximately, but with great accuracy, for the gas bank ahead of the foam. The fractional-flow solutions are more accurate that finite-difference simulations on a comparable grid and can be used to calibrate simulators. Fractional-flow methods also allow one to calculate changing mobility near the injection well to a greater accuracy than simulation. For polymer and SAG (alternating-slug) foam injection, characteristics and shocks collide, making the fractional-flow solution complex. Nonetheless, one can solve exactly for changing mobility near the well, again to greater accuracy than with conventional simulation. For polymer solutions that are Newtonian at high and low shear rates but non-Newtonian in between, the non-Newtonian nature of the displacement depends on injection rate. The fractional-flow method extended to non-Newtonian flow can be useful both for its insights for scale-up of laboratory experiments and to calibrate computer simulators involving non-Newtonian IOR. It can also be an input to streamline simulations.

This paper focuses on flux-continuous pressure equation approximation for strongly anisotropic media. Previous work on families of flux continuous schemes for solving the general geometry-permeability tensor pressure equation has focused on single-parameter families e.g. [1]. These schemes have been shown to remove the O(1) errors introduced by standard two-point flux reservoir simulation schemes when applied to full tensor flow approximation. Improved convergence of the schemes has also been established for specific quadrature points [2]. However these schemes have conditional M-matrices depending on the strength of the cross-terms [1]. When applied to cases involving full tensors arising from strongly anisotropic media, these schemes can fail to satisfy the maximum principle. Loss of solution monotonicity then occurs at high anisotropy ratios, causing spurious oscillations in the numerical pressure solution. New double-family flux-continuous locally conservative schemes are presented for the general geometry-permeability tensor pressure equation. The new double-family formulation is shown to expand on the current single-parameter range of existing schemes that have an M-matrix. While it is shown that a double family formulation does not lead to an unconditional M-matrix scheme, an analysis is presented that classifies the sensitivity of the new schemes with respect to monotonicity and double-family quadrature. Convergence performance with respect to a range of double-family quadrature points is presented for some well known test cases.

The SAGD (Steam Assisted Gravity Drainage) is an enhanced oil recovery process for heavy oils and bitumens. To simulate this thermal process and obtain precise forecasts of oil production, very small-sized cells have to be used because of the fine flow interface. But the use of fine cells throughout the reservoir is expensive in terms of CPU time. To reduce the computation time one can use an adaptive mesh refinement technique which will use a refined grid at the flow interface and coarser cells elsewhere. The first numerical tests of SAGD modelling, using a criterion based on two threshold temperatures, showed that this adaptive-mesh-refinement technique could reduce significantly the number of cells [LLR03]. Nevertheless, the refined zone remains too wide. In this work, we introduce a new criterion for the use of adaptive mesh refinement in SAGD simulation. This criterion is based on the work achieved in [KO00] on a posteriori error estimators for finite volume schemes for hyperbolic equations. It enables to follow more precisely the flow interface. Through numerical experiments we show that the new criterion enables to further decrease the number of cells (and thus CPU times) while maintaining a good accuracy in the results.

Modern reservoir simulation models often involve geometrically complex (distorted or unstructured) three-dimensional grids populated with full-tensor permeability coefficients. The solution of the flow problem on such grids demands the use of accurate spatial discretizations such as the multipoint flux approximation (MPFA) finite volume schemes or mixed finite element methods (MFEM). The extension of the streamline method to modern grids therefore requires a streamline tracing algorithm adapted to these advanced discretizations. In this paper, we present a new algorithm to trace streamlines from MPFA or MFEM on general tetrahedral or hexahedral grids and in the presence of full-tensor permeabilities. Our approach was already used successfully in 2D and this paper presents the extension to 3D of the previous work by the authors. The method is based on the mathematical framework of MFEM. Since MPFA schemes have recently been interpreted as MFEM, our streamline tracing algorithm is also applicable for the MPFA finite volume methods. Using the mixed finite element velocity shape functions, the velocity field is reconstructed in each grid cell by direct interpolation of the MFEM velocity unknowns or MPFA subfluxes. An integration of the velocity field to arbitrary accuracy yields the streamlines. The method is the natural extension of Pollock’s (1988) tracing method to general tetrahedral or hexahedral grids. The new algorithm is more accurate than the methods developed by Prévost (2003) and Haegland et al. (2007) to trace streamlines from MPFA solutions and avoids the expensive flux postprocessing techniques used by both methods. After a description of the theory and implementation of our streamline tracing method, we test its performance on a full-field reservoir model discretized by an unstructured grid and populated with heterogeneous tensor coefficients.

Run-time performance of a reservoir simulator is significantly impacted by the selection of the linear solver preconditioner, iterative method and their adjustable parameters. The choice of the best solver algorithm and its optimal parameters is a difficult problem that even the experienced simulator users cannot adequately solve by themselves. The typical user action is to use the default solver settings or a small perturbation of them that are frequently far from optimal and consequently the performance may deteriorate. There has been extensive research to develop automatic performance tuning and self-adaptive solver selection systems. For example Self-Adapting Large-scale Solver Architecture (SALSA) developed at the University of Tennessee requires running of large number of problems to initialize the system before using it. In contrast we propose an adaptive control on-line system to optimize the simulator performance by dynamically adjusting the solver parameters during the simulation. We start with a large set of parameters and quickly choose the best combinations that are continuously adapted during the simulation using the solver runtime performance measurements (e.g. solver CPU time) to guide the search. This software system, called the Intelligent Performance Assistant (IPA), has been successfully integrated into ExxonMobil’s proprietary reservoir simulator and deployed with it worldwide. The system can handle a large number of combinations of solver parameters, currently in the order of 108, and consistently improves run time performance of real simulation models, frequently by 30% or more, compared to the performance with the default solver settings. Moreover, IPA includes a persistent memory of solver performance statistics. The runtime statistics from these individual runs can be gathered, processed using data mining techniques and integrated in the IPA system, thus allowing its continuous improvement.

Geological media exhibit permeability fields and porosities that differ by several orders of magnitude across highly varying length scales. Computational methods used to model flow through such media should be capable of treating rough coefficients and grids. Further, the adherence of these methods to basic physical properties such as local mass balance and continuity of fluxes is of great importance. Both discontinuous Galerkin (DG) and mixed finite element (MFE) methods satisfy local mass balance and can accurately treat rough coefficients and grids. The appropriate choice of physical models and numerical methods can substantially reduce computational cost with no loss of accuracy. MFE is popular due to its accurate approximation of both pressure and flux but is limited to relatively structured grids. On the other hand, DG supports higher order local approximations, is robust and handles unstructured grids, but is very expensive because of the number of unknowns. To this end, we present DG-DG and DG-MFE domain decomposition couplings for slightly compressible single phase flow in porous media. Mortar finite elements are used to impose weak continuity of fluxes and pressures on the interfaces. The sub-domain grids can be non-matching and the mortar grid can be much coarser making this a multiscale method. The resulting nonlinear algebraic system is solved via a non-overlapping domain decomposition algorithm, which reduces the global problem to an interface problem for the pressures. Solutions of numerical experiments performed on simple test cases are first presented to validate the method. Then, additional results of some challenging problems in reservoir simulation are shown to motivate the future application of the theory.

Families of flux-continuous, locally conservative, finite-volume schemes have been developed for solving the general tensor pressure equation on structured and unstructured grids [1,2]. A family of flux-continuous schemes is quantified by a quadrature parameterization [3]. Improved convergence has been observed for certain quadrature points [4]. In this paper a q-family of flux-continuous (CVD-MPFA) schemes are used as a part of numerical upscaling procedure for upscaling fine scale permeability fields on to coarse grid scales. A series of data sets [5] are tested where the upscaled permeability tensor is computed on a sequence of grid levels. The upscaling sequence is repeated for three distinct quadrature points (q=0.1, q=0.5, q=1) belonging to the family of schemes. Three types of refinement study are presented: 1. Refinement study with invariant permeability distribution with respect to grid level, involving a classical mathematical convergence study. The same coarse scale underlying permeability map (and therefore the problem) is preserved on all grid levels including the fine scale reference solution. 2. Refinement study with renormalized permeability. In this study, the local permeability is upscaled to the next grid level hierarchically, so that permeability values are renormalized to each coarser level. Hence, showing the effect of locally upscaled permeability, compared to that obtained directly from the fine scale solution. 3. Reservoir field refinement study, in this study the permeability distribution for each grid level is obtained by upscaling directly from the fine scale permeability field as in standard simulation practice. The study is carried out for the discretization of the scheme in physical space. The benefit of using specific quadrature points is demonstrated for upscaling in the study and superconvergence is observed. REFERENCES 1. M.G. Edwards, C.F.Rogers. Finite volume discretization with imposed flux continuity for the general tensor pressure equation. Comput. Geo. (1998). 2. M.G.Edwards, Unstructured, control-volume distributed, full-tensor finite volume schemes with flow based grids. Comput.Geo. (2002). 3. M. Pal, M.G. Edwards and A.R. Lamb., Convergence study of a family of flux-continuous, finite-volume schemes for the general tensor pressure equation. IJNME, 2005. 4. M. Pal, Families of Control-Volume Distributed CVD (MPFA) Finite Volume Schemes for the Porous Medium Pressure Equation on Structured and Unstructured Grids. PhD Thesis, Swansea University, 2007. 5. M.A.Christie, SPE, Herriot-Watt University, and M.J Blunt, Imperial College, Tenth SPE comparative Solution Project: A Comparison of Upscaling Techniques, 2001.

We propose an upscaling method that is based on dynamic simulation of a given model in which the accuracy of the upscaled model is continuously monitored via indirect measures of the error. If the indirect error measures are larger than a specified tolerance, the upscaled model is dynamically updated with approximate fine-scale information that is reconstructed by a multi-scale finite volume method (Jenny et al., JCP 217: 627--641, 2006). Upscaling of multi-phase flow entails detailed flow information from the underlying fine scale. We apply adaptive prolongation and restriction operators for the flow and transport equations in constructing an approximate fine-scale solution. This new method reduces the inaccuracy associated with traditional upscaling methods, which rely on prescribed boundary conditions in computing the upscaled variables. This dynamic upscaling algorithm is validated for incompressible two-phase flow in two dimensional heterogeneous domains. We demonstrate that the dynamically upscaled model achieves high numerical efficiency compared with fine-scale computations and also provides excellent agreement with reference fine-scale solutions.

In reservoir simulations, we use an ensemble of fine-scale geological models that are upscaled to coarser representations for flow simulations. The primary objective of the upscaling process is to meet the computational limits of the simulator. However, we argue that from a system-theoretic point of view, a more fundamental underlying reason for upscaling is that the complexity level of a model should be adjusted to the amount of available information from measurements and the extent of control (input) exercised by adjusting the well parameters. That is because for a given configuration of wells, a large number of combinations of state variables (pressure and saturation values in the gridblocks) are not actually controllable and observable and accordingly they are not affecting the input-output behavior of the model. Therefore, we propose a “control-relevant-upscaling” (CRU) approach that determines equivalent coarse-scale permeabilities based on the actual system’s input-output behavior. The coarse-scale parameters are obtained as the solution of an optimization problem that minimizes the distance between the input-output behaviors of the fine- and coarse-scale models. This distance is measured by using Hankel- or energy norms, in which we use Hankel singular values and Markov parameters as a measure of the combined controllability and observability, and response of the system, respectively. This work focuses on single-phase flow upscaling, where we develop a CRU algorithm for reservoir systems. Moreover, we address the potential benefits of using proper orthogonal decomposition (POD) in combination with our CRU method to obtain a reduced-order CRU algorithm that accelerate the upscaling procedure. In the cases considered, the CRU algorithm shows superior input-output behavior as compared to upscaling algorithms commonly used in reservoir simulators.

In reservoir simulation the pressure is the solution of an elliptic equation. It follows from the maximum principle that this equation satisfies a monotonicity property. This property should be preserved when discretising the pressure equation. If this is not the case, the pressure solution may have false internal oscillations and extrema on no-flow boundaries. Thus, there is a need for sufficient conditions for the discretization methods to be monotone. These conditions will depend on the permeability tensor and the grid. Previously monotonicity for control volume methods has been studied on quadrilateral grids and on hexagonal grids. The former was done for general methods which reproduces linear potential fields, while the latter was done in a Control Volume Finite Element setting. These analyses have given sharp sufficient conditions for the discretisation methods to be monotone. However, for methods whose cell stencils include cells which do not have any edges common with the central cell in the discretisation scheme, no work has been done on unstructured grids. In this work, we study monotonicity on triangular grids for control volume methods which are exact for linear potential fields. We derive sufficient conditions for monotonicity of the MPFA-O and -L methods. The found monotonicity regions for the MPFA methods are also tested numerically. The tests are done both on uniform grids in homogeneous media, and on perturbed grids, which corresponds to heterogeneous media. The investigations are done for single phase flow only. However, the results are relevant for multiphase simulations. The results obtained in this work may be utilised in grid generation. In this way we can construct grids where the discretisation of the pressure equation is guaranteed to be monotone.

Today's geological models are increasingly complex. The use of unstructured grids is an efficient way to account for this complexity especially when dealing with large network of fractures and faults. Flux based approximations are commonly used in reservoir simulation community to discretize the flow equations. The two-point flux approximation (TPFA) is the simplest and the most robust discretization technique. Application of TPFA requires an orthogonal grid. When applied to an unstructured model, Perpendicular Bisector (PEBI) grids are usually used as they are orthogonal by construction. But the construction of PEBI grids can become challenging for fully 3D models containing a large number of fractures and faults. In this work we propose to use an unstructured grid which matches exactly all discrete features (fractures and faults) but may not be orthogonal. The objective is to relax the orthogonality criteria to simplify the grid generation step and to deal with the orthogonality at the flux approximation step. Typical discretization techniques for nonorthogonal grids are based on multiple-point flux approximation (MPFA) where more than two pressure points are used to evaluate the flux. Although substantial progress has been made in improving MPFA, they remain more complex and less robust than TPFA. In this work we present a simple methodology to apply a TPFA to unstructured nonorthogonal grids. The idea is to determine the location of the pressure node inside each control-volume to make the connections with the surrounding cells as orthogonal as possible. Although it is not always possible to make all connections perfectly orthogonal, this optimization procedure improves systematically the grid quality. In addition, this purely geometrical optimization is done at the preprocessing level and does not require any flow information and can be used with any connectivity based flow solver.

In reservoir simulation, a proper treatment of wells is crucial for the performance of fluid flow models. Important well parameters such as the well flow rate and the wellbore pressure are highly sensitive to the computational accuracy of near-well flow models. Near-well regions are characterised as high flow density regions, and the dominating flow pattern exhibits a radial-like nature with large pressure gradients and, for multiphase flow, large saturation gradients. Skew or horizontal wells will also in general imply a strong effect of anisotropy and heterogeneities due to geological layering of the reservoirs. Due to difference in the reservoir scale and the wellbore radius, reservoir flow models do not capture the true flow behavior in the well vicinity. Thus, more flexible models with local grid refinement in near-well regions are needed. The models should also be adapted to handle complex geological near-well structures and multiphase flow simulations. Existing near-well models are in general based on homogeneous media in the well vicinity. Anisotropies and strong heterogeneities are less accounted for in near-well flow simulations. In this work, we construct analytical solutions for near-well flow which is not aligned with a radial inflow pattern. These solutions resemble strongly heterogeneous, possible anisotropic media. We compare different control volume discretisation schemes and radial-type grids for such cases and give their convergence behavior. The objective is to obtain a clearer view on the accuracy of the near-well grids and discretisation schemes for large contrast in permeability, and hence, to determine which is the preferable grid and a suitable numerical scheme given certain near-well conditions. The simulations are performed for singlephase flow, however, the results will also have relevance for multiphase flow simulations.

Operator splitting time discretization techniques are getting more and more popular for the numerical simulation of non-linear problems in porous media. The main idea is to split complex operators in evolution equations into simpler ones which are successively solved in each time step. For porous media flow a natural splitting is given by the diffusive and the advective part of the flux. Then for each part, optimal solvers for the particular evolution equation can be used, i.e., a parabolic solver for the diffusion equation and a hyperbolic solver (e.g., characteristics methods) for the advection equation. In this talk, we first present an operator splitting discretization using an implicit finite volume scheme for the diffusive part and a semi-Lagrangian method for the advective part. Hence, the computational cost of the semi-Lagrangian method can be neglected compared to the one of the implicit finite volume scheme. In the second part, we introduce a two-scale operator splitting method where the diffusion part is solved on a coarse grid while the advection part is considered on a fine grid. This is motivated by the fact that the diffusion and the advection act on different scales. Then, the resulting operator splitting discretization is equivalent in performance to a fast hyberbolic solver. In numerical examples, our proposed methods are compared with standard implicit finite volume methods. We will demonstrate that by using the operator splitting time discretization, the number of iterations in the implicit finite volume scheme may be significantly reduced.Moreover, for the two-scale method, it turns out that despite of the poor approximation in the diffusive part, the solution is still of equal quality compared to that of standard methods on fine grids. To summarize, operator splitting techniques provide a very flexible and powerful framework which makes it quite attractive for an efficient solution strategy for flow problems in porous media.

The typical scales of pore pressure and fluid saturations variations are different both for space and temporal variables. The multiscale method is based on using the coarse grid for pore pressure equation and fine grid for saturation equations. The time step for these equations may be different too. The essential feature of the method is the basic function construction by solving the one phase stationary equations. These basis functions have the same peculiarities determined by the fine grid structures. The solution of pore pressure equations may be constructed as linear span of these basic functions for Buckley-Leverett equations. In general case it is possible to construct the tensor total permeability coefficients and to up-scale the relative permeability. The method of calculating these variables is the dissipative energy integral approximation. As a result we have the non-linear parabolic equation for pore pressure with tensor coefficients. The implicit schemes on the fine grid are considered for saturation equation. The general method is high resolution method. It is essential this method is a high resolution one. The results of modeling some problems are represented.

If oxygen is in intimate contact with fuel Arrhenius law says that reaction will occur, even at low temperatures. Heat losses can become equal to the small reaction heat generated, so that the system remains trapped in a slow reaction mode. Such a mode is indistinguishable from extinction. On the other hand, if heat losses remain smaller than the heat generated by the reaction the temperature increases and spontaneous ignition occurs. Heat losses are strongly dependent on the geometry of the heat generating region. We therefore distinguish three idealized geometries, viz. linear, cylindrical and spherical infinite domains. We analyze the long time behavior of a basic heat diffusion problem that incorporates an Arrhenius heat generation term; the resulting model is temperature controlled. Only in spherical geometry the model recovers the results from the conventional heat transfer concept. Indeed, we come to the following results. As reactant depletion is ignored, spontaneous ignition always occurs, in linear and cylindrical coordinates, even if it takes very long time. For the spherical case, the occurrence of ignition or extinction depends on the process conditions.

This paper presents an efficient numerical scheme for multiphase compositional flow coupled with subsurface heat transfer. The flow equations are first presented followed by a brief discussion of the equation of state (EOS) and a description of the two-phase flash algorithm. An implicit-pressure, explicit-concentrations (IMPEC) sequential algorithm is then applied iteratively to enforce the non-linear volume balance saturation) constraint. The pressure equation is solved using mixed FEM, while the concentrations are updated consistent with requiring local mass balance of every component. Thermal effects also play an important role in such problems since they effect the phase properties and hence, stable and accurate locally conservative methods are desirable to model the thermal energy balance equation. To this end, we present also a time-split scheme for modeling the energy balance equation that is sequentially coupled to the flow solution rendering it cheap yet accurate for the complex problems being modeled. Results of benchmark problems in compositional flow modeling are presented and validated where possible. Some large-scale parallel tests are performed on more challenging applications such as those with highly heterogeneous permeability fields on very fine grids. Efficient parallel scalability of the code on upto 512 processors is also demonstrated. Finally, some test cases simulating "real world" problems of CO2 sequestration in deep saline aquifers are presented. Initial results show reasonably good agreement of CO2 plume shapes, arrival times, leakage rates and production curves with semi-analytic solutions, wherever available. All computations were carried out within the IPARS-CO2 code base framework.

High-resolution geological models of subsurface reservoirs typically contain much more details than can be used in conventional reservoir simulators. Geomodels are therefore upscaled before flow simulation is performed. Here, we present a non-uniform coarsening strategy to upscale geomodels, where the coarse grid is generated by grouping cells in the fine-scaled geomodel selectively into connected coarse blocks, with some minimum volume and with the total flow through each coarse block bounded above. Transport is then modeled on this simulation grid. For purely advective flow, the coarsening strategy has been shown to be robust, allowing for both accurate and fast transport simulations of highly heterogeneous and fractured reservoirs. Here, we consider advection-dominated two-phase flow, where the capillary diffusion is discretized separately from the advective and gravity terms using operator splitting. In particular, we investigate different damping strategies of the diffusive term to counteract overestimation of the diffusion operator on the coarse grid, this to ensure accurate diffusion modelling in the transport solver.

Numerical solution of conservation laws with large variations in the local CFL number, requires impractically small time steps for stable, explicit time integration. Implicit time stepping with large time steps is therefore necessary but is computationally expensive and consequently limited to first-order discretization schemes. A numerical scheme that is accurate and at the same time allows large time steps is clearly desirable. We develop a numerical formulation based on the Adaptive Implicit Method (AIM) that achieves both objectives by treating implicitly only a few grid cells that restrict the time step size, and by using high-order spatial discretization in both the implicit and explicit cells. Our novel approach solves the fully coupled transport and flow equations with high order-accuracy within the AIM framework to reduce the computational cost and improve the solution accuracy. We show that high-order discretization is computationally efficient when combined with AIM and leads to significant improvements in accuracy compared to the first-order solution.

Geostatistical simulation of permeabilities on a geologically detailed resolution will account for permeabilities in the cells containing wells through kriging. These permeabilities are typically based on porosity logs and core plug measurements of both porosity and permeability. No direct measurement of permeability on the geomodelling grid scale exists. However, well tests give information on the effective permeability in an area close to the well region, covering several grid cells, and are therefore data on an aggregate scale. Given an assumption of radial flow into a vertical well, the effective permeability becomes a convolution of the permeabilities in the well test region. Downscaling of the well test data is possible by co-kriging the aggregate scaled permeability field with that of the geomodelling grid scale. Thereby, the geomodelling grid permeabilities will honour data on both scales. The effective permeability is downscaled through inverse block kriging, which implies a deconvolution procedure. Keeping the computation costs low when introducing a new conditioning parameter is ensured by transforming into the Fourier domain, since the Fast Fourier Transform algorithm is an efficient method for solving the inverse block co-kriging. Initial testing of the implemented algorithm has been made on real data from a StatoilHydro operated field on the Norwegian continental shelf in a proof of concept test. Three wells with well test data were chosen, and their derived effective permeabilities were included in the permeability simulations. Cases both with and without well test conditioning were run, and compared in a well test simulation software. All three near well regions showed a significant improvement. These tests indicate that the conditioning method can be a useful contribution in bringing dynamic data into the reservoir characterization.

This work presents new Kriging-based algorithms that allow incorporating nonlinear constraints such as a given average and variance into the kriged field. The first of these constraints allows one to obtain, for example, permeability fields whose average matches well-test-derived permeability; the second one enables one to generate fields that have a desired variability. Therefore, a well-known drawback of Kriging, namely excessive smoothness, is overcome. As these constraints are applied, the Kriging estimates of all blocks become mutually dependent and must be solved simultaneously. This fact led to the development of a new concept, that of Simultaneous Kriging. Furthermore, because permeability is not an additive variable, the problem becomes nonlinear; an optimization procedure based on Newton's method is used to solve it. Both average and variance constraints are applied through Lagrange multiplier technique. Power averaging with a ω exponent that varies from -1 to 1 is used to estimate field-generated permeability within the well-test drainage area. The optimization algorithm searches for an optimum ω value as well as a set of gridblock values such that the desired field average and variance are met. Starting from an initial guess and tolerance, the algorithm will reach the closest minimum within the solution space. Unlike other Kriging methods, the solution is not unique; in fact, infinite, equiprobable solutions can be found. From this point of view, the proposed algorithms become more like conditional simulation. A number of permeability fields were generated and used as input of a numerical simulator to calculate their actual well-test permeability. The simulation results proved that the fields generated by the algorithm actually matched well-test permeability and field data variance within a reasonable tolerance. The same procedure can be used to incorporate other nonlinear constraints, such as facies proportion.

The Ensemble Kalman Filter (EnKF) is derived under the assumption of Gaussian probability distributions. Thus, the successful application of the EnKF when conditioning to dynamic data depends on how well the stochastic reservoir state can be approximated by a Gaussian distribution. For facies models, the distribution becomes non-Gaussian and multi-modal, and the EnKF fails to preserve the qualitative geological structure. To apply the EnKF with models where the Gaussian approximation fails, we propose to map the ensemble to a higher dimensional feature space before the analysis step. By careful selection of the mapping, moments of arbitrary order in the canonical reservoir parameterization can be embedded in the 1st and 2nd order moments in the feature space parametrization. As a result, the parameterization in the feature space might be better approximated by a Gaussian distribution. However, the mapping from the canonical parameterization to the feature space parameterization does not have an inverse. Thus, finding the analyzed ensemble of canonical states from the analyzed ensemble represented in the feature space is an inverse problem. By using the kernel trick, the mapping to the feature space is never explicitly computed, and we solve the inverse problem efficiently by minimizing a cost function based on the feature space distance. As a result, the computational efficiency of the EnKF is retained, and the methodology is applicable for large scale reservoir models. The proposed methodology is evaluated as an alternative to other approaches for estimation of facies with the EnKF, such as the truncated pluri-Gaussian approach. Results from a field case are shown.

Over the last decade, the ensemble Kalman filter (EnKF) method has become an attractive tool for history matching reservoir simulation models and production forecasting. Recently, EnKF has been successfully applied to real field studies (e.g. Bianco et al., 2007). Nonetheless, Lorentzen et al. (2005) observed a consistency problem. They showed, using the Kolmogorov-Smirnov test, that for the PUNQ-S3 model, the posterior cdfs of the total cumulative oil production for 10 ensembles were not coming from the same distribution, as was expected since the 10 initial ensembles were generated using a common distribution. The forecasts from these initial ensembles gave consistent cdfs. We investigate if this issue is related to the inbreeding of the ensemble members when applying EnKF. Houtekamer and Mitchell (1998) proposed to use a paired EnKF with covariance localization to improve atmospheric data assimilation. This method was developed to hinder ensemble collapse due to inbreeding of ensemble members and spurious correlations problems far from the observation points. We show that using a paired EnKF, where each Kalman gain is used to update each other's ensemble, do not in fact prevent inbreeding. A new approach, coupled EnKF, that do prevent inbreeding is presented. The present work first review the issue of consistency encountered with EnKF and investigate the use of paired and coupled EnKFs without covariance localization as a possible remedy to this problem. The method is tested on a simple nonlinear model for which the posterior distribution can be analytically solved. The obtained results are compared with the traditional EnKF and the analytic solution. Next, a hierarchical ensemble filter approach inspired by Anderson (2007) acting as a covariance localization technique is proposed and tested on the PUNS-S3 model against the traditional EnKF. This approach seems better than paired or coupled EnKF to help solve the observed inconsistencies.

Uncertainty quantification is probably one of the most important and challenging aspects of reservoir modeling, engineering and management. Traditional approaches have failed to provide manageable solution models. In most cases, uncertainty is represented through either a multi-variate distribution model (posterior) in a Bayesian modeling context or through a set of multiple realizations in a more frequentist view. The traditional Monte-Carlo simulation approach, whereby multiple alternative high-resolution models are generated and used as input to reservoir simulation or optimization codes, has not been proven to be practical mainly due to CPU limitation, nor is it effective or flexible enough in terms of addressing a varying degree of uncertainty assessment issues. In this paper, we propose a reformulation of reservoir uncertainty in terms of distances between any two reservoir model realizations. Instead of considering a reservoir model realization as some point in a high-dimensional space on which a complex posterior distribution is defined, we consider a distance matrix which contains the distances between any two model realizations. The latter matrix is of size NR×NR, NR being the number of realizations, much less for example that an N×N covariance matrix, N being the number of grid-blocks. Note that, unlike a co-variance matrix or probability function, a distance can be tailored to the specific problem at hand, for example water-breakthrough, or OOIP, even though the realizations remain the same. Next, the classical Karhunen-Loeve expansion of a Gaussian random field, based on the eigenvalue decomposition of an N×N covariance table, can now be formulated as function of the NR×NR distance matrix. To achieve this, we construct an NR×NR kernel matrix using the classical radial basis function, which is function of the distance and perform eigenvalue decomposition of this kernel matrix. In kernel space, new realizations can be generated or adjusted to new data. As application we discuss a new technique, termed metric EnKF yo simulataneously update multiple non-Gaussian realizations with production data. Using this simple dual framework to probabilistic-based uncertainty, I show how many types of reservoir uncertainty (spatial, non-spatial, scenario-based etc..) can be modeled in this fashion. I show various applications of this framework to the problem of history matching, model updating and optimization.

This paper presents non-intrusive, efficient stochastic approaches for predicting uncertainties associated with petroleum reservoir simulations. The Monte Carlo simulation method, which is the most common and straightforward approach for uncertainty quantification in the industry, requires to perform a large number of reservoir simulations and is thus computationally expensive especially for large-scale problems. We propose an efficient and accurate alternative through the collocation-based stochastic approaches. The reservoirs are considered to exhibit randomly heterogeneous flow properties. The underlying random permeability field can be represented by the Karhunen-Loeve expansion (or principal component analysis), which reduces the dimensionality of random space. Two different collocation-based methods are introduced to propagate uncertainty of the reservoir response. The first one is the probabilistic collocation method that deals with the random reservoir responses by employing the orthogonal polynomial functions as the bases of the random space and utilizing the collocation technique in the random space. The second one is the sparse grid collocation method that is based on the multi-dimensional interpolation and high-dimensional quadrature techniques. They are non-intrusive in that the resulting equations have the exactly the same form as the original equations and can thus be solved with existing reservoir simulators. These methods are efficient since only a small number of simulations are required and the statistical moments and probability density functions of the quantities of interest in the oil reservoirs can be accurately estimated. The proposed approaches are demonstrated with a 3D reservoir model originating from the 9th SPE comparative project. The accuracy, efficiency, and compatibility are compared against Monte Carlo simulations. This study reveals that, compared to traditional Monte Carlo simulations, the collocation-based stochastic approaches can accurately quantify uncertainty in petroleum reservoirs and greatly reduce the computational cost.

First, we present an efficient method for integrating dynamic data in high-resolution subsurface models. The method consists of two key technologies: (i) a very fast multiscale-streamline flow simulator, and (ii) a fast and robust 'generalized travel-time inversion' method. The travel-time inversion is based on sensitivities computed analytically along streamlines using only one forward simulation. The sensitivities are also used to selectively reduce the updating of basis functions in the multiscale mixed finite-element pressure solver. Second, we discuss extensions of the methodology to grids with large differences in cell sizes and unstructured connections. To this end, we suggest to use rescaled sensitivities (average cell volume multiplied by local sensitivity density) in the inversion and propose a generalized smoothing operator for the regularization to impose smooth modification on reservoir parameters. Two numerical examples demonstrate that this reduces undesired grid effects. Finally, we show a slightly more complex example with two faults and infill drilling.

Variational data assimilation techniques (automatic history matching) can be used to adapt a prior permeability field in a reservoir model using production data. Classical variational data assimilation requires, however, the implementation of an adjoint model, which is an enormous programming effort. Moreover, it requires the results of one complete simulation of forward and adjoint models to be stored, which is a serious problem in real-life applications. Therefore, we propose a new approach to variational data assimilation that is based on model reduction, where the adjoint of the tangent linear approximation of the original model is replaced by the adjoint of a linear reduced model. The Proper Orthogonal Decomposition approach is used to determine a reduced model. Using the reduced adjoint the gradient of the objective function is approximated and the minimization problem is solved in the reduced space. If necessary, the procedure is iterated with the updated estimate of the parameters. We evaluated the model-reduced method for a simple 2D reservoir model. We compared the method with variational data assimilation where the gradient is approximated by finite differences and we found that the reduced-order method is about 50 % more efficient. We foresee that the computational efficiency will significantly increase for larger model size and our current research is focused on quantifying this computational benefit.

Uncertainty in reservoir models can be quantified by generating large numbers of history-matched models, and using those models to forecast ranges of hydrocarbons produced. The need to run large numbers of simulations inevitably drives the engineer to compromises in either the physics represented in the reservoir model, or in the resolution of the simulations run. These compromises will often introduce biases in the simulations, and the unknown reservoir parameters are estimated using the biased simulations, which can lead to biases in the parameter estimates. Solution error models can be used to correct for the effects of the biases. Solution error models work by building a statistical model for the differences between fine and coarse simulations (or between full physics and reduced physics simulations) using data from simulations at a limited number of locations in parameter space. The statistical model then produces estimates of the error elsewhere in parameter space; these estimates are used to correct the effects of the coarse model biases. In this work, we apply a solution error model to material balance calculations. Material balance is frequently used in reservoir engineering to estimate the initial oil in place. However such models are very simple, treating the reservoir as a tank and allowing instantaneous equilibration of fluids within the tank. The results of material balance simulations will therefore not be consistent with multi-cell reservoir simulations. We use a model based on Teal South Reservoir in the Gulf of Mexico to demonstrate how an error model can correct a material balance model to the accuracy of a reservoir simulation.

Popularity of application of the ensemble Kalman filter (EnKF) and other ensemble methods to history matching of reservoir models has been growing for a number of years. Recent applications demonstrated that ensemble based methods are efficient for producing models with a good match of the history. However, they don’t necessarily keep consistency of the resulting models with the initial geostatistical description. One of the causes of this problem is an imperfection of the estimation of the model state error covariance from ensembles of a finite size. First we illustrate a spurious stochastic behaviour of the covariance when applying the EnKF to a simple "well-test like" synthetic example. Experimental covariance estimated from the ensemble of a practical size of about a hundred members is clearly quite different from one that can be expected from an analytical solution. When the size of the ensemble is increased by an order of magnitude the covariance function becomes significantly smoother. Nevertheless, in all the cases it is not difficult to see that the main disturbance of the covariance happens beyond the area of the pressure front propagation. Thus, results of the EnKF potentially can be improved by application of a flow based covariance localisation. Then we show an application of the EnKF with a streamline based covariance localisation (Devegowda et al. 2007) to a real field problem. The forward model is solved by a conventional finite-difference reservoir simulator and at every update step streamlines are traced using a separate routine (Jimenez et al. 2007). Traced streamlines are used for determination of an influence zone for each well and covariance of the production data from each well and model parameters is localised accordingly. We compare the results in case of localisation based on the zones determined by the full length of the streamlines and in case of localisation zones limited by the flow fronts propagation. In both cases there is an improvement of the production data match and less disturbance of the initial geostatistical realisations compared to the EnKF without localisation.

The context of this work is a dynamic data assimilation workflow where a relatively small number of simulator inputs (usually around 10) has been selected in order to calibrate the reservoir simulation model. The probabilistic approach is used to solve this complex ill-posed inverse problem where the objective is to obtain a posterior distribution of the selected parameters. The remaining uncertainty on these simulator inputs is then propagated on the output of interest such as oil, water and gas productions to obtain confidence intervals for future forecasts. To reduce the number of required simulations the objective function is approximated using a non parametric response surface method based on Gaussian Process regression (kriging). To obtain a predictive response surface a sequential design of experiments is adopted that aims at discovering the minima of the objective function and also to accurately reproduce the basins associated to these minima. At each step of the design an MCMC method is used to explore the minima and to select the new simulations to perform. Also, differently than in previous works the response surface estimated variance is included in the posterior computation. As a result, the response surface accuracy improvement obtained by the sequential design produces also a reduction of uncertainty in the obtained estimates of the input posterior distribution. This uncertainty reduction effect together with the accuracy improvement of the response surface are monitored at run time at each iteration of the sequential design. An application is presented on the PUNQS field case with seven uncertain parameters. Less than two hundreds fluid flow simulations were used to produce a reliable sample of the posterior distribution and to propagate the remaining uncertainty on future forecasts.

Significant progress has been made recently in the numerical simulation of heterogeneous reservoir media. One of the fundamental reasons for the hysteretic nonlinear behavior of porous reservoir media is that heterogeneous or damaged materials contain an enormous number of mesoscopic features such as microcracks and macrocracks, joints, and grain to grain contacts containing multiple phases. Each of these mesoscopic units exhibits a hysteretic behavior which dominates the macroscopic reservoir response. Based on the work of Preisach and Mayergoyz (P-M space model), Coleman and Hodgdon, a generalized phenomenological model has been developed to describe the hysteretic nonlinear response of capillary pressure and relative permeability. The proposed approach enables the description of active hysteresis by the solution of either a differential or an integral equation. The model focuses on the correct representation of the primary drainage (forced), the imbibition (spontaneous and forced), the secondary drainage (spontaneous and forced) curves and scanning curves. The functions and parameters used in the model can be fine-tuned to match different experimental data or can be used as history matching parameters. It is shown that the proposed model incorporates the Killough type hysteresis as one analytical solution. The differential form of the proposed model allows a smooth transition of both relative permeability and capillary pressures from drainage-to-imbibition or imbibition-to-drainage states and requires minimum storage of parameters during the simulation. Validation of the model indicates that the proposed hysteresis model is stable and robust. The results are presented and discussed and future studies outlined.

Simulation of multi-phase flow processes for enhanced oil recovery requires accurate descriptions of the capillary pressures and relative permeabilities as functions of the phase saturations Pore-scale network modelling is useful tool to estimate these flow parameters. Some of the main network model characteristics are the pore and throat conductances and capillary entry pressures. Both parameters strongly depend on the pore and throat geometries. Because the shapes of the real pore cross-sections are generally highly irregular, it is important to use idealized shapes that lead to accurate approximations of the above parameters. The most common approach has been to choose a circle/ (irregular) triangle / square (C-T-S) pore geometry with a shape factor that matches that of the real pore shape. For these shapes, simple correlations between the flow parameters and the shape factor are available. However, it is well known that the parameters for these very regular convex shapes are often inaccurate compared to the real pore shapes. Here, we propose to represent the shapes by the regular, but generally non-convex, n-cornered star shapes. A new n-corner star shape characterization technique has been developed, which takes shape factor and dimensionless hydraulic radius as input parameters. A novel numerical technique has been used to derive the real pore entry pressures, whereas analytical expressions are used for the star shape. A set of 70 individual pores has been extracted from high-resolution 2D images of a Bentheim rock sample. The real shapes have been approximated by C-T-S, as well as by n-corner star shapes. The comparison results between predicted and real shape parameters for the entire set of pores show that the accuracy of the entry radius prediction for both approaches is approximately the same and quite good. The single phase conductance estimation is much better for n-corner star approximation than for C-T-S. Finally, a capillary bundle model has been constructed from the 70 pores to test the predictive capabilities of the characterization approaches in terms of relative permeabilies and capillary pressures. It has been shown that the n-corner star provides a better approximation of the "real shape" capillary pressure curve than C-T-S. The real shape relative permeabilities are in a good agreement with curves predicted by both shape approximations.

Faults in reservoir simulation models are typically represented as discontinuities between grid cells with homogeneous transmissibility multipliers used to represent their effect on fluid flow. Transmissibility multipliers may be manually altered in a trial-and-error fashion to condition the model to production data, often at the expense of geological and physical realism. Methods developed in the past decade incorporate prior geological constraints on fault permeability and thickness to produce more realistic heterogeneous transmissibility multipliers. As a typical simulation grid may have several thousand faulted cell connections, conditioning of these properties to production data is not routine. We present a method to condition geologically-derived estimates of fault transmissibility to production data. Appropriate probability density functions are estimated from databases of fault permeability and thickness. Together with parameters for grid properties, these are used as input to a geostatistical property simulation workflow in a 3D geomodelling tool. Multiple property realisations are generated, with fault permeability being averaged into grid permeability and exported for use in numerical fluid flow simulations. Conditioning to production data is achieved using the Ensemble Kalman Filter (EnKF), which has been proved as an efficient approach for such problems. Since the EnKF can handle a large number of uncertain parameters and requires the forward model only as a ‘black box’, it allows for consideration of geological features, such as faults, as objects with spatially continuous heterogeneous properties that can be conditioned to production data. The method is illustrated using a producing North Sea reservoir. Grid and fault properties can be conditioned to production data while honouring the input probability distributions and spatial continuity model. Case-specific and more general conclusions can be drawn about the advantages of the method; potential improvements and extensions are discussed.

Permeability samples for performance forecasting should be consistent with all available information. Markov chain Monte Carlo (McMC) samples correctly from the posterior distribution, but is computationally extremely intensive. RML is a good approximation to McMC. RML involves history matching of each permeability sample, and although RML is computationally less intensive than McMC, the computational effort to generate a sufficiently large number of samples is huge for all but very small reservoir models. In this work, we have investigated if the sampling procedure can be made more efficient by using a predictor-corrector approach in the history matching step of RML. The predictor applies sequential parameter estimation to obtain an estimate with few degrees of freedom, utilizing only part of the available information. The corrector downscales the predictor estimate in a two-step procedure involving all available information, including the estimate obtained with the predictor. The first corrector step is a variant of Kriging. The second corrector step is parameter estimation, again involving few degrees of freedom, with basis functions derived from the results of the predictor.

Recent years the fractured reservoir characterization is of special interest since the dual medium simulation has wide application in field studies. In practice the main difficulties that impede the simulation are evaluation of fracture properties and including of their continuous analog into reservoir model. The research focuses on characterization of deformable fractured reservoir. A new approach to evaluate fracture properties from well test has been developed. Theoretical basis of the approach is compressible discrete fracture network (CDFN) model which is an extension of the DFN concept to take into account compressibility of fractures in terms of both porosity and permeability. The fracture compressibility is related to matrix (rock) compressibility that generates change in matrix block size and fracture aperture during pressure evolution. The new approach consists of two components. The first is numerical analysis of pressure drawdown and build-up to discover behaviour that is appropriate to the compressible fracture network. The second is well test interpretation to characterize fracture density; aperture, permeability, porosity and their dependence on pressure. The interpretation consists in solution of an inverse problem using analytical procedures and can be applied to identify properties of purely fractured and fractured porous reservoirs. Simultaneously the numerical simulation of flow in CDFN is considered. The simulation of fluid flow in near well bore area allows to calibrate properties of CDFN. The approach has been tested per analysis of an actual transient pressure well test. The flow simulation allowed to compare actual and synthetic transient pressure response and hence to estimate correctness and accuracy of the approach. The developed model, approach and flow simulation allow to identify fracture properties and to calibrate them using well test simulation. Obtained fracture permeability and porosity, their dependence on pressure and matrix block size are used in dual (single) medium reservoir simulation.

A new approach for early stage reserves estimation for dual porosity fractured reservoir is presented. When assessing the volumetrical reservoir content by Monte-Carlo simulation, rock volume, porosity and saturations are treated as stochastic variables. When estimating the reserves the recovery factor must be also handled as stochastic variable. The recovery factor for a given depletion mechanism depends on already mentioned parameters and additionally on the shape factor, the matrix permeability and the wettability. This dependency will be assessed by means of fine scaled numerical simulation resulting in so-called recovery curves. This paper discusses how the recovery curves will be created and how they will be used for reserve estimation. The recovery curves are generated using single porosity small scale models and are in terms of recovery versus time or dimensionless time. The block size of this small scale models corresponds to a given shape factor in the shape factor distribution. The matrix depletion processes are basically dependent on time and matrix parameters. Therefore, matrix parameters that have an influence on the depletion are incorporated in the dimensionless time. These are the shape factor, the matrix permeability, the oil viscosity and a characteristic pressure in the system. The paper also contains the application of the method to a field case.

Local parameterization techniques, such as the gradual deformation method and the perturbation method, have already proved to be efficient for the integration of production data in a history matching process, considering regions of the reservoir associated to wells and defined according to physical or geometrical criteria. The aim of this study is to apply local parameterization techniques for the integration of 4D seismic data, using a new definition of regions based on the error on the seismic data. Each time seismic data are available, streamlines arriving in the parts of the reservoir with the highest error on these data are computed. The aim is to identify “influence areas” which should contribute to the behavior in the badly matched regions. The successive partitions can be used sequentially. Several methods are also proposed to combine them in a single partition. We intend in this study to improve the matching of the saturation distribution in the reservoir obtained at different times from seismic data interpretation. Considering for instance a constant porosity, an increase (decrease) of the permeability in the influence areas associated to a delay (advance) of the saturation front may improve the match. This can be achieved by modifying locally the mean of the permeability, using for instance a non-stationary mean to simulate the realization. We also propose to modify linear volume averages of the simulated petrophysical properties distributions, in one or several regions, using a kriging based methodology. Using this methodology, history matching processes of production and saturation data have been performed, optimizing local averages of the permeability and porosity realizations in the influence areas. The results show a strong improvement of the saturation match. They are also compared to the ones obtained with local gradual deformation based optimizations considering the same partitions of the reservoir.

We used the representer method to perform data assimilation (automatic history matching) of production data and 4D seismic data to estimate the permeability field in numerical reservoir models. The 4D seismic data were incorporated in the form of interpreted time-lapse density changes. The representer method requires two numerical simulations per measurement which becomes infeasible if seismic data in all grid blocks would be taken into account. Therefore, only the most informative data at the saturation front moving over time were used in the assimilation. The method was tested in a twin-experiment using synthetic data. The results were assessed in terms of the quality of the history match (mismatch between ‘true’ and simulated measurements) and in terms of predictions (mismatch between ‘true’ and simulated predictions after the history matching period). We found a clear improvement of the joint assimilation inversion over individual assimilation of production and seismic data.

We present a general framework for thermodynamic equilibrium calculations of multi-phase, multi-component mixtures. We use the fact that the compositional space can be represented as a high-dimensional simplex. For given values of pressure and temperature, the phase behavior of a particular system can be described using a low-dimensional tie-simplex, for example, tie-lines and tie-triangles for two and three phase systems, respectively. In general, a high-dimensional compositional space can be parameterized using a lower dimensional tie-simplex sub-space. Tie-simplex computation and interpolation procedures complement the parameterization to complete the mathematical framework. The robustness and efficiency of the method is demonstrated using several multi-phase equilibrium problems of practical interest. One type of problem is the equilibrium flash calculation of systems with a large number of phases. The complexity and strong nonlinear behaviors associated with such systems pose serious difficulties for standard techniques. The tie-simplex representation of the equilibrium data, which may be obtained using a particular equation-of-state for example, can be pre-processed. The parameterized space can be used to obtain the phase compositions (i.e., flash results) or used as an initial guess to accelerate convergence of standard Equation of State (EoS) based procedures. The second type of application is multi-phase multi-component displacement problems. In the standard compositional simulation approach, an EoS is used to describe the phase behavior. For each gridblock, given the temperature, pressure and overall compositions, the EoS is used to detect the phase state (e.g., one, two, or more phases), and to calculate the phase compositions, if multiple phases are present. These EoS computations can dominate the overall simulation cost. Adaptive computation of the tie-simplex space can be used to speed up the EoS computations of large-scale problems of practical interest by an order of magnitude, or more.

Integrated inversion of production history data and 4D seismic data for reservoir model characterization leads to a nonlinear inverse problem that is usually cumbersome to solve : the associated forward problem based, on one hand, on fluid flow simulation in the reservoir for production data modeling, and on the other hand, on a petro-elastic model for 4D time lapse seismic data modeling, is usually computationally time consuming, the number of measurements to be inverted is large (up to 500 000), the number of model parameters to be determined is up to 100. Moreover, all the derivatives of the modeled data with respect to those parameters are usually not available. We propose an optimization method based on a Sequential Quadratic Programming algorithm which uses gradient approximation coupled with a BFGS approximation of the Hessian. In addition, the proposed method allows to handle equality and inequality nonlinear constraints. Some realistic applications are presented to illustrate the efficiency of the method.

ulation are increasingly included in best-practice workflows in the Oil & Gas industry. Most optimization methods applied to model validation in reservoir simulation, including so-called Evolutionary Algorithms like Genetic Algorithms (GA) and Evolution Strategies (ES), use an objective function definition based on the overall simulation period. The integration of a sequential data assimilation process is conceptually not embedded in those optimization methods. The Ensemble Kalman Filter (EnKF) has entered this field for its appealing features. Sequential data assimilation allows the implementation of real-time model updates where classical optimization techniques require simulating the complete history period. This may have negative effects on efficiency and use of computing time. In contrast, EnKF sequentially assimilates information streams into a set of numerical models. While being a special case of a fully fledged particle filter the EnKF method with application to reservoir simulation has proven to generate results with a reasonable amount of ensemble members. The similarities between a particle filter (Monte Carlo Filter) and an Evolutionary Algorithm (Generic Algorithm) have been previously pointed out from a rather theoretical point of view (1,2). In this work we present a concise overview of Evolutionary Algorithms and the Ensemble Kalman Filter in such a way that the cross-relations become apparent. Similarities are highlighted and the potential for hybrid couplings is discussed. Practical implications for the implementation of these methods are derived. 1. Higuchi, Tomoyuki. Monte carlo filter using the genetic algorithm operators. Journal of Statistical Computation and Simulation. 1997, Vol. 1, 59, pp. 1-23. 2. —. Self-organizing Time Series Model. [book auth.] A. Doucet, J.F.G. de Freitas and N.J. Gordon. Sequential Monte Carlo Methods in Practice. s.l. : Springer, 2001, pp. 429-444.

The determination of the main petrophysical parameters is normally based on the log data. During the log interpretation process suitable mathematical algorithms taking into account all the available data and information are applied in the attempt to reliably relate log measurements to mineralogy, porosity and water saturation. Uncertainty is always associated to the results of the interpretation due to possible errors in the measurements, in the selected petrophysical models, and/or in the input parameters required by the models. The possibility to estimate the reliability of the petrophysical characterization of the reservoir rocks has a strong impact on the evaluation of the hydrocarbon originally in place (HOIP) and, thus, on the technical and economical exploitation strategies. The evaluation of the uncertainties associated to the results of the well log interpretation process can be performed only by applying a methodology that couples a robust optimization process to a representative statistical approach. On the basis of previous studies and applications to real cases, a methodology for log inversion and uncertainty estimation was formulated. According to this methodology, log interpretation was performed using the iterative solution of the Lagrangian relaxed problem with the Gauss-Newton algorithm, in which the constrains were managed with the active set method; the Monte Carlo statistical approach was applied to the log interpretation routine in order to assess the associated uncertainty. The use of a fast iterative inversion method proved fully compatible with the use of the Monte Carlo approach to estimate the range of uncertainty associated to the reservoir characterization. The rigorous formulation of the methodology and a discussion of the applicability limits and convergence requirements of the inversion method are presented in this paper. Results of the analyses that were carried out in the study showed that the validity limits were perfectly consistent with the domain of the petrophysical interpretation. The results obtained by the application of the methodology to a real case, a deep-water exploration well data complicated by a poor characterization of the reservoir fluids, is also presented the paper.

This work is devoted to developing of comprehensive mechanistic model of heat transfer between an ESP motor and multiphase fluid. The model takes into account not only such parameters as flow pattern, gas/oil ratio, flow regime, etc. but also solid phase. Solid phase was included based on modern experimental investigations in the area of heat transfer in multiphase systems containing solid phase and flowing in annuli. The developed model was verified by means of temperature data taken using telemetry system during steady-state production from the wells with ESPs. Based on the model sensitivity analysis for the change of such parameters as solid size, solid phase concentration, and motor shroud size was performed. Impact of each parameter on motor cooling is discussed. It is shown that optimization of motor shroud size can lead to significant heat transfer enhancement between the ESP motor and solid-liquid mixture as compared to heat transfer between the ESP motor and liquid. This effect is achieved due to solid particles which have sufficient freedom in bombarding motor wall so that they thin viscous boundary layer. However, if the motor shroud size selection is inconsistent to the given solid phase concentration and solid particle size, this can cause deterioration of thermodynamic conditions for motor cooling.

During development of gas field one of the most serious problem is a well production decline due to borehole zone properties deterioration. For well production recovery there are a lot of various methods. And hydrochloric acid treatment received a large area of implementation. However it is a common practice that actual well production increase is lower then expected. This fact has a negative influence on the economic indexes and may bring to loss of profit. Estimation of treatment efficiency demand knowledge of such impotent value as relative increase of gas production. Previously believed that after borehole zone treatment the borehole zone permeability can change only trough cleaning and seam dissolution. But special research showed that after borehole zone treatment the compressibility of the seam may change too. In connection with above we build a mathematical model of gas filtration toward the well with complex borehole zone witch take into account both blockage and compressibility. In case when we consider borehole zone blockage only the filtration equation has analytical solution in terms of integral function. But when we account compressibility of the seam a numerical methods are needed. We considered relative permeability of gas as product of two functions k(r,p)=A(r)*f(p), where function A(r) simulate a borehole zone blockage and f(p) simulate compressibility of the seam. Previously function A(r) was considered as fully modeling. There are various formulation of its view such as linear function, quadratic function, root function and constant function. But all of these formulations haven’t a physical validation. We suggest using probable convolutions method with Gaussian kernel for well-logging data interpretation and definition redial borehole zone filtration properties. In our previous papers showed that a view of A(r) function greatly influence on well production. That is why an adequate definition of this function is very impotent. Using this model we researched dependence of hydrochloric acid treatment efficiency from such parameters as radius of acid penetration, permeability on the well, coefficient of compressibility. And we offered technological scheme of treatment efficiency prediction.

The multiphase flow in fractural media with different fractal structures is considered. The homogenized model of processes depends on the structure pattern. One of them leads to the known double porosity model, the other may be described by tensor total and pseudo phase permeability. The influence of the capillary effects is taken into consideration. The algorithms of the tensor permeability calculation are developed by using the self-similar properties of the fractal structure. We have used the support operator method for difference scheme construction. This method allows approximating both the dissipative energy and phase flux. The results of modeling some problems are represented.

The issue of systematic approach to well test design does not leave any doubts. It is particularly relevant for test design in an explored area with the complex geology where we usually possess the minimum of data about the reservoir and are required to derive the latter with sufficient confidence for the shortest time. Deconvolution technique allows us to decrease duration of well testing to the minimum albeit demands of linearity of the system investigated. In this work the thorough particular well test design is performed in the attempt to take into account all the issues concerning application of deconvolution technique and conventional test data interpretation. For design purposes the integration of simulations in Eclipse and Saphir was done in order to built unsophisticated but representative model in Saphir. Also the influence of aquifer and well interference were regarded. The interference was taken into account to linearize the system for deconvolution technique application. The method for time distribution between drawdown and build-up in the test sequence was suggested. Sensitivity analysis is implemented to reservoir and well parameters. Also the conventional test and test for deconvsolution are compared and economics issue are adduced.

Irreducible oil (IO) is researched as hydrodynamic system, consisting of separate types. IO of micro- and macrolevel, connate and conventionally mobile IO were defined. We identified that such technological factors as pressure gradient, wettability, surface tension, etc. effect IO significantly. Irreducible oil saturation (IOS) is one of the most important parameter characterizing development efficiency. Traditionally IOS was associated with formation microconstruction and considered through displacement efficiency using its relationship with filtration capacity properties of formation. IOS was considered as one of producing formation properties. Technology influence on final oil recovery used to be determined through volumetric efficiency. It was determined that conventionally mobile irreducible oil becomes movable under the critical pressure gradient to capillary pressure ratio. When the mobility threshold is reached the values of IOS depend on technological displacement condition. After achieving the second threshold value almost all conventionally mobile oil is withdrawn and only connate oil remains in formation. Present geo-technological models do not take the IOS influence of these factors into consideration. Therefore the dynamic models of IOS were developed with the attention paid on threshold mobility character, operating practices and mechanisms of IOS formations. Distribution research of IOS in croswell area was based on special hydrodynamic models. Results of investigation technology effect on IOS (e.g. well spacing system, pattern arrangement, disturbances and modernization of designed system) are presented. As shown, the compound anthropogenic heterogeneity of IOS is formed in crosswell area even for homogeneous formations. Depending on well spacing system IOS can be heterogeneous regarding displacement efficiency and phase permeability. We consider the possibility of additional IO recovery using intensive technologies (flow control, sidetracks, horizontal wells). Taking different well spacing systems, analysis of relative influence of formation microconstruction and of technological factors on IOS was given.

Nowadays nonconventional wells (horizontal, multilateral, Rad Tech and others) are being extensively drilled throughout the world for the development of low-profit fields. Construction of these wells enables to reduce filtration resistivity resulting in productivity index increase and costs reduction. To select the optimal well design with regard to reservoir characteristics, effective well operation and determination of filtration characteristics one should possess calculation methods for steady and unsteady liquid flow in reservoir. Few related papers have been published so far. However analytical methods for steady flowing are suitable for homogeneous beds with simple geometry and equal length laterals. Available approaches for description of pressure build up allow to account for various lateral trajectories but FEA or semi-analytical decisions methods are too labor-intensive for practical application. Commercial program for interpretation of pressure build up is lacking. Therefore simple methods of productivity index determination and pressure transient test interpretation are suggested for nonconventional wells. These methods are suitable for low thickness beds. The basis of these methods is the superposition of filtration resistivity for two plane problems. Trajectory of laterals is simulated as a number of closely spaced vertical wells or nodes. The suggested method allows determining the field of application and regularities for nonconventional wells. Dimensionless fluid-movement profile calculated from steady fluid flow and a superposition method for pressure builds up in the nodes are used for determination of pressure build up. For description of build up in a node we recommend a diffusion equation in Laplace space and Stephest numerical algorithm. The problem is solved for porous and porous-fractured reservoirs. Numerical calculations show that cross-flows occur after the horizontal or multilateral well shut-down. Pressure derivative maximum testifies to low effective length of the borehole or positive skin-effect. Knowledge of effective intervals length is critical to pressure curve interpretation.