Buckley-Leverett Research Articles

Stochastic porous media modeling and high-resolution schemes for numerical simulation of subsurface immiscible fluid flow transport

This paper proposes stochastic petroleum porous media modeling for immiscible fluid flow simulation using Dykstra–Parson coefficient (VDP) and autocorrelation lengths to generate 2D stochastic permeability values which were also used to generate porosity fields through a linear interpolation technique based on Carman–Kozeny equation. The proposed method of permeability field generation in this study was compared to turning bands method (TBM) and uniform sampling randomization method (USRM). On the other hand, many studies have also reported that, upstream mobility weighting schemes, commonly used in conventional numerical reservoir simulators do not accurately capture immiscible displacement shocks and discontinuities through stochastically generated porous media. This can be attributed to high level of numerical smearing in first-order schemes, oftentimes misinterpreted as subsurface geological features. Therefore, this work employs high-resolution schemes of SUPERBEE flux limiter, weighted essentially non-oscillatory scheme (WENO), and monotone upstream-centered schemes for conservation laws (MUSCL) to accurately capture immiscible fluid flow transport in stochastic porous media. The high-order schemes results match well with Buckley Leverett (BL) analytical solution without any non-oscillatory solutions. The governing fluid flow equations were solved numerically using simultaneous solution (SS) technique, sequential solution (SEQ) technique and iterative implicit pressure and explicit saturation (IMPES) technique which produce acceptable numerical stability and convergence rate. A comparative and numerical examples study of flow transport through the proposed method, TBM and USRM permeability fields revealed detailed subsurface instabilities with their corresponding ultimate recovery factors. Also, the impact of autocorrelation lengths on immiscible fluid flow transport were analyzed and quantified. A finite number of lines used in the TBM resulted into visual artifact banding phenomenon unlike the proposed method and USRM. In all, the proposed permeability and porosity fields generation coupled with the numerical simulator developed will aid in developing efficient mobility control schemes to improve on poor volumetric sweep efficiency in porous media.

Doctrines vs. Realities in Reservoir Engineering

This article, written by Special Publications Editor Adam Wilson, contains highlights of paper SPE 185633, “Doctrines and Realities in Reservoir Engineering,” by Euthymios Vittoratos, Petroprognostica, and Anthony Kovscek, SPE, Stanford University, prepared for the 2017 SPE Western Regional Meeting, Bakersfield, California, USA, 23–27 April. The paper has not been peer reviewed. The appreciation of empirical realities permits the improvement of commercial depletion planning and enables a greater number of projects. The concepts articulated here are applicable in almost all heavy-oil reservoirs and may be applicable also for lighter oils that possess oil chemistry typical of heavier oils (e.g., high acid content). This paper reviews the evidence for and against three doctrines in current use to develop depletion plans: (1) optimal recovery is obtained using a voidage-replacement ratio (VRR) of unity, (2) the Buckley-Leverett (BL) formulation (phase slippage) applies uniformly to heavier oils, and (3) viscous fingering dominates unstable multiphase flows. VRR=1 Doctrine: Optimal Waterflood Response Occurs With VRR=1 Reality: Periods of VRR<1 Increase Oil Recovery Though initially viewed with some suspicion, waterflooding has demonstrated multifold greater recovery than primary depletion for many reservoirs. It is, by far, the most important oil-recovery process. One visualization of the process is that of a leaky, deformable water piston displacing a fraction of the oil to the production wells. Maximum recovery occurs when the voidage created by the produced fluids (oil plus water) at reservoir conditions (pressure and temperature) is replaced fully by an equal volume of the injected water also at reservoir conditions (i.e., the waterflood VRR=1). In this way, reservoir pressure is maintained. The waterflood initial production is pure oil until water breakthrough. Subsequently, water is produced along with the oil, resulting in an increasing water/oil ratio (WOR). This simple visualization is intuitively attractive and has dominated the thinking of the industry. This “leaky, deformable piston” visualization is simplistic and neglects or sidelines complex physical and chemical properties of the reservoir. If the waterflood is operated below the oil bubblepoint pressure through a continued VRR<1, a gaseous phase forms, creating three-phase-flow conditions and more-complex displacement possibilities. Furthermore, reservoir oils can have widely different physical and chemical properties, varying by many orders of magnitude in viscosity and substantially in interfacial activity as described by metrics such as the total acid number. Despite the full recognition of these complications, over time, the now al-most universal practice developed of operating waterfloods with VRR=1. Analytical and Numerical VRR Simulations. VRR<1 activates many recovery mechanisms. Some mechanisms can be enumerated and examined, and others may exist that are yet to be identified. What is the possible upside of operating a waterflood with VRR<1? Making a rough estimate is not too difficult if we limit the number of mechanisms. With VRR<1, two dominant recovery mechanisms exist, solution-gas drive (VRR=0) and water displacement (VRR=1). As a first approximation, at an intermediary VRR<1, the two mechanisms may be additive: Oil Recovery (at Optimal VRR)=Oil Recovery by Pure Waterflood (VRR=1)+Pure Solution-Gas Drive (VRR=0).

Significance of Geochemistry in Single-Well Chemical-Tracer Tests by Coupling a Multiphase-Flow Simulator to the Geochemical Package

Summary Single-well chemical tracer (SWCT) is the most commonly used field method to determine oil or water saturation in one-well enhanced-oil-recovery (EOR) pilots. Because hydrolysis of an ester, which is the main part of the method, leads to forming acid as well as alcohol, the equilibrium state of the reservoir is disturbed, and thus the pH changes. It is generally accepted that the hydrolysis-reaction rate is mainly dependent on the pH and temperature. Therefore, it is required to know the extent to which this dependency might affect the shape of the product-tracer profiles and the numerical interpretation of the field-test data for computing residual oil saturation (Sor). In this study, this notion has been investigated by coupling a multiphase-flow simulator to the geochemistry package PHREEQC (Parkhust and Appelo 2013). In this study, the PHREEQC geochemical simulator has been used to illustrate the extent to which different parameters might affect the pH variation during the test. The PHREEQC database has been modified to take the ester-hydrolysis reaction into account by adding the ester, alcohol, and acid-product species. The hydrolysis-reaction mechanisms of ester have also been programmed to account for the dependency of the hydrolysis reaction on the pH. Also, because ester partitions into the oil phase and travels behind the water phase (i.e., Darcy velocity), performing two-phase flow would be necessary to highlight the significance of the pH dependency of the hydrolysis-reaction rate on the tracer profiles. For doing that, a multiphase Buckley-Leverett (BL) flow simulation is coupled with IPhreeqc, which is an open-source module of the PHREEQC geochemical package. Then, a California Turbidite SWCT test has been re-evaluated to verify the approach. At the end, the geochemistry of a reservoir with an almost weak resistance (high temperature and weak buffer capacity) against pH variation in the SWCT test has been studied using the geochemical-based approach. The results show that the variation of the hydrolysis rate with pH could affect mainly the tail edge of the predicted tracer profiles, and it could marginally affect the apex of the profiles; however, it might affect the interpreted value of the Sor measurement as the resistance against pH variation becomes weaker. In these conditions, adapting the SWCT-test designs (i.e., shut-in time and injecting lower concentration of ester) could diminish the pH variation. The pH dependency of the hydrolysis-reaction rate is recommended for the numerical interpretation of the field SWCT-test data. The results of this study can be used to minimize the uncertainties of the SWCT tests and to improve the reliability of the Sor measurements.

Adjoint analysis of Buckley-Leverett and two-phase flow equations

This paper analyzes the adjoint equations and boundary conditions for porous media flow models, specifically the Buckley-Leverett equation, and the compressible two-phase flow equations in mass conservation form. An adjoint analysis of a general scalar hyperbolic conservation law whose primal solutions include a shock jump is initially presented, and the results are later specialized to the Buckley-Leverett equation. The non-convexity of the Buckley-Leverett flux function results in adjoint characteristics that are parallel to the shock front upstream of the shock and emerge from the shock front downstream of the shock. Thus, in contrast to the behavior of Burgers’ equation where the adjoint is continuous at a shock, the Buckley-Leverett adjoint, in general, contains a discontinuous jump across the shock. Discrete adjoint solutions from space-time discontinuous Galerkin finite element approximations of the Buckley-Leverett equation are shown to be consistent with the derived closed-form analytical solutions. Furthermore, a general result relating the adjoint equations for different (though equivalent) primal equations is used to relate the two-phase flow adjoints to the Buckley-Leverett adjoint. Adjoint solutions from space-time discontinuous Galerkin finite element approximations of the two-phase flow equations are observed to obey this relationship.

Group classification applications for analysis of discrete models of flow in porous media

This work is about examples of applications of group classifications for continuous (differential equations) and discrete models of flow in porous media. Difference schemes and discrete dynamical systems are considered as discrete models. Group classifications has been carried out on the basis of continuous groups (Lie groups) of transformations in n-dimensional Euclidean spaces. The example of explicit difference scheme for Buckley-Leverett equation has been demonstrated where the correspondence to one of classes of the obtained earlier classifications of discrete dynamical system has been shown. The example of fractal capillary network on the basis of Sierpinski triangle has been demonstrated for the group classification of discrete dynamical systems. The approach of continuous symmetry absence or presence proof has been demonstrated using classifications of discrete dynamical systems. Difference schemes with continuous symmetry groups for partial differential equations of parabolic type, which correspond to the equation of gas flow in porous media and Rapoport-Leas equation (generalized Buckley-Leverett equation), are considered. Examples of difference schemes for these cases are given and related to classes of the obtained earlier classifications. Applications of continuous symmetry groups for numerical solutions generations are demonstrated for these cases.

Application of the CABARET scheme for calculation of discontinuous solutions of the scalar conservation law with nonconvex flux

The monotonicity of the CABARET scheme approximating the quasi-linear scalar conservation law with nonconvex flux is analyzed. The conditions of monotonicity of this scheme are obtained in the case when the speed of propagation of characteristics of the approximated equation is limited and nonnegative. As an example, the test calculations by the CABARET scheme of discontinuous solutions of the Buckley–Leverett equation are presented.

A Moving Mesh Finite Difference Method for Non-Monotone Solutions of Non-Equilibrium Equations in Porous Media

AbstractAn adaptive moving mesh finite difference method is presented to solve two types of equations with dynamic capillary pressure effect in porous media. One is the non-equilibrium Richards Equation and the other is the modified Buckley-Leverett equation. The governing equations are discretized with an adaptive moving mesh finite difference method in the space direction and an implicit-explicit method in the time direction. In order to obtain high quality meshes, an adaptive monitor function with directional control is applied to redistribute the mesh grid in every time step, then a diffusive mechanism is used to smooth the monitor function. The behaviors of the central difference flux, the standard local Lax-Friedrich flux and the local Lax-Friedrich flux with reconstruction are investigated by solving a 1D modified Buckley-Leverett equation. With the moving mesh technique, good mesh quality and high numerical accuracy are obtained. A collection of one-dimensional and two-dimensional numerical experiments is presented to demonstrate the accuracy and effectiveness of the proposed method.

Generalizing the modified Buckley–Leverett equation with TCAT capillary pressure

The Buckley–Leverett partial differential equation has long been used to model two-phase flow in porous media. In recent years, the PDE has been modified to include a rate-dependent capillary pressure constitutive equation, known as dynamic capillary pressure. Previous traveling wave analysis of the modified Buckley–Leverett equation uncovered non-classical solutions involving undercompressive shocks. More recently, thermodynamically constrained averaging theory (TCAT) has generalized the capillary pressure equation by including additional dependence on fluid properties. In this paper, the model and traveling wave analysis are updated to incorporate TCAT capillary pressure as a generalization of dynamic capillary pressure. Solutions of the corresponding Riemann problem are similar to previous results except in the physically relevant situation in which both phases are pure fluids. The results presented here shed new light on the nature of the interface between one pure fluid displacing another pure fluid, in accordance with TCAT.

Totally Volume Integral of Fluxes for Discontinuous Galerkin Method (TVI-DG) I-Unsteady Scalar One Dimensional Conservation Laws

The volume integral of Riemann flux in the discontinuous Galerkin (DG) method is introduced in this paper. The boundaries integrals of the fluxes (Riemann flux) are transformed into volume integral. The new family of DG method is accomplished by applying divergence theorem to the boundaries integrals of the flux. Therefore, the (DG) method is independent of the boundaries integrals of fluxes (Riemann flux) at the cell (element) boundaries as in classical (DG) methods. The modified streamline upwind Petrov-Galerkin method is used to capture the oscillation of unphysical flow for shocked flow problems. The numerical results of applying totally volume integral discontinuous Galerkin method (TVI-DG) are presented to unsteady scalar hyperbolic equations (linear convection equation, inviscid Burger's equation and Buckley-Leverett equation) for one dimensional case. The numerical finding of this scheme is very accurate as compared with other high order schemes as the weighted compact finite difference method WCOMP.

Simulation of incompressible two-phase flow in porous media with large timesteps

Multiphase flow in porous media occurs in several disciplines including petroleum reservoir engineering, petroleum systems' analysis, and CO2 sequestration. While simulations often use a fully implicit discretization to increase the time step size, restrictions on the time step often exist due to non-convergence of the nonlinear solver (e.g. Newton's method). Here this problem is addressed for the Buckley–Leverett equations, which model incompressible, immiscible, two-phase flow with no capillary potential. The equations are recast as a gradient flow using the phase-field method, and a convex energy splitting scheme is applied to enable large timesteps, even for high degrees of heterogeneity in permeability and viscosity. By using the phase-field formulation as a homotopy map, the underlying hyperbolic flow equations can be solved with large timesteps. For a heterogeneous test problem, the new homotopy method allows the timestep to be increased by more than six orders of magnitude relative to the unmodified equations while maintaining convergence.

Adaptive WENO Methods Based on Radial Basis Function Reconstruction

We explore the use of radial basis functions (RBF) in the weighted essentially non-oscillatory (WENO) reconstruction process used to solve hyperbolic conservation laws, resulting in a numerical method of arbitrarily high order to solve problems with discontinuous solutions. Thanks to the mesh-less property of the RBFs, the method is suitable for non-uniform grids and mesh adaptation. We focus on multiquadric radial basis functions and propose a simple strategy to choose the shape parameter to control the balance between achievable accuracy and the numerical stability. We also develop an original smoothness indicator which is independent of the RBF for the WENO reconstruction step. Moreover, we introduce type I and type II RBF-WENO methods by computing specific linear weights. The RBF-WENO method is used to solve linear and nonlinear problems for both scalar and systems of conservation laws, including Burgers equation, the Buckley–Leverett equation, and the Euler equations. Numerical results confirm the performance of the proposed method. We finally consider an effective conservative adaptive algorithm that captures moving shocks and rapidly varying solutions well. Numerical results on moving grids are presented for both Burgers equation and the more complex Euler equations.

Solution of the equations for one-dimensional, two-phase filtration in porous medium with heat exchange by geometric methods

Buckley-Leverett equations describe non viscous, immiscible, two-phase filtration, which is often of interest in modelling of oil production. For many parameters and initial conditions, the solutions of these equations exhibit non-smooth behaviour, namely discontinuities in form of shock waves. In this paper we propose an asymptotic method of obtaining the solution and means of resolving the singularities and shock waves. This method is implemented as Maple code.

Solution of the equations for one-dimensional, two-phase, immiscible flow by geometric methods

Buckley–Leverett equations describe non viscous, immiscible, two-phase filtration, which is often of interest in modelling of oil production. For many parameters and initial conditions, the solutions of these equations exhibit non-smooth behaviour, namely discontinuities in form of shock waves. In this paper we obtain a novel method for the solution of Buckley–Leverett equations, which is based on geometry of differential equations. This method is fast, accurate, stable, and describes non-smooth phenomena. The main idea of the method is that classic discontinuous solutions correspond to the continuous surfaces in the space of jets - the so-called multi-valued solutions (Bocharov et al., Symmetries and conservation laws for differential equations of mathematical physics. American Mathematical Society, Providence, 1998). A mapping of multi-valued solutions from the jet space onto the plane of the independent variables is constructed. This mapping is not one-to-one, and its singular points form a curve on the plane of the independent variables, which is called the caustic. The real shock occurs at the points close to the caustic and is determined by the Rankine–Hugoniot conditions.

Computing numerical solutions of the pseudo-parabolic Buckley–Leverett equation with dynamic capillary pressure

We present numerical approaches for solving a pseudo-parabolic partial differential equation, which models incompressible two phase flow in porous media taking into account dynamic effects in the capillary pressure. First, we briefly discuss two numerical schemes based on the operator splitting technique. Our numerical experiments show that the standard splitting, widely used to solve parabolic problems, may fail when applied to pseudo-parabolic models. As an illustration, we give an example for this case. So we present an operator splitting scheme based on a dispersive-like character that obtains correct numerical solutions. Then, we discuss an unsplit efficient numerical modelling, locally conservative by construction. This framework is based on a fully coupled space–time mixed hybrid finite element/volume discretization approach in order to account for the delicate local nonlinear balance between the numerical approximations of the hyperbolic flux and the pseudo-parabolic term, but linked to a natural dispersive-like character of the full pseudo-parabolic equation. We compare our numerical results with approximate solutions constructed with methods recently introduced in the specialized literature, in order to establish that we are computing the expected qualitative behaviour of the solutions.

Travelling wave solution of the Buckley–Leverett equation

A two-dimensional Buckley–Leverett system governing motion of two-phase flow is considered. Travelling-wave solutions for these equations are found. Wavefronts of these solutions may be circles, lines and parabolae. Values of pressure and saturation on the wave fronts are found.

Stochastic Galerkin framework with locally reduced bases for nonlinear two-phase transport in heterogeneous formations

The generalized polynomial chaos method with multiwavelet basis functions is applied to the Buckley–Leverett equation. We consider a spatially homogeneous domain modeled as a random field. The problem is projected onto stochastic basis functions which yields an extended system of partial differential equations. Analysis and numerical methods leading to reduced computational cost are presented for the extended system of equations.The accurate representation of the evolution of a discontinuous stochastic solution over time requires a large number of stochastic basis functions. Adaptivity of the stochastic basis to reduce computational cost is challenging in the stochastic Galerkin setting since the change of basis affects the system matrix itself. To achieve adaptivity without adding overhead by rewriting the entire system of equations for every grid cell, we devise a basis reduction method that distinguishes between locally significant and insignificant modes without changing the actual system matrices.Results are presented for problems in one and two spatial dimensions, with varying number of stochastic dimensions. We show how to obtain stochastic velocity fields from realistic permeability fields and demonstrate the performance of the stochastic Galerkin method with local basis reduction. The system of conservation laws is discretized with a finite volume method and we demonstrate numerical convergence to the reference solution obtained through Monte Carlo sampling.

Control of displacement front in a model of immiscible two-phase flow in porous media

For the Buckley–Leverett equation describing the flow of two immiscible fluids in porous media, an exact parametric representation of the solution is constructed with the help of the Backlund transformation. As a result, the advance of the displacement front can be controlled to a high degree of accuracy. The method is illustrated using an example of a typical oil well with actual parameters.

Explicit strong stability preserving multistep Runge–Kutta methods

High-order spatial discretizations with strong stability properties (such as monotonicity) are desirable for the solution of hyperbolic PDEs. Methods may be compared in terms of the strong stability preserving (SSP) time-step. We prove an upper bound on the SSP coefficient of explicit multistep Runge–Kutta methods of order two and above. Order conditions and monotonicity conditions for such methods are worked out in terms of the method coefficients. Numerical optimization is used to find optimized explicit methods of up to five steps, eight stages, and tenth order. These methods are tested on the advection and Buckley-Leverett equations, and the results for the observed total variation diminishing and positivity preserving time-step are presented.

Modeling process of embolization arteriovenous malformation on the basis of two-phase filtration model

Arteriovenous malformation is a chaotic disordered interlacement of very small diameter vessels, performing reset of blood from the artery into the vein. In this regard it can be adequately modeled using porous medium. In this model process of embolization described as penetration of non-adhesive substance ONYX into the porous medium, filled with blood, both of these fluids are not mixed with each other. In one-dimensional approximation such processes are well described by Buckley-Leverett equation. In this paper Buckley-Leverett equation is solved numerically by using a new modification of Cabaret scheme. The results of numerical modeling process of embolization of AVM are shown.

On the generalized Buckley-Leverett equation

In this paper we study the generalized Buckley-Leverett equation with nonlocal regularizing terms. One of these regularizing terms is diffusive, while the other one is conservative. We prove that if the regularizing terms have order higher than one (combined), there exists a global strong solution for arbitrarily large initial data. In the case where the regularizing terms have combined order one, we prove the global existence of solution under some size restriction for the initial data. Moreover, in the case where the conservative regularizing term vanishes, regardless of the order of the diffusion and under a certain hypothesis on the initial data, we also prove the global existence of a strong solution, and we obtain some new entropy balances. Finally, we provide numerics suggesting that, if the order of the diffusion is 0 < α < 1, a finite time blow up of the solution is possible.