 Home
 Conferences
 Conference Proceedings
 Conferences
ECMOR XI  11th European Conference on the Mathematics of Oil Recovery
 Conference date: 08 Sep 2008  11 Sep 2008
 Location: Bergen, Norway
 ISBN: 9789073781559
 Published: 08 September 2008
1  50 of 105 results


Continuous Darcyflux Approximations with Full Pressure Support for Structured and Unstructured Grids in 3D
Authors H. Zheng and M.G. EdwardsThis paper presents the development of families of controlvolume distributed fluxcontinuous schemes with full controlvolume surface pressure continuity for the general three dimensional pressure equation. Most effort to date has focused on schemes with pointwise continuity where the local continuity conditions are relatively compact. While full surface pressure continuity requires an increase in the number of local continuity conditions the new schemes prove to be relatively robust. Families of full pressure continuity schemes are developed here for structured and unstructured grids in three dimensions. The resulting Darcy flux approximations are applied to a range of three dimensional test cases that verify consistency of the schemes. Convergence tests of the threedimensional families of schemes are presented, for a range of quadrature points. Mmatrix conditions are presented and the schemes are tested for monotonicity. The full pressure continuity schemes are shown to be beneficial compared with earlier pointwise continuity schemes [1] both in terms improved monotonicity and convergence. The schemes are applied to challenging three dimensional test cases including heterogeneity for both structured and unstructured grids. TECHNICAL CONTRIBUTIONS Extension of full controlvolume surface pressure continuity flux continuous schemes to structured and unstructured grids in three dimensions. Mmatrix and monotonicity analysis of three dimensional methods. Benefits are demonstrated for challenging three dimensional test cases including heterogeneity on structured and unstructured grids.



A General Framework for Non Conforming Approximations of the Single Phase Darcy Equation
Authors D. Di Pietro, L. Agelas, R. Eymard and R. MassonIn this paper we present an abstract analysis framework for nonconforming approximations of heterogeneous and anisotropic elliptic operators on general 2D and 3D meshes. The analysis applies to discontinuous Galerkin methods, to the popular MPFA Omethod, as well as to hybrid finite volume methods. A number of examples is provided in the paper. The guidelines of the analysis can be summarized as follows: (i) each method is rewritten in weak formulation by introducing two discrete gradient operators. Penalty terms may be required for non $H^1$conformal spaces to ensure coercivity; (ii) a discrete $H^1$norm is introduced such that a discrete version of the Rellich theorem holds; (iii) an a priori estimate on the discrete solution in the discrete $H^1$norm is derived using the coercivity of the bilinear form; (iv) the convergence of the discrete problem to the continuous one is eventually deduced, thus proving the convergence of the approximation. In order for the above analysis to hold, the gradient operator applied to the discrete solution should be consistent, whereas the one applied to the test function should be weakly convergent in $L^2$. According to the method, further assumptions on the mesh may be required. This analysis framework easily extends to nonlinear problems like the ones encountered in the study of multiphase Darcy flows, and it allows to weaken regularity assumptions on the coefficients of the problem as well as on the exact solution. From a practical viewpoint, its main interest lies in the possibility to analyze and to design new methods as well as to improve existing ones.



On Efficient Implicit Upwind Schemes
Authors J.R. Natvig and K.A. LieAlthough 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 nontrivial 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 highorder 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.



Spatiallyvarying Compact Multipoint Flux Approximations for 3D Adapted Grids with Guaranteed Monotonicity
Authors J.V. Lambers and M.G. GerritsenWe propose a new singlephase local transmissibility upscaling method for adapted grids in 3D domains that uses spatially varying and compact multipoint flux approximations (MPFA), based on the VCMP method previously introduced for 2D Cartesian grids. For each cell face in the coarse upscaled grid, we create a local fine grid region surrounding the face on which we solve three generic local flow problems. The multipoint stencils used to calculate the fluxes across coarse grid cell faces involve up to ten neighboring pressure values. They are required to honor the three generic flow problems as closely as possible while maximizing compactness and ensuring that the flux approximation is as close as possible to being twopoint. The resulting MPFA stencil is spatially varying and reduces to twopoint expressions in cases without fulltensor anisotropy. Numerical tests show that the method significantly improves upscaling accuracy as compared to commonly used local methods and also compares favorably with a localglobal upscaling method. We also present a corrector method that adapts the stencils locally to guarantee that the resulting pressure matrix is an Mmatrix. This corrector method is needed primarily for cases where strong permeability contrasts are misaligned with the grid locally. The corrector method has negligible impact on the flow accuracy. Finally, we show how the computed MPFA can be used to guide adaptivity of the grid, thus allowing rapid, automatic generation of grids that can account for difficult geologic features such as faults and fractures and efficiently resolve finescale features such as narrow channels.



Analytical Estimation of Advanced IOR Performance
Authors R.A. Berenblyum, A.A. Shchipanov, E.M. Reich, L. Surguchev and K. PotschIn 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 inhouse analytical prescreening 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 fullscale reservoir models. The analytical methods can significantly decrease the decision making time by either: a) decreasing a number of computationally expensive fullscale simulations for mature developed reservoirs; or b) quickly screening new prospect where reservoir knowledge is limited.



Shortest Path Algorithm for Porescale Simulation of Wateralternatinggas Injection
Authors M.I.J. van Dijke, J.E. Juri and K.S. SorbieFor simulation of threephase flow processes, such as wateralternatinggas injection (WAG), accurate descriptions of the threephase capillary pressures and relative permeabilities as functions of the phase saturations are essential. Instead of deriving these functions from experiments, which is very difficult, they may be obtained from porescale network simulations. During multiple injection processes, especially under mixedwet conditions, isolated clusters of the three phases occur. These clusters may still be mobilized, even under the assumption of capillarydominated flow, as part of socalled multiple displacement chains of adjacent clusters stretching from inlet to outlet of the model. To determine the most favourable displacement chain requires implementation of an efficient shortest path algorithm. For the present problem the distances or costs are the capillary entry pressures between the various phase clusters. Because these entry pressures can also be negative, socalled negative cost cycles arise, which invalidate traditional algorithms. Instead, we have implemented an efficient shortest path algorithm with negative cycle detection If negative cycles arise, the corresponding cyclic displacement chains are carried out before any linear displacement from inlet to outlet. Negative cost cycles correspond to spontaneous displacements, for which the prevailing invading phase pressure is higher than required. Additionally, and probably even more significant, the shortest path algorithm determines in principle the pressures of all phase clusters in the network. These pressures are in turn used in the accurate calculation of the threephase entry pressures, as well as the film and layer volumes and conductances in pore corners. Some WAG network simulations in mixedwet media have been carried out using the new algorithm to demonstrate the occurrence of long displacement chains and the much improved efficiency, also for increasing network sizes. A significant number of negative cycles occurred, but the corresponding displacement chains were short. Additionally, we have analysed the obtained values of the cluster pressures.



Fractionalflow Theory Applied to NonNewtonian IOR Processes
Authors W.R. Rossen, R.T. Johns, K.R. Kibodeaux, H. Lai and N. Moradi TehraniThe method of characteristics, or fractionalflow 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 nonNewtonian rheology in polymer and foam IOR has been a serious limitation. We extend fractionalflow methods for twophase flow to nonNewtonian fluids in onedimensional cylindrical flow, where rheology changes with radial position r. The fractionalflow 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 nonNewtonian polymer and foam. For continuous injection of foam, one can map the 1D displacement in time and space exactly for the nonNewtonian foam bank, and approximately, but with great accuracy, for the gas bank ahead of the foam. The fractionalflow solutions are more accurate that finitedifference simulations on a comparable grid and can be used to calibrate simulators. Fractionalflow methods also allow one to calculate changing mobility near the injection well to a greater accuracy than simulation. For polymer and SAG (alternatingslug) foam injection, characteristics and shocks collide, making the fractionalflow 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 nonNewtonian in between, the nonNewtonian nature of the displacement depends on injection rate. The fractionalflow method extended to nonNewtonian flow can be useful both for its insights for scaleup of laboratory experiments and to calibrate computer simulators involving nonNewtonian IOR. It can also be an input to streamline simulations.



Fractionalflow Theory of Foam Displacements with Oil
Authors M. NamdarZanganeh, T. La Force, S.I. Kam, T.L.M. van der Heijden and W.R. RossenFractionalflow theory has proven useful for understanding the factors that control foam displacements and as a benchmark against which to test foam simulators. Most applications of fractionalmethods to foam have excluded oil. Recently, Mayberry and Kam (SPE 100964) presented out fractionalflow solutions for foam injection with a constant effect of oil on the foam. We extend fractionalflow methods to foam displacements with oil, using the effects of oil and water saturations on foam as represented in the STARS simulator. There can be abrupt shifts in the composition paths at the limiting water and oil saturations for foam stability. In the immiscible threephase SAG displacements examined, if foam collapses at the initial oil saturation in the reservoir, there is a verysmallvelocity shock from the injected condition to complete foam collapse. The displacement is nearly as inefficient as if no foam were present at all. It does not matter in these cases whether foam is weakened by low water saturation. The displacement is efficient, however, if foam is unaffected by oil but weakened at low water saturation. These results may reflect our foam model, where foam is only partially destroyed at low water saturations but is completely destroyed by high oil saturation. Two idealized models for threephase displacements can be represented at twophase displacements with chemical or miscible shocks: A model for a firstcontact miscible gas flood with foam suggests an optimal water fraction in foam that puts the gas front just slightly ahead of the foam (surfactant) front. An idealized model of a surfactant flood with foam for mobility control suggests it is important to inject a sufficiently high water fraction in the foam that the gas front is behind the surfactant front as the flood proceeds. We present simulations to verify the solutions obtained with fractionalflow methods.



Doublefamilies of QuasiPositive Fluxcontinuous Finitevolume Schemes on Structured and Unstructured Grids
Authors M.G. Edwards and H. ZhengThis paper focuses on fluxcontinuous pressure equation approximation for strongly anisotropic media. Previous work on families of flux continuous schemes for solving the general geometrypermeability tensor pressure equation has focused on singleparameter families e.g. [1]. These schemes have been shown to remove the O(1) errors introduced by standard twopoint 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 Mmatrices depending on the strength of the crossterms [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 doublefamily fluxcontinuous locally conservative schemes are presented for the general geometrypermeability tensor pressure equation. The new doublefamily formulation is shown to expand on the current singleparameter range of existing schemes that have an Mmatrix. While it is shown that a double family formulation does not lead to an unconditional Mmatrix scheme, an analysis is presented that classifies the sensitivity of the new schemes with respect to monotonicity and doublefamily quadrature. Convergence performance with respect to a range of doublefamily quadrature points is presented for some well known test cases.



AdaptivelylocalizedcontinuationNewton; Reservoir Simulation Nonlinear Solvers that Converge All the Time
Authors R.M. Younis, H. Tchelepi and K. AzizGrowing interest in understanding, predicting, and controlling advanced oil recovery methods emphasizes the importance of numerical methods that exploit the nature of the underlying physics. The Fully Implicit Method offers unconditional stability in the sense of discrete approximations. This stability comes at the expense of transferring the inherent physical stiffness onto the coupled nonlinear residual equations which need to be solved at each timestep. Current reservoir simulators apply safeguarded variants of Newton’s method, and often can neither guarantee convergence, nor provide estimates of the relation between convergence rate and timestep size. In practice, timestep chops become necessary, and are guided heuristically. With growing complexity, such as in thermally reactive compositional models, this can lead to substantial losses in computational effort, and prohibitively small timesteps. We establish an alternate class of nonlinear iteration that both converges, and associates a timestep to each iteration. Moreover, the linear solution process within each iteration is performed locally. By casting the nonlinear residual for a given timestep as an initialvalueproblem, we formulate a solution process that associates a timestep size with each iteration. Subsequently, no iterations are wasted, and a solution is always attainable. Moreover, we show that the rate of progression is as rapid as a standard Newton counterpart whenever it does converge. Finally, by exploiting the local nature of nonlinear waves that is typical to all multiphase problems, we establish a linear solution process that performs computation only where necessary. That is, given a linear convergence tolerance, we identify the minimal subset of solution components that will change by more than the specified tolerance. Using this a priori criterion, each linear step solves a reduced system of equations. Several challenging largescale simulation examples are presented, and the results demonstrate the robustness of the proposed method as well as its performance.



Study of a New Refinement Criterion for Adaptive Mesh Refinement in Sagd Simulation
Authors M. Mamaghani, C. Chainais and G. EncheryThe 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 smallsized 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 adaptivemeshrefinement 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.



Mimetic MPFA
Authors A.F. Stephansen and R.A. KlausenThe analysis of the multi point flux approximation (MPFA) method has so far relied on the possibility of seeing it as a mixed finite element method for which the convergence is then established. This type of analysis has been successfully applied to triangles and quadrilaterals, lately also in the case of rough meshes. Another well known conservative method, the mimetic finite difference method, has also traditionally relied on the analogy with a mixed finite element method to establish convergence. Recently however a new type of analysis proposed by Brezzi, Lipnikov and Shashkov, cf. [3], permits to show convergence of the mimetic method on a general polyhedral mesh. We propose to formulate the MPFA Omethod in a mimetic finite difference framework, in order to extend the proof of convergence to polyhedral meshes. The general nonsymmetry of the MPFA method and its local explicit flux, found trough a splitting of the flux, are challenges that require special attention. We discuss various qualities needed for a structured analysis of MPFA, and the effect of the lost symmetry.



Streamline Tracing on General Tetrahedral and Hexahedral Grids in the Presence of Tensor Permeability Coefficients
Authors S.F. Matringe, R. Juanes and H.A. TchelepiModern reservoir simulation models often involve geometrically complex (distorted or unstructured) threedimensional grids populated with fulltensor 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 fulltensor 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 fullfield reservoir model discretized by an unstructured grid and populated with heterogeneous tensor coefficients.



On Front Tracking for Compressible Flow
Authors H.M. Nilsen and K.A. LieStreamline simulation is particularly efficient for incompressible twophase flow, for which one can use front tracking to solve the 1D transport problems along streamlines. Here we investigate the extension of this method to compressible (and immiscible) flow and discuss some of the difficulties involved, and in particular the choices one has in writing the 1D transport equation(s). Our study is motivated by the simulation of CO2 injection, and we therefore also develop methods that are particularly suited for solving compressible flow where one phase is incompressible. Altogether, we present four fronttracking methods that are based on a combination of solving ordinary Riemann problems and Riemann problems with discontinuous flux.



Adaptive Control for Solver Performance Optimization in Reservoir Simulation
Authors I.D. Mishev, N. Fedorova, S. Terekhov, B.L. Beckner, A.K. Usadi, M.B. Ray and O. DiyankovRuntime 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 selfadaptive solver selection systems. For example SelfAdapting Largescale 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 online 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.



A Novel Percolative Aggregation Approach for Solving Highly IllConditioned Systems
Authors H. Klie and M.F. WheelerA key aspect in any algebraic multilevel procedure is to be able to reliably capture the physical behavior behind system coefficients. However, the modeling of complex reservoir scenarios generally involves the computation of extremely discontinuous and nonlinear coefficients that, in turn, compromise the robustness and efficiency of smoothing and coarsening strategies. In addition, the need of dealing with large discretized domains leads to highly illconditioned systems that represent a true challenge to any linear solver technology known today. In this work, we exploit the fact that flow trend information can be used as a basis for developing a robust percolative aggregation (PA) twostage preconditioning method. To this end, we identify coefficient aggregates by a means of an efficient version of the HoshenKopelman percolation algorithm suitable for any flow network structure. By partitioning and reordering unknowns according to these aggregates, we obtain a set of highconductive blocks followed by a lowconductive block. Diagonal scaling allows for weakening the interaction between these high and low conductive blocks plus improving the overall conditioning of the system. The solution of the highconductive blocks approximates well the solution of the original problem. Remaining sources of errors from this approximation are mainly due to small eigenvalues that are properly eliminated with a deflation step. The combination of the algebraic solution of the highconductive blocks with deflation is realized as a twostage preconditioner strategy. The deflation stage is carried out by further isolating the aggregate blocks with a matrix compensation procedure. We validate the performance of the PA approach against ILU preconditioning. Preliminary numerical results indicate that the PA twostage preconditioning can be used as a promising alternative to employ existing algebraic multilevel methods in a more effective way.



Mortar Coupling of Discontinuous Galerkin and Mixed Finite Element Methods
Authors G.V. Pencheva, S.G. Thomas and M.F. WheelerGeological 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 DGDG and DGMFE 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 subdomain grids can be nonmatching and the mortar grid can be much coarser making this a multiscale method. The resulting nonlinear algebraic system is solved via a nonoverlapping 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.



Multipoint Flux Approximations via Upscaling
Authors C.L. Farmer, A.J. Fitzpatrick and R. PotsepaevControlvolume discretizations for fluid flow in anisotropic porous media are investigated. The method uses numerical solutions of the simplest model equation, –div[K(x) grad p(x)] = f(x). The permeability tensor, K(x), is allowed to have discontinuities. Multipoint Flux Approximations (MPFA) are used, and transmissibility coefficients are obtained, in the usual way, from local numerical flow experiments (transmissibility upscaling) for each cell face. For regular Korthogonal grids, with a uniform permeability tensor, the scheme reduces to the standard twopoint flux approximation. Monotonicity of the solution matrix is discussed and a version of the method that provides an Mmatrix is described. This discretization scheme is applied to reservoir simulation on 3D structured grids with distorted geometry, highly anisotropic media, and discontinuities in the permeability tensor. Simulation results are presented and compared with results from other twopoint and multipoint flux approximations.



The Competing Effects of Discretization and Upscaling – A Study Using the qFamily of CVDMPFA
Authors M. Pal and M.G. EdwardsFamilies of fluxcontinuous, locally conservative, finitevolume schemes have been developed for solving the general tensor pressure equation on structured and unstructured grids [1,2]. A family of fluxcontinuous schemes is quantified by a quadrature parameterization [3]. Improved convergence has been observed for certain quadrature points [4]. In this paper a qfamily of fluxcontinuous (CVDMPFA) 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, controlvolume distributed, fulltensor 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 fluxcontinuous, finitevolume schemes for the general tensor pressure equation. IJNME, 2005. 4. M. Pal, Families of ControlVolume 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, HerriotWatt University, and M.J Blunt, Imperial College, Tenth SPE comparative Solution Project: A Comparison of Upscaling Techniques, 2001.



Ensemble Level Upscaling of 3D Welldriven Flow for Efficient Uncertainty Quantification
Authors Y. Chen, K. Park and L.J. DurlofskyUpscaling is commonly applied to account for the effects of finescale permeability heterogeneity in coarsescale simulation models. Ensemble level upscaling (EnLU) is a recently developed approach that aims to efficiently generate coarsescale flow models capable of reproducing the ensemble statistics (e.g., P50, P10 and P90) of finescale flow predictions for multiple reservoir models. Often the most expensive part of standard coarsening procedures is the generation of upscaled twophase flow functions (e.g., relative permeabilities). EnLU provides a means for efficiently generating these upscaled functions using stochastic simulation. This involves the use of coarseblock attributes that are both fast to compute and representative of the effects of finescale permeability on upscaled functions. In this paper, we establish improved attributes for use in EnLU, namely the coefficient of variation of the finescale singlephase velocity field (computed during computation of upscaled absolute permeability) and the integral range of the finescale permeability variogram. Geostatistical simulation methods, which account for spatial correlations of the statistically generated upscaled functions, are also applied. The overall methodology thus enables the efficient generation of coarsescale flow models. The procedure is tested on 3D welldriven flow problems with different (Gaussian) permeability distributions and high fluid mobility ratios. EnLU is shown to capture the ensemble statistics of finescale flow results (flow rate and oil cut as a function of time) with similar accuracy to full flowbased methods but at a small fraction of the computational cost.



Indirecterrorbased, Dynamic Upscaling of Multiphase Flow in Porous Media
Authors S.H. Lee, H. Zhou and H. TchelepiWe 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 finescale information that is reconstructed by a multiscale finite volume method (Jenny et al., JCP 217: 627641, 2006). Upscaling of multiphase 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 finescale 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 twophase flow in two dimensional heterogeneous domains. We demonstrate that the dynamically upscaled model achieves high numerical efficiency compared with finescale computations and also provides excellent agreement with reference finescale solutions.



Vertical Aggregation of Reservoir Simulation Models for Numerical Well Testing
Authors R.A. Archer and J. KohUpscaling reservoir simulation models in an appropriate manner is important for any well test interpretation workflow which involves numerical reservoir simulation. Approaches based on steadystate flow behaviour are not necessarily optimal when trying to capture the details of pressure transient phenomena. This study proposes a new layer aggregation strategy which forms a weighted sum of the permeability differences between cells in neighbouring layers. The weighting uses a kernel function (G) presented by Oliver (1990). This approach was successfully applied in four upscaling experiments using a permeability distribution based on the Tarbert formation. These cases included a standard well test, an interference test and a test from a partially penetrating well. The results showed that coarsened models with as few as five layers can reproduce the pressure transient behaviour of a twenty layer fine grid model.



Upscaling of Large Reservoir Systems by Using a Controlrelevant Approach
Authors S.A. VakiliGhahani, R. Markovinovic and J.D. JansenIn reservoir simulations, we use an ensemble of finescale 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 systemtheoretic 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 inputoutput behavior of the model. Therefore, we propose a “controlrelevantupscaling” (CRU) approach that determines equivalent coarsescale permeabilities based on the actual system’s inputoutput behavior. The coarsescale parameters are obtained as the solution of an optimization problem that minimizes the distance between the inputoutput behaviors of the fine and coarsescale 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 singlephase 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 reducedorder CRU algorithm that accelerate the upscaling procedure. In the cases considered, the CRU algorithm shows superior inputoutput behavior as compared to upscaling algorithms commonly used in reservoir simulators.



The MPFA G Scheme for Heterogeneous Anisotropic Diffusion Problems on General Meshes
Authors L. Agelas, D. Di Pietro, I. Kapyrin and R. MassonFinite volume cell centered discretizations of multiphase porous media flow are very popular in the oil industry for their low computational cost and their ability to incorporate complex nonlinear closure laws and physics. However, the trend towards a more accurate description of the geometry and of the porous medium requires to dispose of flux approximations handling general polyhedral meshes and full diffusion tensors. The convergence on general meshes as well as the robustness with respect to anisotropy and heterogeneity of the diffusion tensor should come at a reasonable computational cost. An important property on which the analysis relies is coercivity, which ensures the stability of the scheme and allows to prove convergence of consistent discretizations. Meeting all these requirements is still a challenge and this is an active field of research. The Lmethod, which has been recently introduced by Ivar Aavatsmark and coworkers, displays enhanced monotonicity properties on distorted meshes and for anisotropic diffusion tensors. In this work, we present a generalization of the Lmethod based on a discrete variational framework and on constant gradient reconstructions. The coercivity of the discretization requires a local condition depending both on the mesh and on the diffusion tensor to be satisfied. The monotonicity and convergence properties of the scheme are assessed on challenging singlephase problems with distorted meshes and anisotropic and heterogeneous diffusion tensors. The proposed method is compared with the standard O and Lmethods as well as with an unconditionally symmetric, coercive cell centered finite volume scheme using a larger stencil.



Monotonicity for Control Volume Methods on Unstructured Grids
Authors E. Keilegavlen and I. AavatsmarkIn 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 noflow 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 MPFAO 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.



Higher Dimensional Upwind Schemes for Flow in Porous Media on Unstructured Grids
Authors M.S. Lamine and M.G. EdwardsStandard reservoir simulation schemes employ singlepoint upstream weighting for approximation of the convective fluxes when multiple phases or components are present. These schemes rely upon upwind information that is determined according to the grid geometry. As a consequence directional diffusion is introduced into the solution that is grid dependent. The effect can be particularly important for cases where the flow is across grid coordinate lines and is known as crosswind diffusion. Novel higher dimensional upwind schemes that minimize crosswind diffusion are presented for convective flow approximation on structured and unstructured grids. The schemes are free of spurious oscillations and remain locally conservative. The higher dimensional schemes are coupled with fulltensor Darcy flux approximations. Benefits of the resulting schemes are demonstrated for classical convective test cases in reservoir simulation including cases with full tensor permeability fields, where the methods prove to be particularly effective. The test cases involve a range of unstructured grids with variations in orientation and permeability that lead to flow fields that are poorly resolved by standard simulation methods. The higher dimensional formulations are shown to effectively reduce numerical crosswind diffusion effect, leading to improved resolution of concentration and saturation fronts. TECHNICAL CONTRIBUTIONS Locally conservative unstructured multidimensional schemes coupled with fulltensor Darcy flux approximations are presented for reservoir simulation. The multidimensional schemes are developed for general unstructured grids. The schemes are proven to be positive subject to conditions on the tracing direction and permit higher CFL conditions than standard schemes. Comparisons with single point upstream weighting scheme are made on a range of unstructured grids for different grid orientation and aspect ratios in cases involving full tensor coefficient velocity fields. The numerical results demonstrate the benefits of the higher dimensional schemes both in terms of improved front resolution and significant reduction in crosswind diffusion.



Grid Optimization to Improve Orthogonality of Twopoint Flux Approximation for Unstructured 3D Fractured Reservoirs
More LessToday'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 twopoint 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 multiplepoint 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 controlvolume 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.



High Performance Dual Mesh Method to Simulated Twophase Gravity Dominated Flows in Porous Media
Authors M.A. Ashjari and B. FiroozabadiThis paper presents a new combined method for accurate upscaling of twophase displacements in highly heterogeneous reservoirs. The method has the capability to retain its high performance for various flow regimes, from viscous to gravity dominant displacements, without the need for further modifications and computational steps. Two different grids are incorporated for simulation. The grid on fine scale is used to recognize the complicated physics of flow which depends on dominated driving forces and their interaction with heterogeneity. However, to achieve a fast simulation, the global flow calculation is performed on the coarse scale grid using upscaled equivalent properties. The communication between two different scale grids is achieved by the dual mesh method (DMM) procedure. Since DMM performance is still dependent on the accuracy of the coarse scale simulation, vorticitybased coarse grid generation technique is also incorporated to limit the upscaling errors. The technique optimizes the coarse grid distribution based on vorticity preservation concept where singlephase vorticity is attempted to be preserved among fine and coarse grid models. To demonstrate accuracy and efficiency, the combined DDMvorticity method is applied to highly heterogeneous systems in two dimensions with and without gravity. The results reveal that the flow regime has only minor impact on the performance of the combined method.



Discretisation Schemes for Anisotropic Heterogeneous Problems on Nearwell Grids
Authors S.S. Mundal, E. Keilegavlen and I. AavatsmarkIn 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 nearwell flow models. Nearwell regions are characterised as high flow density regions, and the dominating flow pattern exhibits a radiallike 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 nearwell regions are needed. The models should also be adapted to handle complex geological nearwell structures and multiphase flow simulations. Existing nearwell models are in general based on homogeneous media in the well vicinity. Anisotropies and strong heterogeneities are less accounted for in nearwell flow simulations. In this work, we construct analytical solutions for nearwell 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 radialtype grids for such cases and give their convergence behavior. The objective is to obtain a clearer view on the accuracy of the nearwell 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 nearwell conditions. The simulations are performed for singlephase flow, however, the results will also have relevance for multiphase flow simulations.



A stochastic trapping model for nonwetting phase ganglias in two phase flow through porous media
Authors M. Tyagi, H.A. Tchelepi and P. JennyResidual trapping is one of the CO2 retention mechanism during CO2 storage in saline aquifers. In a region where the CO2 saturation falls below the residual saturation, isolated ganglias of CO2 are formed which can become immobile (trapped). In this paper, we propose a stochastic model for trapping and release of nonwetting phase ganglias. Opposed to traditional trapping models, which only take into account of the mean quantities, the proposed stochastic model contains additional statistical informations, e.g., correlation time. By selecting appropriate model parameters, both drainage and imbibition can be described without a priori knowledge of the displacement direction.



A Twoscale Operator Splitting for Flow Problems in Porous Media
Authors Y. Cao, A. Weiss and B. WohlmuthOperator splitting time discretization techniques are getting more and more popular for the numerical simulation of nonlinear 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 semiLagrangian method for the advective part. Hence, the computational cost of the semiLagrangian method can be neglected compared to the one of the implicit finite volume scheme. In the second part, we introduce a twoscale 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 twoscale 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.



A General Multiscale FiniteVolume Method for Compressible Multiphase Flow in Porous Media
Authors H. Hajibeygi and P. JennyThe Multiscale FiniteVolume (MSFV) method was originally proposed to solve elliptic problems arising from incompressible multiphase flow in highly heterogeneous porous media at reduced computational cost. However, when phases are compressible, especially when the reservoir contains gas, mathematical formulations lead to a parabolic equation to be solved for pressure. In this paper we introduce a general MSFV method to deal with such parabolic problems. In this scheme, the basis and correction functions are numerical solutions of the full parabolic problems in localized domains. Hence compressibility effects are represented by a stencil in the coarsesystem matrix, i.e. not only by a diagonal entry. Furthermore, to enhance the computational efficiency of the scheme, the basis functions are kept independent of the pressure field. As a result, only the correction functions (1 per dualcoarse cell) have to be updated during the iterative procedure. It is an important property of this approach that it requires no additional simplifications. Finally, its good efficiency is demonstrated for a number of challenging test cases.



Multiscale Asynchronous Algorithms Based on the Superelements Method for Multiphase Flow
Authors A.K. Pergament, V.A. Semiletov and P.Y. TominThe 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 BuckleyLeverett equations. In general case it is possible to construct the tensor total permeability coefficients and to upscale the relative permeability. The method of calculating these variables is the dissipative energy integral approximation. As a result we have the nonlinear 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.



Multiscale/Mimetic Pressure Solvers with Nearwell Grid Adaptation
Authors B. Skaflestad and S. KrogstadMultiscale methods have proved to be an efficient means to achieving fast simulations on large reservoir models. The main characteristic of these methods is high resolution flow fields obtained at relatively low computational cost. We are thus able to resolve large scale flow patterns as well as finescale details that would be impossible to obtain for models for which direct simulation using traditional methods would be prohibitive. However, there are still a number of open problems in applying these methods to reservoir simulation. In particular, we observe some discrepancy in the performance of wells when compared to direct simulation on the finescale grid. To improve the multiscale method's predictive power for individual wells, we consider two direct opposite strategies: First, by resolving the nearwell flow in the coarse grid by adaptive grid refinement in regions near wells. Second, by ensuring that nearwell flow is sufficiently captured in the corresponding multiscale basis functions. For the latter strategy, we consider both adaptive alignment of coarse pressure grid to well trajectories, and an oversampling method for the computation of the multiscale basis functions corresponding to wells. In this paper we will study the effects of such nearwell grid adaptations, and state pros and cons for the approaches considered.



Spontaneous Ignition in Porous Media at Long Times
Authors J. Bruining and D. MarchesinIf 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.



Study of Heavy Oil Recovery from a Fractured Carbonate Reservoir Using in Situ Combustion
Authors H. Fadaei, M. Quintard, G. Debenest, G. Renard and A. KampIn spite of its strategic importance, the topic of recovery of heavy crude oils from fractured carbonate reservoir has not been extensively addressed. Thermal methods seem well suited for this kind of problems, particularly in situ combustion has shown promising results in laboratory experiments. Extensive work has been done on development of thermal process simulator but for the insitu combustion applied specially in fracture reservoirs where one is dealing with multiscale multiprocess problem, many unknowns are still exist. The recovery mechanism, reservoir and operational conditions on which the combustion can propagate in fractured systems are not enough clear. Also due to safety issues, air injection required careful assessment of the reservoir displacement mechanisms in particular the magnitude and the kinetics of matrixfracture transfers. To understand the mechanism of heavy oil recovery from a fractured reservoir we propose the development of a numerical simulation strategy, starting from existing simulation tools that are adapted to this particular problem. This will allow firstly understanding the role of each driving mechanism and physical as well as operational parameters in recovery process and secondly choosing the suitable upscaling method. The study is based on the fine grid, single porosity, multiphase and multicomponent simulation of a core surrounded by two parallel air invaded fractures using the thermal simulator. Firstly the simulator is validated for different processes: one and two dimensional diffusion and onedimensional combustion are compared with the corresponding analytical solutions. The twodimensional combustion is validated using experimental results available in the literature. The simulation results predict the conditions on which the combustion is sustained in the fractured reservoir as a function of oxygen diffusion coefficient, injection rate, and the permeability of the matrix. Oil production via natural drainage, hot fluid injection and insitu combustion are compared to address the importance of different driving mechanisms. At the block scale the effect of fracture spacing, heterogeneity in the matrix and the grid size on the propagation of combustion and the oil production are studied and then a suitable upscaling procedure is proposed.



Modeling CO2 Sequestration Using a Sequentially Coupled 'IterativeIMPECtimesplitthermal' Compositional Simulator
Authors M. Delshad, S.G. Thomas and M.F. WheelerThis 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 twophase flash algorithm. An implicitpressure, explicitconcentrations (IMPEC) sequential algorithm is then applied iteratively to enforce the nonlinear 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 timesplit 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 largescale 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 semianalytic solutions, wherever available. All computations were carried out within the IPARSCO2 code base framework.



Densitydriven Natural Convection in Dual Layered and Anisotropic Porous Media with Application for CO2 Injection Project
Authors R. Farajzadeh, F. Farshbaf Zinati, P.L.J. Zitha and J. BruiningIn this paper we investigate the mass transfer of CO2 injected into a layered and anisotropic (sub)surface porous formation saturated with water. Solutions of carbon dioxide in water and oil are denser than pure water or oil. We perform our analysis to a rectangular part of the porous medium that is impermeable at the sides except at the top, which is exposed to CO2. For this configuration densitydriven natural convection enhances the mass transfer rate of CO2 into the initially stagnant liquid. The analysis is done numerically using mass and momentum conservation laws and diffusion of CO2 into the liquid. This configuration leads to an unstable flow process. Numerical computations do not show natural convection effects for homogeneous initial conditions. Therefore a sinusoidal perturbation is added for the initial top boundary condition. It is found that the development of fingers is fastest for mass transfer enhancement by natural convection is largest for large anisotropy ratio’s and smaller for small ratio's. It is found that the mass transfer increases and concentration front moves faster with increasing Rayleigh number if the high permeability layer is on top. Of particular interest is the case when the Rayleigh number for the high permeable layer is above the critical Rayleigh number (Racr = 40) and smaller than Racr for the low permeable layer. The results of this paper have implications in enhanced oil recovery and CO2 sequestration in aquifers.



Operator Splitting of Advection and Diffusion on Nonuniformly Coarsened Grids
Authors V.L. Hauge, J.E. Aarnes and K.A. LieHighresolution 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 nonuniform coarsening strategy to upscale geomodels, where the coarse grid is generated by grouping cells in the finescaled 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 advectiondominated twophase 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.



Asynchronous Time Integration of Fluxconservative Transport
Authors B.T. Mallison, M.G. Gerritsen and G. KreissWe investigate the potential of a fluxconservative, asynchronous method for the explicit time integration of subsurface transport equations. This method updates each discrete unknown using a local time step, chosen either in accordance with local stability conditions, or based on predicted change of the solution, as in Omelchenko and Karimabadi (Selfadaptive time integration of fluxconservative equations with sources. J. Comput. Phys. 216 (2006)). We show that the scheme offers the advantage of avoiding the overlyrestrictive global CFL conditions. This makes it attractive for transport problems with localized time scales, such as those encountered in reservoir simulation where localized time scales can arise due to well singularities, spatial heterogeneities, moving fronts, or localized kinetics, amongst others. We conduct an analysis of the accuracy properties of the method for onedimensional linear transport and first order discretization in space. The method is found to be locally inconsistent when the temporal step sizes in two adjacent cells differ significantly. We show numerically, however, that these errors do not destroy the order of accuracy when localized. The asynchronous time stepping compares favorably with a traditional first order explicit method for both linear and nonlinear problems, giving similar accuracy for much reduced computational costs. The computational advantage is even more striking in two dimensions where we combine our integration strategy with an IMPES treatment of multiphase flow and transport. Global time steps between pressure updates are determined using a strategy often used for adaptive implicit methods. Numerical results are given for immiscible and miscible displacements using a structured grid with nested local refinements around wells.



High Order Adaptive Implicit Method for Conservation Laws
Authors A. Riaz, R. de Loubens and H. TchelepiNumerical 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 firstorder 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 highorder spatial discretization in both the implicit and explicit cells. Our novel approach solves the fully coupled transport and flow equations with high orderaccuracy within the AIM framework to reduce the computational cost and improve the solution accuracy. We show that highorder discretization is computationally efficient when combined with AIM and leads to significant improvements in accuracy compared to the firstorder solution.



Effects of Sand Dunes on Full Wavefields in Exploration
More LessSeismic pressure point sources have been simulated near the surface of sand dunes in Oman. Real topographic elevation data was assembled from recorded field data in the region, infilled with necessary interpolated data, to form the surface topography implemented. Recent results describing seismic sand curves for the region are employed in the nearsurface medium definition, and a constant half space of realistic seismic velocities and densities are used for the simulated subsand medium. Topographic effects in the form of scattering and reflections from prominent sand dune topography are confirmed, showing differences between scattering from elongated and normally oriented sand dunes to the excited source orientation. 8th order FiniteDifferences (FDs) are used to discretize full elastic wave equations in the velocitystress formulation, reducing gradually the FD order when approaching all numerical grid boundaries. At the free surface, I employed boundary conditions for arbitrary surface topography derived for the particle velocities, discretized by 2nd order FDs, and at the side and bottom boundaries exponential damping is applied for nonreflective boundaries. A rectangular grid is used to model the Cartesian wave equations from a curved grid by adding in the equations extra terms that depend on the arbitrary surface topography and its spatial derivatives. The curved grid is adapted to an arbitrary surface topography with gradually less curvature with depth towards the bottom of the grid, which is planar. The code is parallellized by domain decomposition through MPI (Message Passing Interface) parallellization between processors. This makes possible the time domain modelling of higher frequencies (even relevant for exploration) and larger areas than could be done otherwise. Only in the last 35 years has such complex modelling been commercially viable in hydrocarbon exploration.



Reservoir Characterization Improvement by Accounting for Well Test Extracted Effective Permeabilities
Authors A. Skorstad, F. Georgsen, P. Abrahamsen and E. SmørgravGeostatistical 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 cokriging 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 cokriging. 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.



Anisotropic Permeability in Fractured Reservoirs from Joint Inversion of Seismic AVAZ and Production Data
Authors M. Jakobsen and A. ShahrainiAnisotropic effective medium theory (for homogenization of fractured/composite porous media under the assumption of complete scaleseparation) can be used to map the effects of subgrid (or subseismic) fractures onto the gridscale. The effective poroelastic stiffness and hydraulic permeability tensors that determines the seismic amplitude versus angleazimuth (AVAZ) and production data (provided that one has information about the overburden and the porous matrix, as well as suitable tools for seismic forward modelling and fluid flow simulation) can therefore be viewed as functions of a relatively small number of parameters related with the fractures (e.g., fracture density and azimuthal orientation). In this paper, we develop a rock physicsbased Bayesian method for joint inversion of seismic AVAZ and production data with respect to these parameters of the fractures. Within this stochastic framework, the expectation value of the effective permeability tensor is given by an integral of an effective medium approximation (derived by using rigorous integral equation or Green’s tensor function methods) weighted by the posterior probability density function over the allowed values of the parameters of the fractures. The present work complements those of Will et al. (2005), Jakobsen et al. (2007) and Jakobsen and Shahraini (2008) in the sense that we are using different inversion methods, seismic attributes and/or (rock physics) models of the fractured reservoir. Our interdisciplinary method for characterization of subgrid fractures can in principle be applied to rather complicated models of fractured reservoirs (e.g., involving multiple sets of fractures characterized by different fracture densities and orientations, etc.). However, the initial (synthetic) modelling and inversion results presented here are associated with a relatively simple model of a fractured reservoir in which a single set of vertical fractures are located within subset of the (piecewise constant) reservoir model (that can in principle be identified using complementary methods). We have analysed synthetic AVAZ and production data contaminated with different amounts of normal distributed noise, and managed to recover the (unknown parameters of the fractures) effective permeability tensors with good accuracy. The results show clearly that the incorporation of seismic anisotropy attributes into the history matching process helps to reduce the uncertainties of the estimated permeability tensors.



A Novel Kriging Approach for Incorporating Nonlinear Constraints
Authors F.P. Campozana and L.W. LakeThis work presents new Krigingbased 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 welltestderived permeability; the second one enables one to generate fields that have a desired variability. Therefore, a wellknown 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 fieldgenerated permeability within the welltest 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 welltest permeability. The simulation results proved that the fields generated by the algorithm actually matched welltest permeability and field data variance within a reasonable tolerance. The same procedure can be used to incorporate other nonlinear constraints, such as facies proportion.



History Matching of Truncated Gaussian Models by Parallel Interacting Markov Chains on a Reduced Dimensional Space
More LessIn oil industry and subsurface hydrology, geostatistical models are often used to represent the spatial distribution of different lithofacies in the reservoir. Two main model families exist: multipoint and truncated Gaussian models. We focus here on the latter. In history matching of lithofacies reservoir model, we attempt to find multiple realizations of lithofacies configuration, that are conditional to dynamic data and representative of the model uncertainty space. This problem can be formalized in the Bayesian framework. Given a truncated Gaussian model as a prior and the dynamic data with its associated measurement error, we want to sample from the conditional distribution of the facies given the data. A relevant way to generate conditioned realizations is to use Markov Chains Monte Carlo (MCMC). However, the dimension of the model and the computational cost of each iteration are two important pitfalls for the use of MCMC. In practice, we have to stop the chain far before it has scanned the whole support of the posterior. Further more, as the relationship between the data and the random field is nonlinear, the posterior can be multimodal. Hence, the chain may stay stuck in one of the modes. In this work, we first show how to reduce drastically the dimension of the problem by using a truncated KarhunenLoève expansion of the Gaussian random field underlying the lithofacies realization. Then we show how we can improve the mixing properties of classical single MCMC, without increasing the global computational cost, by the use of parallel interacting Markov chains at different temperatures. Applying the dimension reduction and this innovative sampling method lowers drastically the number of iterations needed to sample efficiently from the posterior. We show the encouraging results obtained when applying the methodology to a synthetic history matching case.



Using the EnKF with Kernel Methods for Estimation of NonGaussian Variables
Authors J.G. Vabø, G. Evensen, J. Hove and J.A. SkjervheimThe 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 nonGaussian and multimodal, 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 pluriGaussian approach. Results from a field case are shown.



History Matching Using a Multiscale Ensemble Kalman Filter
Authors W. Lawniczak, R.G. Hanea, A.W. Heemink, D. McLaughlin and J.D. JansenSince the first version of Kalman Filter was introduced in 1960 it received a lot of attention in mathematical and engineering world. There are many successful successors like for example Ensemble Kalman Filter (Evensen 1996) which has been applied also for reservoir engineering problems. The method proposed in [Zhou et al. 2007] draws together the ensemble filtering ideas and an efficient covariance representation, and is expected to perform well in history matching for reservoir engineering. It is the Ensemble Multiscale Filter. The EnMSF is a different way to represent the covariance of an ensemble. The computations are done on a tree structure and are based on an ensemble of possible realizations of the states and/or parameters of interest. The ensemble consists of replicates that are the values of states per pixel. The pixels in the grid are partitioned between the nodes of the finest scale in the tree. A construction of the tree is led by the eigenvalue decomposition. Then, the state combinations with the greatest corresponding eigenvalues are kept on the higher scales. The updated states/parameters using the EnMSF are believed to keep geological structure due to localization property. It comes from the filter’s characterization where the pixels from the grid (e.g. permeability field) are distributed (in groups) over the finest scale tree nodes. We present a comparison of covariance matrices obtained with different setups used in the EnMSF. This sensitivity study is necessary since there are many parameters in the algorithm which can be adjusted to the needs of an application; they are connected to the tree construction part. The study gives the idea of how to efficiently use the EnMSF. The localization property is discussed based on the example where the filter is run with a simple simulator (2D, 2 phase) and a binary ensemble is used (the pixels in the replicates of permeability have two values only). Several possible patterns for ordering the pixels are applied.



Comparing different ensemble Kalman filter approaches
Authors B.V. Valles and G. NaevdalOver 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 KolmogorovSmirnov test, that for the PUNQS3 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 PUNSS3 model against the traditional EnKF. This approach seems better than paired or coupled EnKF to help solve the observed inconsistencies.



Channel Facies Estimation Based on Gaussian Perturbations in the EnKF
Authors D. Moreno, S.I. Aanonsen, G. Evensen and J.A. SkjervheimThe ensemble Kalman filter (EnKF) method has proven to be a promising tool for reservoir model updating and history matching. However, because the EnKF requires that the prior distributions for the parameters to be estimated are Gaussian, or approximately Gaussian, the application to facies models, and channel systems in particular, has been a challenge. In this paper we suggest two different approaches for parameterization of the facies models in terms of Gaussian perturbations of an existing "best guess" model. Method 1 is inspired by level set methods, where surfaces (here facies boundaries) are implicitly modelled through a level set function, normally defined as the signed distance from the nearest surface. Model realizations are generated by adding a Gaussian random field to the initial level set function. Method 2 is based on a similar idea, but the realizations are generated by adding a Gaussian random field to a smoothed indicator function. The smoothing is performed with a Shapiro filter, which is fast and simple to use for any number of dimensions. The performance of the methods is illustrated using a 3D model inspired by a real North Sea fluvial reservoir. It is shown that realistic facies model realizations may be generated from realizations of a Gaussian random field. Based on one realization from the prior for each of the methods, two sets of synthetic production data were generated. Then the prior model ensemble generated with each of the methods were conditioned to each of the two data sets using EnKF. Reasonably good matches were obtained in all cases, including those where the true model is not a realization from the same statistical model as is used to generate the prior realizations.
