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

In 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 O-method, 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 re-written 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 non-linear problems like the ones encountered in the study of multi-phase 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.

We propose a new single-phase local transmissibility upscaling method for adapted grids in 3-D domains that uses spatially varying and compact multi-point flux approximations (MPFA), based on the VCMP method previously introduced for 2-D 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 multi-point 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 two-point. The resulting MPFA stencil is spatially varying and reduces to two-point expressions in cases without full-tensor anisotropy. Numerical tests show that the method significantly improves upscaling accuracy as compared to commonly used local methods and also compares favorably with a local-global upscaling method. We also present a corrector method that adapts the stencils locally to guarantee that the resulting pressure matrix is an M-matrix. This corrector method is needed primarily for cases where strong permeability contrasts are mis-aligned 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 fine-scale features such as narrow channels.

For simulation of three-phase flow processes, such as water-alternating-gas injection (WAG), accurate descriptions of the three-phase 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 pore-scale network simulations. During multiple injection processes, especially under mixed-wet conditions, isolated clusters of the three phases occur. These clusters may still be mobilized, even under the assumption of capillary-dominated flow, as part of so-called 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, so-called 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 three-phase entry pressures, as well as the film and layer volumes and conductances in pore corners. Some WAG network simulations in mixed-wet 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.

Fractional-flow 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 fractional-methods to foam have excluded oil. Recently, Mayberry and Kam (SPE 100964) presented out fractional-flow solutions for foam injection with a constant effect of oil on the foam. We extend fractional-flow 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 three-phase SAG displacements examined, if foam collapses at the initial oil saturation in the reservoir, there is a very-small-velocity 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 three-phase displacements can be represented at two-phase displacements with chemical or miscible shocks: A model for a first-contact 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 fractional-flow methods.

Growing 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 time-step. Current reservoir simulators apply safe-guarded variants of Newton’s method, and often can neither guarantee convergence, nor provide estimates of the relation between convergence rate and time-step size. In practice, time-step 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 time-steps. We establish an alternate class of nonlinear iteration that both converges, and associates a time-step to each iteration. Moreover, the linear solution process within each iteration is performed locally. By casting the nonlinear residual for a given time-step as an initial-value-problem, we formulate a solution process that associates a time-step 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 large-scale simulation examples are presented, and the results demonstrate the robustness of the proposed method as well as its performance.

The 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 O-method 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 simulation is particularly efficient for incompressible two-phase flow, for which one can use front tracking to solve the 1-D 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 1-D 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 front-tracking methods that are based on a combination of solving ordinary Riemann problems and Riemann problems with discontinuous flux.

A 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 ill-conditioned 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) two-stage preconditioning method. To this end, we identify coefficient aggregates by a means of an efficient version of the Hoshen-Kopelman percolation algorithm suitable for any flow network structure. By partitioning and reordering unknowns according to these aggregates, we obtain a set of high-conductive blocks followed by a low-conductive 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 high-conductive 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 high-conductive blocks with deflation is realized as a two-stage 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 two-stage preconditioning can be used as a promising alternative to employ existing algebraic multilevel methods in a more effective way.

Control-volume 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 K-orthogonal grids, with a uniform permeability tensor, the scheme reduces to the standard two-point flux approximation. Monotonicity of the solution matrix is discussed and a version of the method that provides an M-matrix 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 two-point and multipoint flux approximations.

Upscaling is commonly applied to account for the effects of fine-scale permeability heterogeneity in coarse-scale simulation models. Ensemble level upscaling (EnLU) is a recently developed approach that aims to efficiently generate coarse-scale flow models capable of reproducing the ensemble statistics (e.g., P50, P10 and P90) of fine-scale flow predictions for multiple reservoir models. Often the most expensive part of standard coarsening procedures is the generation of upscaled two-phase flow functions (e.g., relative permeabilities). EnLU provides a means for efficiently generating these upscaled functions using stochastic simulation. This involves the use of coarse-block attributes that are both fast to compute and representative of the effects of fine-scale permeability on upscaled functions. In this paper, we establish improved attributes for use in EnLU, namely the coefficient of variation of the fine-scale single-phase velocity field (computed during computation of upscaled absolute permeability) and the integral range of the fine-scale 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 coarse-scale flow models. The procedure is tested on 3D well-driven flow problems with different (Gaussian) permeability distributions and high fluid mobility ratios. EnLU is shown to capture the ensemble statistics of fine-scale flow results (flow rate and oil cut as a function of time) with similar accuracy to full flow-based methods but at a small fraction of the computational cost.

Upscaling reservoir simulation models in an appropriate manner is important for any well test interpretation workflow which involves numerical reservoir simulation. Approaches based on steady-state 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.

Finite 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 L-method, which has been recently introduced by Ivar Aavatsmark and co-workers, displays enhanced monotonicity properties on distorted meshes and for anisotropic diffusion tensors. In this work, we present a generalization of the L-method 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 single-phase problems with distorted meshes and anisotropic and heterogeneous diffusion tensors. The proposed method is compared with the standard O- and L-methods as well as with an unconditionally symmetric, coercive cell centered finite volume scheme using a larger stencil.

Standard reservoir simulation schemes employ single-point 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 cross-wind diffusion. Novel higher dimensional upwind schemes that minimize cross-wind 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 full-tensor 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 cross-wind diffusion effect, leading to improved resolution of concentration and saturation fronts. TECHNICAL CONTRIBUTIONS Locally conservative unstructured multi-dimensional schemes coupled with full-tensor Darcy flux approximations are presented for reservoir simulation. The multi-dimensional 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 cross-wind diffusion.

This paper presents a new combined method for accurate upscaling of two-phase 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, vorticity-based 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 single-phase vorticity is attempted to be preserved among fine and coarse grid models. To demonstrate accuracy and efficiency, the combined DDM-vorticity 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.

Residual 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 non-wetting 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.

The Multiscale Finite-Volume (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 coarse-system 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 dual-coarse 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 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 fine-scale 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 fine-scale grid. To improve the multiscale method's predictive power for individual wells, we consider two direct opposite strategies: First, by resolving the near-well flow in the coarse grid by adaptive grid refinement in regions near wells. Second, by ensuring that near-well 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 near-well grid adaptations, and state pros and cons for the approaches considered.

In 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 in-situ combustion applied specially in fracture reservoirs where one is dealing with multi-scale multi-process 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 matrix-fracture 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 up-scaling method. The study is based on the fine grid, single porosity, multi-phase and multi-component 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 one-dimensional combustion are compared with the corresponding analytical solutions. The two-dimensional 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 in-situ 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 up-scaling procedure is proposed.

In 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 density-driven 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.

We investigate the potential of a flux-conservative, 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 (Self-adaptive time integration of flux-conservative equations with sources. J. Comput. Phys. 216 (2006)). We show that the scheme offers the advantage of avoiding the overly-restrictive 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 one-dimensional 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.

Seismic 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, in-filled with necessary interpolated data, to form the surface topography implemented. Recent results describing seismic sand curves for the region are employed in the near-surface medium definition, and a constant half space of realistic seismic velocities and densities are used for the simulated sub-sand 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 Finite-Differences (F-Ds) are used to discretize full elastic wave equations in the velocity-stress formulation, reducing gradually the F-D 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 F-Ds, and at the side and bottom boundaries exponential damping is applied for non-reflective 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 3-5 years has such complex modelling been commercially viable in hydrocarbon exploration.

Anisotropic effective medium theory (for homogenization of fractured/composite porous media under the assumption of complete scale-separation) can be used to map the effects of sub-grid (or sub-seismic) fractures onto the grid-scale. The effective poroelastic stiffness and hydraulic permeability tensors that determines the seismic amplitude versus angle-azimuth (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 physics-based 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 sub-grid 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 sub-set 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.

In 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 non-linear, 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 Karhunen-Loè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.

Since 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.

The 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.