Scholarly article on topic 'Simulation-based Determination of Relative Permeability in Laminated Rocks'

Simulation-based Determination of Relative Permeability in Laminated Rocks Academic research paper on "Earth and related environmental sciences"

Share paper
Academic journal
Energy Procedia
{"Relative permeability determination" / "Laminated rocks" / "Numerical simulation"}

Abstract of research paper on Earth and related environmental sciences, author of scientific article — Mohammad H. Sedaghat, Kirill Gerke, Siroos Azizmohammadi, Stephan K. Matthai

Abstract Reservoir simulation using the extended Darcy's law approach requires relative permeability curves derived either via analytic saturation functions (Corey models etc.) or from special core analysis (SCAL). Since such experimental exploration of the space of influential parameters (pore geometry and wettability) is costly and time consuming, establishing ways to extract ensemble relative permeability from numerical simulation, k ri , over the entire range of water saturation is highly desirable. In this work, a Steady State Saturation Variation (SSSV) technique is proposed. It computes ensemble k ri (s w ) for layered rocks when both capillary and viscous forces are strong.

Academic research paper on topic "Simulation-based Determination of Relative Permeability in Laminated Rocks"


Available online at


Energy Procedia 97 (2016) 433 - 439

European Geosciences Union General Assembly 2016, EGU Division Energy, Resources & Environment, ERE

Simulation-Based Determination of Relative Permeability in

Laminated Rocks

Mohammad H. Sedaghata*, Kirill Gerkea, Siroos Azizmohammadib, Stephan K. Matthaia

aDepartment of Infrastructure Engineering, The University of Melbourne, Grattan Street 175, Parkville 3010, Melbourne, Australia bChair of Reservoir Engineering, Montan University of Leoben, Parkstrasse 27, 8700 Leoben, Austria


Reservoir simulation using the extended Darcy's law approach requires relative permeability curves derived either via analytic saturation functions (Corey models etc.) or from special core analysis (SCAL). Since such experimental exploration of the space of influential parameters (pore geometry and wettability) is costly and time consuming, establishing ways to extract ensemble relative permeability from numerical simulation, kn, over the entire range of water saturation is highly desirable. In this work, a Steady State Saturation Variation (SSSV) technique is proposed. It computes ensemble kri(sw) for layered rocks when both capillary and viscous forces are strong.

© 2016 The Authors.Publishedby ElsevierLtd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.Org/licenses/by-nc-nd/4.0/).

Peer-review under responsibility of the organizing committee of the General Assembly of the European Geosciences Union (EGU) Keywords: Relative permeability determination; Laminated rocks; Numerical simulation;

1. Introduction

The relative permeability kri(sw) quantifies the flow of a fluid phase i through a medium in presence of another fluid phase [1]. In a water-wet medium at the capillary limit, water is restricted to the small pores and oil therefore flows with less resistance. Hence, krw has a smaller value than kro as compared with an oil-wet medium where this relationship is reversed. It follows that the shape of relative permeability curves is a strong indicator of wettability. Capillary pressure (Pc) is usually defined as the difference between wetting and non-wetting phase pressure. It is influenced by the interfacial tension, contact angle, and the mean radius of curvature of non-wetting phase menisci

"Corresponding author. Tel.: +61401687786 ; fax.: +61 3 834 46084 Email:

1876-6102 © 2016 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (

Peer-review under responsibility of the organizing committee of the General Assembly of the European Geosciences Union (EGU) doi: 10.1016/j.egypro.2016.10.041

constrained by the diameter of pores and throats occupied as a function of water saturation, Sw [2]. As the contact angle changes from water-wet to oil-wet conditions, capillary pressure assumes different values in water-wet and oil-wet media: where strongly water-wet Pc is positive and declines with water saturation; where strongly oil-wet Pc is negative and its absolute value gets greater as Sw increases. So, both kri and Pc are key indicators of wettability [3]. Layers of a distinct mineralogy may have different a wettability. Precipitation of minerals in the void space of rocks such as chlorite or feldspar also can create different wettability characteristics. Two adjacent rocks may thus assume a different wettability and rocks below and above the water-oil contact in an oil reservoir, source rock or along hydrocarbon pathways may experience a change of their primary wettability. Wettability may also vary across compositional laminations (mm-cm scale layering) in porous sedimentary rocks.

Relative permeability is fundamental to reservoir engineering, and accurate reservoir simulation critically depends on the availability of kri data. Thus, reservoir characterisation typically involves Special Core Analysis (SCAL) based on steady or unsteady (core-flooding) laboratory experiments. However, for larger cores, such experiments are known to be time-consuming and a determination of kri over the full range of water saturations from Swi to Sor may not be possible. In the classic steady state method, both ends of the core sample are kept at a constant pressure while the two phases are flowing simultaneously at a rate that is measured [5]. In the unsteady approach, one end of the core sample is held at a constant pressure while fluid is injected at the other at a constant rate [4]. In both cases kri is inferred from phase outflow rates. If one replicates these setups in a numeric model, kri measurements cannot be made before breakthrough of the injected phase to the other end occurs. Thus, to compute kri for larger multi-layer systems where breakthrough may take a long time, a different simulation analysis method is required. For stratified reservoirs, Hearn [6] recognised this issue and calculated pseudo relative permeability curves. His model is based on arithmetic averaging. Although it scales up relative permeability in thickly layered rocks, it does not work for thinly laminated systems because it neglects capillary transfer between the layers precluding consideration of related effects arising from wettability variations.

In this work, we introduce a novel way to extract ensemble relative permeability from numerical simulations. First, it is shown that this approach works significantly better than the direct numeric replication of laboratory methods and, in contrast with the lab, it permits to determine ensemble relative permeability curves over the entire range of water saturation. Second, the method is applied to a model of laminated rock with mm-scale laminations as characterized in the lab. Wettability is prescribed as discontinuous across layer boundaries to maximize cross flow. The evolving flow behaviour near the material interfaces is investigated and discussed.

2. Methodology

Darcy's law for multiple fluid phases in the absence of gravity is [1]:

n = (1)

where k, P, v, ^ and kr are respectively absolute permeability, pressure, velocity, viscosity, and relative permeability. Index "i" could be water, oil, gas or any flowing fluid. While ^ and k are known, kri at every pressure gradient is easily calculated after is determined.

The following experimental and simulation approaches are normally used to obtain relative permeability: (1) Experiment-based Unsteady State (US): Water injection rate and production pressure are set constant and pressure is measured at the inlet. Water flooding is run once and the water cut, water saturation and total mobility are calculated using material balance and Darcy law [4]. This method is readily replicated numerically for highly heterogeneous models where the displacing front is not uniform and water breakthrough occurs early. This includes discrete fracture and matrix models with high fracture-matrix flux ratio promoting an early breakthrough time [7]. However, for large (viscous-force) dominated models with a compact displacement front that is not levelled out by capillary spreading, calculation of the kri of the injected phase prior to the propagation of its shock front to the outlet boundary is not possible. This should also apply to data post processing for large-scale physical experiments since JBN (Johnson-Bossler-Naumann) or similar calculations require the determination of water saturations at the sample outlet [8]. (2) Experiment-based Steady State (SS): Pressure at both boundaries is fixed and input rates for

each phase are measured. The test is run at a constant pressure gradient. The two phases are injected simultaneously but at different proportions. The combined injection rates have to be adjusted at each time step to maintain the constant pressure gradient. Measuring injection rates for each phase, kri is calculated at the given sample saturation. Saturations are determined by X-Ray scanning. This method is practical for all types of sample geometries, however, it is impractical and time consuming because phase velocities have to be adjusted after each saturation increment to maintain the constant pressure gradient [5]. (3) Simulation-based Unsteady State (US): Similar to experimental US water flooding is performed, however, volume averaging of result properties is performed over the entire model to get pressure gradients and phase velocities [9]. (4) Simulation-based Steady State (SS): It is applicable when the sample state is near the capillary limit (CL) in which capillary pressure gradient is zero (VPC = 0) or viscous limit (VL) in which fractional flow gradient is zero (Vfw = 0). This implies that saturation must be uniformly distributed over the grid [10]. Since these conditions are never fulfilled in numerical models of laminated rock with wettability contrasts among layers, SS is not suitable for the analysis presented here.

Thus, a Steady State Saturation Variation (SSSV) method is proposed to obtain kri values over the entire saturation range from Swi to Sor: first, the model is initialized to a constant value of saturation. Since the layers exchange fluids along their interfaces by capillary forces, the ensuing saturation distributions alter ensemble kri over time. There is a characteristic equilibration time when the saturation pattern stabilizes. However, in contrast to the simulation-based steady state method VPC and V/w remain non zero in the entire region. To measure kri at this saturation (equilibrium state), the initial saturation also needs to be applied as a boundary condition for inflow on the injection boundary. Then the fractional inflow rate of each phase can be measured for the given VP. As the saturation distribution is frozen in the whole model, velocity of each phase is constant, and material balance gives the steady state relative permeability at the corresponding Sw(t). This procedure is similar to that applied in physical experiments. The horizontal velocity may vary vertically along the inflow boundary when the model/sample is stratified horizontally and the inflow cells are located in different material domains. In this case, the inflow velocity is averaged along the inflow boundary. In the FECFVM simulator used in this study [13], this is accomplished by finite volume surface integration of piecewise constant finite element velocities at the inflow boundary held at constant pressure.

To verify the SSSV method, it was first applied to a homogenous 10 x 10 meters isotropic model. Since there is no discontinuity in Pc, the saturation distribution ensuing from different assigned initial water saturation values is always uniform, equilibrium is established at the first time step, and there is no necessity for an initialisation computation. Also, the determined kri is insensitive to the applied pressure gradient. To conduct the numeric experiment a rectangular domain was meshed and assigned with a single set of petrophysical properties. General Corey saturation functions were used. For the second more realistic verification experiment, a 10-cm laminated rock sample was chosen as characterised by Crocker Data Processing Pty Ltd research [11]. The petrophysical properties of this sample have been constrained on the mm-scale. Data relevant for the simulation was extracted from the corresponding core log and assigned to the discretized model. As shown in Figure 1, porosity is constant and vertical and horizontal permeability are equal. Although the lithology of the layers is not the same and silt and clay are concentrated in the upper layer while coarse sand dominates the lower one, the horizontal permeability of these two layers is very close. This is why these two narrow 5x20 cm layers with same porosity and close and isotropic permeability but possible different SCAL data were selected. Using the correlation proposed by Huet [12] (Equations 2 and 3), entry pressure, Pd, and Corey parameter, y, were derived from the known porosity and permeability taken from the log (Fig. 1.). Results are shown in Table 1. Corey saturation functions are shown in Equations 4 to 7. The upper layer was assumed to be oil-wet. So, since its contact angle is assumed 180°, its Pd is negative.

Pd = 751.3360.0°-8469.^-° S166. (1 -5wi)0 0489 (2)

Y = 0.0 0 0 84.0-10485 ^0.5498^ (1 _ swi)~22790. P°"39 (3)

Where effective water saturation (Swe) is defined as:

k - S Y

"■rw Jwe

k-ro = (1 _ $we)2- (1 _ )

Pc - Pd-Sje

(6) (7)

The simulations benefit from the capability of the Finite Element Centered Finite Volume (FECFV) to model capillary transfer processes across layer boundaries accurately [13]. Time marching was performed using an IMPES (sequential implicit solution of pressure equation followed by explicit solution of the saturation transport equation) approach. For the simulations presented here in an incompressible two-phase flow formulation was used, eliminating the possibility for the emergence of pressure transients.

Table 1: Properties of the layers.

Parameters Upper Layer Lower Layer

Average Porosity 0.112 0.112

Average Permeability (mD) 737 868

Lithology Fine Sand, Silt, Clay Coarse Sand, Fine Sand

Pd (Pa) -26482.87 24477.15

Y 1.99 1.54

Swi 0.1 0.2

Sor 0.2 0.1

Contact angle (0) 180° 0°

3. Results and Discussions

In Figure 2, the kri results obtained with the SSSV approach are contrasted with the other approaches (experimental and simulation based unsteady-state approach). There is an acceptable match between the assigned kri values and the computed ones. The experimental-based unsteady state (US), JBN [8] modification and simulation-based unsteady state methods failed to predict kri before front breakthrough. Only the SSSV results fit the expected outcome exactly. Thus, the large homogeneous isotropic test outcome verifies the SSSV method while highlighting the limitations of the classical approaches at a sample scale where capillary pressure equilibrium may never be attained in the presence of simultaneous flow through the block. The numeric experiment shows, however, that all the approaches work well after breakthrough, i.e. when fully saturated conditions have been established.

■ Kro - SSSV

■ Krw - SSSV

Kro-US- Experimental Way Krw-US- Experimental Way

> Kro-JBN - Experimental Way

> Krw - JBN - Experimental Way Kro-US-Simulation Way Krw-US-Simulation Way J Kro - Brooks Corey Krw - Brooks Corey

0 0.2 0.4 0.6 0.!

0.2 0.4 0.6 0.8

Fig. 2. Relative permeability curves computed using various methods on the (a) homogenous and (b) layered sample (AP=1000 Pa).

Figure 2b is a plot of ensemble relative permeability curves computed using the SSSV approach for the layered system comparing kn values computed on the equilibrated model (solid lines) with those obtained directly from the initially prescribed saturation values of the layers (stippled lines) for a given flow rate through the layers. The equilibrium ensemble relative permeability curve plots exactly between the disequilibrium (initial state) water-wet and oil-wet pair of curves. It is not impacted by water breakthrough to the model boundary and gives kn values for all saturation states. Moreover, water saturation in the kri crossover point (the point krw = kra) would normally be taken as an indication of an neutral wettability between the ones prescribed to the model. While these results appear logical in hindsight, they highlight an interpretational issue similar to the one which arises with Amott indices in the case of fractional wettability: there is a composite effect that cannot be directly untangled from the volume averaged results. In addition, the position of the solid curve is injection rate dependent. The behavior of SSSV ensemble curves at different wettability conditions and the role of layer sorting and assigned kri models on the dynamic ensemble relative permeability curves as obtained by SSSV approach is subject of another article which addresses the effects of gravity as well.

It should also be noted that the degree of smoothness of the curves obtained from the numerical models is a function of the number of equilibration-measurement steps. Performing more steps increases the computational effort. This is similar to physical experiments.

Figures 3 displays saturation, velocity and total mobility distributions when the layered model reaches a steady state. As is obvious from Fig. 3a, saturation varies across the model and the wetting phase is accumulating directly below the layer interface. Fig. 3b reveals that the change in wettability across the interface gives rise to vigorous flow along it because the total mobility is maximized. The wetting-non-wetting layer interface behaves like a high-permeability zone, yet this is a total mobility effect due to non-equal changes of wetting and non-wetting phase relative permeability values (Fig. 3c). Over a range of layer-parallel flow velocities this total mobility effect is sufficiently prominent to warrant consideration in the upscaling of saturation functions in layered porous media. Since the saturation gradients causing this effect are flow-velocity dependent (Fig. 3d) our observations challenge classic volume-based upscaling approaches that have no provisions for such a rate-dependent behavior.

(c) (d)

Fig. 3. Steady-state saturation (a) flow velocity magnitude (b) total mobility (c) and total mobility when AP is four times higher (d)

distributions in the layered model [14].

4. Conclusions

Standard experimental relative-permeability determination methods were replicated in numeric experiments and then applied on the 10-m scale where saturation patterns are no longer dominated by capillary forces. Then, curves were determined numerically for models of homogenous uniform isotropic and layered sedimentary rocks. The results highlight issues arising from non-uniform saturation distributions that require a different experimental setup.

A method is proposed that involves the computation of saturation distributions that reflect a dynamic equilibrium between viscous and capillary forces. With this input dynamic relative permeability curves are obtained that are dynamic, capturing inter-layer phase exchanges. This rate-dependent behaviour cannot be captured by conventional relative permeability upscaling approaches. In the capillary limit, the proposed analysis method gives results identical to standard methods. When individual laminations of a layered system have different wettabilities, interface saturation configurations can arise that strongly enhance total mobility. This makes such interfaces behave like high-permeability streaks.


[1] Muskat M, Wyckoff RD, Botset HG, Meres MW. Flow of gas-liquid mixtures through sands. Trans Amer Inst. Mining and Metallurgical Engineers 1937;123:69-98.

[2] Muskat M. Calculation of initial fluid distributions in oil reservoirs. Trans of AIME 1949;179:119-127.

[3] Anderson WG. The effects of wettability on relative permeability. J Pet Tech 1987;1453-1987.

[4] Lake LW. Enhanced oil recovery. New Jersey: Prentice Hall;1989.

[5] Honarpour M, Koederitz L, Harvey H. Relative permeability of petroleum reservoirs. Boca Raton Fla.: C.R.C. Press;1986.

[6] Hearn CL. Simulation of stratified waterflooding by pseudo relative permeability curves. Journal of Petroleum Technology 1971;23:805-813.

[7] Matthai SK, Mezentsev A, Belayneh M. Finite element-node-centred finite-volume two-phase-flow experiments with fractured rock represented by unstructured hybrid-element meshes. SPE Reservoir Evaluation & Engineering 2007;10: 740-756.

[8] Johnson EF, Bossler DP, Naumann VO. Calculation of relative permeability from displacement experiments. Trans. AIME 1959;216:370-372.

[9] Durlofsky JL. Numerical calculation of equivalent grid block permeability tensors for heterogenous porous media. Water Resources Research 1991;27:699-708.

[10] Jonoud S, Jackson MD. New criteria for the validity of steady-state upscaling. Transport in Porous Media 2008;71: 53-73

[11] Isherwood M, Harvey N. Integration of conventional petrophysical interpretation and borehole images. Back to Exploration CSPG CSEG CWLS Convention 2008.

[12] Huet CC, Newsham K, Rushing JA, Blasingame TA. A modified purcell/burdine model for estimating absolute permeability from mercury-injection capillary pressure data. International Petroleum Technology Conference 2005.

[13] Bazr-Afkan S, Matthai SK. A new hybrid simulation method for multiphase flow on unstructured grids with discrete representations of material interfaces. IAMG 2011.

[14] Ayachit U. The ParaView Guide: A Parallel Visualization Application. USA:Kitware;2015.