ECMOR XII - 12th European Conference on the Mathematics of Oil Recovery

Upscaling is often applied to coarsen highly-detailed geological descriptions. The upscaling of multiphase flow functions is challenging due to their strong dependency on global flow effects. In this work, we present two approaches to generate upscaled two-phase flow functions with global flow effects incorporated. In global two-phase upscaling, the upscaled two-phase functions are directly computed from global fine-scale two-phase solutions. In a new adaptive local-global two-phase upscaling approach, local boundary conditions (for both pressure and saturation) are determined from time-dependent global coarse-scale solutions. This avoids solving global fine-scale two-phase flow. In both approaches, the upscaled two-phase flow functions are adapted to a specific flow scenario, thus providing more accurate coarse models than the use of generic flows. The methods are applied to heterogeneous permeability fields with full-tensor anisotropy. We demonstrate the impact of full-tensor effects on the global flow dependency in two-phase upscaling. It is also shown that the full-tensor effects in two-phase flow are appropriately accounted for by the upscaled two-phase functions with the global flow effects incorporated. The use of those upscaled functions in conjunction with a two-point flux approximation can effectively capture the two-phase crossflow due to the full-tensor anisotropy on the coarse scale.

One of the challenges for reservoir simulation is numerical dispersion. For waterflooding applications the effect is controlled due to the self-sharpening nature of a Buckley-Leverett shock. In this paper, we show that for augmented waterflooding – due to the coupling of compositional dispersion with fractional flow – the trailing shock is no longer self-sharpening. Thus, the simulation of such processes suffers from even higher numerical dispersion effects compared to pure waterflooding. Rather than implementing a higher-order discretization method, we propose a simple scheme based on segregation within a grid block. Compared to current mixing schemes, it differs in that segregation not only affects fluid properties but the transport, too. The scheme is shown to reestablish self-sharpness across the trailing shock. Simulations are performed that demonstrate the effectiveness of the scheme. The simulations also illustrate the possible effects of numerical dispersion on the predictions of augmented waterflooding. Finally, we discuss the extension of this technique to compositional simulation through a coupled limited-flash/segregation-in-flow assumption. Preliminary results demonstrate the potential of the approach as a heuristic method to control numerical dispersion for the simulation of miscible and near-miscible gas injection processes.

Cavities and fractures can significantly affect the flow paths of carbonate reservoirs and should be accurately accounted for during flow simulation. Herein, our goal is to compute the effective permeability of rock samples based on high-resolution 3D CT-scans containing millions of voxels. Hence, we need a flow model that properly accounts for the effects of Darcy flow in the porous material and Stokes flow in the void volumes on relevant scales. The presence of different length scales and large contrasts in the petrophysical parameters lead to highly ill-conditioned linear systems that make such a flow model very difficult to solve, even on large-scale parallel computers. To identify simplifications that render the problem computationally tractable, we analyze the relative importance of the Stokes and Darcy terms for a wide variety of parameter ranges on an idealized 2D model. We find that a system with a through-going free flow region surrounded by a low permeable matrix can be accurately modeled by ignoring the Darcy matrix and simulating only the Stokes flow. Using this insight, we are able to compute the effective permeability of a specific model from a CT-scan that contains more than eight million voxels.

A new analytical upscaling method is described. It is applicable for stratified reservoirs or with vertically distributed properties. The inter-layer communication is assumed to be perfect, corresponding to viscous or gravity dominant regimes. We apply asymptotic analysis to 2D flow equations of two incompressible immiscible fluids and reduce 2D Buckley-Leverrett system to 1D flow problem. Flow velocities in the layers depend on saturations of other layers, which reflects interaction between layers; the exchange terms between the layers may be expressed explicitly. The resulting quasilinear hyperbolic equations allow a self-similar solution, which makes it possible to express the pseudo-fractional flow function as a function of average saturation. The system is solved by the finite difference method. Both cases of discontinuous and continuous (lognormal) distributed layer permeabilities are studied. As comparison, the complete 2D waterflooding problems with very good vertical communication are solved in COMSOL. Generally, saturation profiles of the pseudo 1D model are only slightly different from the 2D simulation result by Comsol. The difference between recovery curves is marginal. This proves validity of the proposed method. Our method is much more advantageous over the classical Hearn-Kurbanov procedure. It is better at predicting displacement fronts and oil recovery.

Carbonate reservoirs contain a large fraction of the remaining oil reserves, but their highly heterogeneous structures, at the pore scale, make them difficult to characterize and this is one of the main reasons for their low recovery factors. Pore space imaging and reconstruction techniques are often used to characterize the pore space topology and geometry. However, the wide range of pore sizes encountered in carbonates, including submicron microporosity, usually renders these techniques unsuitable. A classic technique to approximate the pore size distribution (Ritter and Drake) is based on inversion of the mercury intrusion capillary pressure (MICP) curve. However, this technique is inaccurate due to the lack of accessibility of the large pores before the percolation threshold is reached. In this paper, we use a simple lattice network, characterized by an average co-ordination number, a pore volume-vs-radius power law and a pore (radius) size distribution (PSD) with an arbitrarily large range of radii, to match the MICP curve. Efficient optimization methods have been implemented to estimate these principal parameters, as outlined below. The overall aim of this work is to predict rock flow functions for different floods or saturations paths that are difficult and expensive to obtain from laboratory experiments, based on the PSD estimations. In the present approach the PSD has a given number of non-constant radius classes (bins), for which the probabilities need to be determined. As an initial guess, we take the Ritter and Drake PSD, based on which the non-fixed bin widths are determined, using principal theorems derived in information theory and data compression. Then, the pore-network parameters including the PSD class probabilities are optimized to match the experimental MICP, using two methods that are based on the link between the Boltzmann - Maxwell equation describing the molecular state probabilities in statistical physics and the Bayes method for posterior distribution. The first method is the classic Metropolis Markov Chain Monte Carlo method and the second, hybrid, method adds Hamiltonian Dynamics to avoid the random nature of the classic Metropolis algorithm, thus enhancing the convergence of the classical approach. The methods have successfully been applied to two synthetic cases and three carbonate samples from a Middle East reservoir presenting three different carbonate facies. The hybrid method improves the Markov Chain Monte Carlo efficiency in terms of the accepted system parameters states by a factor of 4 to 5.

It is commonly understood that carbonate reservoirs hold much of the worlds' oil and gas reserves. In this study we examine a heterogeneous vuggy pore class. The vuggy pore space is imaged using x-ray computed tomography (CT). Micro-CT ($mu$CT) is used to better image parts of the interparticle (matrix) porosity. Image segmentation is applied to obtain a three-dimensional (3D) map of the vuggy pore space, which is used to determine its connectivity. We evaluate the performance of selected image segmentation techniques applied for segmenting CT-scans of vuggy carbonate rocks at varying resolutions. Computations are performed for both two-dimensional (2D) and three-dimensional (3D) image segmentation. The computations show that the methods perform well on high resolution $mu$CT scans. For a CT scan with low and non-uniform resolution they wrongly classify a large portion of the voxels. This makes it challenging to determine the connectivity. The computed porosity is, due to image resolution, underestimated when compared to measured porosity from laboratory experiments.

Accurate predictions of gas-injection processes usually require a compositional formulation based on an Equation of State (EoS). Because the thermodynamic behavior of multi-component multiphase systems is highly nonlinear and coupling to the flow equations is quite complex, EoS-based simulations can become computationally prohibitive. For immiscible displacements an efficient approximation of the general compositional problem can be used if the phase-equilibrium ratios (K-values) are assumed to be weak functions of composition; however, this assumption is questionable for near-miscible displacements. The standard K-value approach suffers from significant difficulties, both in terms of robustness and accuracy, in the critical region. Here, we describe an extended K-values method that takes advantage of the solution-path invariance in the compositional space with respect to the hydrodynamic properties. Specifically, we propose an additional degree-of-freedom, which captures the composition dependence of the phase behavior, for use in the tabulation and interpolation of the K-values. An important aspect of the method is the use of the so-called Minimal Critical Pressure criteria (MCP), which indicates when a given composition becomes super-critical. We compare results obtained with EoS- and K-values based simulations for several (isothermal) compositional problems, and we demonstrate the efficiency and accuracy of the proposed method.

Modeling the phase behavior associated with compositional flow simulation of systems that form more than two phases (e.g., steam, or CO2 injection) is a challenging problem. In addition to the coupling of thermodynamics with the nonlinear equations of flow and transport, accurate phase-state identification for mixtures that form three (or more) phases, as function of pressure, temperature and composition, is the subject of active research. We describe a general negative-flash method for multi-component thermal systems that can form three, or more, fluid phases, including a convergence proof. This extended negative-flash approach integrates quite nicely with flow simulation based on adaptive tie-simplex compositional space parameterization. We prove that tie-simplexes change continuously with pressure, temperature, or along a continuous trajectory in compositional space. Continuation of the tie-simplex space provides theoretical justification for the tie-simplex based flow modeling, in which interpolation in pressure and temperature for a limited number of tie-simplexes is performed during a simulation run. We also show how a tie-simplex associated with a given composition can degenerate as a function of pressure and temperature. We focus on the complex behaviors of the tie-triangles and tie-lines associated with three-phase, steam-injection problems. The algorithms that capture the degeneration of the tie-triangles into tie-lines are described in detail. Interpolation in the parameterized compositional space is used to identify the phase-state and proves to be as reliable as a 'conventional' phase-stability test for three-phase mixtures. This extended negative-flash algorithm has been integrated into our tie-simplex based flow simulation framework. We demonstrate the effectiveness of the framework using several compositional steam-injection problems with complex behaviors.

Temperature modeling and analysis of horizontal wells has become increasingly important due to the extensive deployment of downhole temperature monitoring devices. Such sensors allow continuous monitoring of zonal production with a sufficiently high resolution. Published descriptions of the necessary analysis tools have concentrated on well temperature models and convective heat transfer in porous media. They rarely describe transient, reservoir temperature changes during conventional production. Here, fluid expansion effects are superimposed on heat conduction to the surrounding layers coupled to a transient inflow regime. This paper will discuss why these already available solutions to related problems cannot be employed directly. It will propose a comprehensive model of the reservoir temperature around a horizontal, liquid production well. The above task is treated as a multidimensional, heat transfer problem. Both an asymptotic, early time, solution in real space and an accurate, integral solution can be derived. The asymptotic, analytical solution has been tested against a scientific, simulation softwares. The applicability and practical importance of the solution will be discussed. It forms the basis of a Temperature Transient Analysis tool for horizontal, oil producing wells where both the flow properties of the reservoir fluid and the inflow distribution are to be evaluated.

We consider flows possessing a combustion front when a gaseous oxidizer (air) is injected into the porous medium, a rock cylinder thermally insulated preventing lateral heat losses and filled with light or medium viscosity oil. When oxygen reacts with hydrocarbons at low temperatures, a series of reactions occur that will convert a part of hydrocarbons to oxygenated hydrocarbons and gaseous product (water, carbon dioxide etc.). The oxygenated hydrocarbons are compounds like ketones, alcohols, aldehydes and acids. This process is termed low temperature oxydation (LTO). Indeed, LTO only involves some 25 % of the possible sites that can react with oxygen. Therefore also the reaction heat per volume of fuel can at most be 25 % of full hydrocarbon combustion and consequently the temperature in the LTO reaction zone is very low. Upstream of the LTO reaction zone evaporization occurs. We formulate conservation laws for liquid oil, gaseous oil, oxygen and inert gas components (combustion products and nitrogen) that includes the reaction terms. Moreover we give the energy conservation equation. We give an analytical solution to the equations. It turns out that the solution consists, from upstream to downstream, of a thermal wave, an LTO wave that combines oxidation and evaporation and a Buckley-Leverett saturation wave. It is shown that in the solution a major role is played by the existence of a resonance line at which the derivatives of the oxygen and evaporated oil flux versus the oil saturation vanish. It means that the derivative of the oil flux versus saturation is positive upstream of the resonance line and negative downstream of the resonance line. The complete solution is described for typical parameters of LTO oil combustion. Computations show that only a small part of the oil vaporizes or reacts. The oil velocity is close to the LTO speed. Thus the LTO wave represents a mechanism of almost complete oil displacement by means of the temperature increase in the LTO wave, which leads to a decrease of oil viscosity and increase of the gas flux in the wave.