InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Engineering » Energy Engineering » "New Technologies in the Oil and Gas Industry", book edited by Jorge Salgado Gomes, ISBN 978-953-51-0825-2, Published: October 31, 2012 under CC BY 3.0 license. © The Author(s).

Chapter 9

Digital Rock Physics for Fast and Accurate Special Core Analysis in Carbonates

By Mohammed Zubair Kalam
DOI: 10.5772/52949

Article top


Simulated vs. experimental permeability
Figure 1. Simulated vs. experimental permeability
Simulated vs. experimental porosity
Figure 2. Simulated vs. experimental porosity
Porosity-FRF correlation
Figure 3. Porosity-FRF correlation
Porosity-Cementation exponent correlation
Figure 4. Porosity-Cementation exponent correlation
Porosity-Saturation exponent correlation
Figure 5. Porosity-Saturation exponent correlation
Porosity-Vs and Vp correlation
Figure 6. Porosity-Vs and Vp correlation
Digital rock analyses on carbonate samples
Figure 7. Digital rock analyses on carbonate samples
Example grey value histogram.
Figure 8. Example grey value histogram.
Practical pore threshold.
Figure 9. Practical pore threshold.
Fontainebleau sandstone.
Figure 10. Fontainebleau sandstone.
P-wave velocity and stress for carbonate samples.
Figure 11. P-wave velocity and stress for carbonate samples.
Water-oil Pc (Porous Plate): DRP vs lab on same core sample
Figure 12. Water-oil Pc (Porous Plate): DRP vs lab on same core sample
Water-oil Pc (Porous Plate): DRP vs lab in different core samples, but same RRT
Figure 13. Water-oil Pc (Porous Plate): DRP vs lab in different core samples, but same RRT
NMR T2 distribution and MICP pore throat distribution, DRP vs Lab – vuggy core
Figure 14. NMR T2 distribution and MICP pore throat distribution, DRP vs Lab – vuggy core
NMR T2 distribution and MICP pore throat distribution, DRP vs Lab – tight core
Figure 15. NMR T2 distribution and MICP pore throat distribution, DRP vs Lab – tight core
Validating water-oil kr of low permeability composite samples: RRT 6 (10-25 mD)
Figure 16. Validating water-oil kr of low permeability composite samples: RRT 6 (10-25 mD)
Validating water-oil kr of high permeability composite samples: RRT 8 (350-560 mD)
Figure 17. Validating water-oil kr of high permeability composite samples: RRT 8 (350-560 mD)

Digital Rock Physics for Fast and Accurate Special Core Analysis in Carbonates

Mohammed Zubair Kalam

1. Introduction

Initiatives for increasing hydrocarbon recovery from existing fields include the capability to quickly and accurately conduct reservoir simulations to evaluate different improved oil recovery scenarios. These numerical simulations require input parameters such as relative permeabilities, capillary pressures, and other rock and fluid porosity versus permeability trends. These parameters are typically derived from Special Core Analysis (SCAL) tests. Core analysis laboratories have traditionally provided SCAL through experiments conducted on core plugs. Depending on a number of variables, SCAL experiments can take a year or longer to complete and often are not carried out at reservoir conditions with live reservoir fluids. Digital Rock Physics (DRP) investigates and calculates the physical and fluid flow properties of porous rocks. In this approach, high-resolution images of the rock’s pores and mineral grains are obtained and processed, and the rock properties are evaluated by numerical simulation at the pore scale.

Comparisons between the rock properties obtained by DRP studies and those obtained by other means (laboratory SCAL tests, wireline logs, well tests, etc.) are important to validate this new technology and use the results it provides with confidence. This article shares a comparative study of DRP and laboratory SCAL evaluations of carbonate reservoir cores.

This technology is a breakthrough for oil and gas companies that need large volumes of accurate results faster than the current SCAL labs can normally deliver. The oil and gas companies can use this information as input to numerical reservoir simulators, fracture design programs, analytic analysis of PTA, etc. which will improve reserve forecasts, rate forecasts, well placement and completion designs. It can also help with evaluating option for improved oil recovery with sensitivity analysis of various options considering the actual pore scale rock fabric of each reservoir zone. Significant investment savings can also be realized using good DRP derived data compared with the conventional laboratory SCAL tests.

2. Background and summary

The main objective is to provide petrophysical and multiphase flow properties, calculated from 3D digital X-ray micro-tomographic images of the selected reservoir core samples. The simulations have been conducted on sub-samples (micro plugs) and then upscaled to cores plug scale for direct comparison with experimental data. Typically the whole core are imaged on dimensions of 11 – 16.5 cm with a resolution of 500 microns, while the core plugs are imaged from to 2 – 4 cm with resolutions of 12 – 19 microns. The micro plugs have dimensions from 1 – 5 mm with resolutions of 0.3 – 5 microns and with Nano-CT one can look at rocks of 50 – 300 microns with resolutions from 50 – 300 nm. In DRP process, the results of these increasingly smaller and smaller investigations are then integrated by an upscaling, either by steady state or geometric methods. Hence, rock properties are computed from Nano and micro scale to plug scale to a whole core scale.

Absolute permeability can be computed using Lattice-Boltzmann simulations, calculation of Formation Resistivity Factor is based on a solution of the Laplace equation with charge conservation (the equations were solved using a random walk algorithm) and elastic properties were calculated by the finite element method. Primary drainage and waterflood capillary pressure and relative permeabilities are determined from multi-phase flow simulations on the pore network representation of the 3D rock model. Flow simulation input parameters were set according to expected wettability conditions.

The first section outlines the basic DRP based results on reservoir properties determined on complex carbonates from giant Middle East reservoirs and compared with similar measurements performed in SCAL laboratories. Section 3 outlines possible details in the calculation process for each of the many reservoir parameters that can be calculated using DRP, while section 4 illustrates snapshots of multi-phase flow results on complex carbonates from the same giant Middle East reservoirs. Capillary pressure, cementation exponent (m), saturation exponent (n) for both primary drainage and imbibition, water-oil and gas-oil relative permeability and elastic properties of carbonates were calculated from DRP. Very good agreement is obtained between DRP derived properties and available experimental data for the studied data set. The results obtained for porosity, absolute permeability, formation resistivity factor, cementation and saturation exponent are shown in Figure 1. through 5. Calculations of elastic properties have been performed on all reconstructed samples. The elastic parameters Vs and Vp are reported in figure 6 and compared to available literature data.

3. Typical DRP Workflow

3.1. Introduction

In order to meet project objectives the workflow illustrated in Figure 7. was implemented. High resolution (19 μm/voxel) micro-CT images of core plugs were first recorded in order to identify rock types and their distribution within the core plugs, and to select locations for thin sections and micro-plugs.


Comparison of experimental and simulated permeability for all samples (straight line is 1:1 and dotted line is factor 2)

Figure 1.

Simulated vs. experimental permeability


Comparison of experimental and simulated porosity for 100+ samples (straight line is 1:1 and dotted line is +/- 3%)

Figure 2.

Simulated vs. experimental porosity


Calculated Formation Resistivity Factor (FRF) as a function of total porosity for the e-Core models compared to available experimental data.

Figure 3.

Porosity-FRF correlation


Calculated cementation exponent (m) as a function of total porosity for the e-Core models compared to available experimental data.

Figure 4.

Porosity-Cementation exponent correlation


Calculated saturation exponent (n) for primary drainage displacement as a function of total porosity for the e-Core models compared to available experimental data.

Figure 5.

Porosity-Saturation exponent correlation


Calculated Vs (open symbols) and Vp (colored symbol) as a function of total porosity compared to available literature data.

Figure 6.

Porosity-Vs and Vp correlation


Project workflow and calculated properties.

Figure 7.

Digital rock analyses on carbonate samples

Thin sections were then prepared and analysed for each core plug. Micro-plugs were drilled and high resolution (8.0 – 0.3 μm/voxel) micro-CT images were recorded, processed and analysed for each micro-plug, generating 3D digital rock models.

Micro-plug properties have been calculated directly on the digital rock (grid calculations) or on network representations of the pore space (network based calculations):

Grid based reservoir core properties

  • Total porosity ϕ

  • Absolute permeability kabs

  • Formation Resistivity Factor FRF and corresponding cementation exponent m

  • Elastic moduli assuming isotropy (bulk modulus K, shear modulus μ, Young’s modulus E, Poisson’s ratio υ and Lamé’s parameter λ) and corresponding acoustic velocities (P-wave velocity Vp and shear-wave velocity Vs)

Network based flow relations

  • Capillary pressure Pc as a function of water saturation Sw

  • Relative permeability kr as a function of water saturation Sw

  • Resistivity Index RI as a function of water saturation Sw and corresponding saturation exponent n for both primary drainage and imbibition cycle

Capillary pressure and relative permeability curves were established for the following flow processes:

2-phase flow

  • Oil/water primary drainage to initial water saturation Swi

  • Water/oil imbibition to residual oil saturation Sorw

3-phase flow

  • Gas/oil drainage at initial water saturation Swi

Resistivity index curves and corresponding saturation exponents were established for water/oil primary drainage and imbibition.

The last step of the workflow is to upscale micro-plug properties (volumes in the mm3 range) to core plug properties (volumes ranging from ~ 40 to100 cm3). Rock types were identified and micro-plugs representing these rock types were selected. Each core plug voxel was then assigned to a given rock type or pore space. Corresponding micro-plug and pore fluid properties were used as input to the calculation of upscaled properties.

The following properties were upscaled:

Petrophysical properties

  • Total porosity ϕ

  • Absolute permeability kabs

  • Formation Resistivity Factor FRF

  • Cementation Exponent m

Flow properties

  • Capillary pressure Pc as a function of water saturation Sw

  • Relative permeability kr as a function of water saturation Sw

  • Resistivity Index RI as a function of water saturation Sw

Capillary pressure, relative permeability and resistivity index curves were generated for the micro-plug flow processes listed on previous page.

3.2. Imaging

The principle of microfocus X-ray Computed Tomography (mCT) is based on Beer’s law, i.e. the intensity of X-rays is attenuated when passing through physical objects. The attenuated X-rays are captured by a detector to compose a projection image. A series of projection images from different angles (0 to 360°) are collected by rotating the object around its axis. These projection images are processed to generate 2D mCT slices, which subsequently are input date to construct 3D images of the objects. Each volume-element in these 3D images corresponds to one voxel. The voxel value of the recorded image is proportional to the attenuation coefficient, which is mainly a function of the density and the effective atomic number Z of the constituents of the object (Alvarez and Macovski, 1976; Pullan et al., 1981).

3.3. Thin section preparation, sub-sampling

The petrographic study of thin sections cut from the core plugs is an important part of the workflow. It allows describing and examining the heterogeneity of the core plugs on the micro- to macro-scale. Important features are the identification of micro-facies, the microscopic fabric in each facies type, and the mineralogy of the rock - including diagenetic features (cementation and secondary porosity). The first step in studying the facies associations is the analysis of the initial core plug mCT images (at 19 micron/voxel resolution) to identify the spatial distribution of the various facies, which is the basis for selecting where to cut the slab for thin section preparation, and also to select sites for drilling sub-plugs.

All samples are classified according to Dunham (1962) and Embry & Klovan (1971) limestone Classification. Several rock types are described with special attention to the grain types in terms of size, shape and sorting, because these parameters are influencing the flow properties. One of the main purposes of the microscopic analysis is to describe the porosity of the sample. The different pore types are identified using the Choquette and Pray (1970) porosity classification scheme. Furthermore, the high resolution images of the thin sections are used to distinguish between carbonate cements and micrite, in addition to identifying the main minerals (calcite or dolomite). Drilling sites for micro-plugs are also decided based on the characteristic microfacies types of the sample. The mCT scanning resolution of the micro plugs is chosen in accordance with the microscopic analysis, with particular emphasis on micrite, pore sizes and cementation features.

3.4. Segmentation

3.4.1. Step 1: Cropping and filtering

Scanned mCT images of plugs are cropped into the largest possible undisturbed rectangular volume avoiding cracks, uneven edges and gradients towards the outer surface of the sample. Cropped images are scaled to a voxel size representing approximately the size of the corresponding micro-plug images. A noise reduction filter is applied.

Scanned images of micro-plugs are cropped to a volume of 12003 grid cells to avoid gradients at the side of the image and to allow further processing (size limitations of applied software).

3.4.2. Step 2: Histogram Analysis

Each image element carries an 8-bit signal (256 grey values) corresponding to the X-ray attenuation experienced within its volume. A voxel completely filled by empty pore space (air) data is approaching black (0, ρair = 0.0012 g/cm3), while a voxel completely filled by high density mineral such as pyrite is approaching white (255, ρpyrite = 4.8 – 5.0 g/cm3). Voxels completely filled by minerals with lower densities, or voxels filled by various proportions of minerals and/or pore space cover the entire grey scale range from black to white. The grey value is consequently not determining the contents of the voxels uniquely; it may represent a single mineral, two or more minerals or porosity and minerals.

An example of a grey value histogram is shown in Figure 8. The histogram itself is shown as blue diamonds and is displayed on logarithmic scale (to the left). The other two data sets (pink and yellow) show the first and second order differentials of the grey value counts.


Figure 8.

Example grey value histogram.

The following porosities can then be calculated:


where x denotes the fraction of the entire image of the grey value and n and m are counting variables in a grey value range (n is grey value specific, m positive integer to distribute micro-porosity evenly).

This analysis underestimates porosity and permeability. Therefore, a practical threshold needs to be found to decide up until which micro-porosity value voxels are considered pore. For this purpose, the grey value (GVp) is found where:


This is illustrated in Figure 9. Thus, the total porosity of the image is calculated according to:


where the second term gives the amount of micro-porosity extracted from the image. Note that the micro-porosity values for individual grey values are the same in eq. 5 as in eq. 2.


Figure 9.

Practical pore threshold.

The grey values GVp and GV2 give the segmentation thresholds for the micro-plug images. Three images are prepared according to this segmentation method:

  1. representing ϕres and matrix (including remaining micro-porosity)

  2. representing ϕres, micro-porosity and matrix

  3. representing a coarser image where ϕres and matrix are represented as one phase and micro-porosity is segmented into 6 different porosity ranges

Image 1 is segmented into a pore network representation (see Bakke and Øren, 1997, Øren and Bakke, 2002, 2003) which is used to calculate capillary pressure, relative permeability, resistivity index, and saturation exponent. Image 2 is used to calculate porosity, permeability, formation resistivity factor, and cementation exponent. Image 3 is used in a steady-state upscaling routine to include properties calculated for generic models for micritic material representing the micro-porosity.

The scaled image of the whole core plug is segmented in the same way to identify фres as vugs. However, the remaining matrix of the core plug is segmented into pure solid and 1-3 porosity classes according to the segmentation of the respective micro-plugs extracted from the core plug.

3.5. Calculations/Simulations

3.5.1. Single-phase properties Porosity

Three different porosities are reported: total porosity (ϕtot), micro-porosity (ϕµ) and percolating porosity (ϕperc). ϕperc is the fraction of ϕres (see above) that is available to flow, i.e. accessible from any side of the micro-plug image, expressed as a fraction of the whole micro-plug volume. The fraction of ϕres that is not available for flow is isolated porosity (ϕiso). ϕµ is all porosity present in voxels between GVp and GV2 expressed as a fraction of the whole micro-plug volume. ϕperc is only reported for the micro-plug images. ϕtot and ϕµ are reported for both the micro-plugs and the whole core plug. Some micro-plug images with the highest resolution can be considered entirely as micro-porosity. The following relations hold for porosities of micro-plug and core plug images:


where f is fraction. Absolute permeability

A Lattice-Boltzmann method is applied to solve Stokes’ equation in the uniform grid model. Flow is driven either by a constant body force or a constant pressure gradient through the model. Permeability is calculated in three orthogonal directions separately. Sides perpendicular to flow directions are closed during each directional calculation (no-flow boundary conditions). In this study, absolute permeability is calculated using a constant body force because this setting delivers more accurate results when model resolution is sufficient and porosity is relatively high. Averages of three directional calculations are reported for each model realization. For further details see Øren and Bakke (2002). Formation resistivity factor

The steady state electrical conductivity, or formation resistivity factor (F), of a brine saturated rock is governed by the Laplace equation


subject to the boundary condition Φn=0on the solid walls (i.e. insulating walls). Here J is the electrical current, σw is the electrical conductivity of the fluid that fills the pore space, Ф is the potential or voltage and n is the unit vector normal to the solid wall. Numerical solutions of the Laplace equation are obtained by a random walk algorithm or by a finite difference method (Øren and Bakke, 2002).

The effective directional conductivities σi, i = x, y, z are computed by applying a potential gradient across the sample in i-direction. The directional formation resistivity factor Fi is the inverse of the effective electrical conductivity Fi = σiw. We define the average formation resistivity factor F as the harmonic mean of direction dependent formation factors. The cementation exponent m is calculated from F and the sample porosity using Archie’s law F = ϕ-m.

Formation factor and cementation exponent were approximately 10-15% greater than experimental data. The F and m values reported in the previous version were calculated as described above, i.e. by assuming that only the resolved porosity contributes to the conductivity. Micro porosity present in the micrite phase was thus treated as insulating solid. In the second version, we accounted for the conductivity of the micrite phase by assigning a finite conductivity σmic to micrite voxels using Archie’s law σµ = σw (ϕσµ)m, where µ and m are the micrite porosity and the cementation exponent of the micrite phase. Elastic moduli

Pore-scale modeling of elastic properties

Numerical code

The finite element method described by Garboczi and Day (1995) has been implemented for calculation of elastic properties. The method uses a variational formulation of the linear elastic equations and finds the solution by minimizing the elastic energy using a fast conjugate-gradient method. The results are valid for quasi-static conditions or at frequencies which are sufficiently low such that the included pore pressures are in equilibrium throughout the pore space (Arns et al., 2002).

The effective bulk and shear moduli are computed assuming isotropic linear elastic behavior. The Vp and Vs are subsequently calculated using the simulated effective elastic moduli and the effective density according to:

Inputs to the calculations are:

  • A three-dimensional representation of the rock microstructure (a digital rock sample)

  • Density ρ and elastic properties (K and μ) of each mineral composing the rock matrix

  • Density ρ and elastic properties (K and μ) of the fluid present in the pore space

Pore-scale versus experimental data

Scale and sample selection

Acoustic measurements are generally performed on samples with volumes ranging from 40 to 100 cm3. Digital rock samples are significantly smaller, generally in the mm3 range. Care must therefore be taken when laboratory measurements and pore-scale derived properties are compared, due to the scale difference. Laboratory samples and plugs for μCT imaging are not only of different volumes, they are also generally sampled at different locations. Due to the spatial variability, it is recommended to compare trends when laboratory measurements are compared with pore-scale derived properties. This is illustrated in Figure 10 for the P-wave velocity. A digital sample with 12.9% porosity has been made. However, samples tested in the laboratory do not have porosities close to this value. By sub-sampling the original digital sample, a relative broad porosity range is obtained (from 9.9 to 16.3%). The pore-scale derived velocity – porosity trend is now overlapping with the measurements performed in the laboratory, and a comparison can be made.


Pore-scale derived P- and S-wave velocities (dark and light blue symbols) are compared with laboratory measurements (green filled symbols). Dark blue points represent a large digital “mother” sample, while the light blue points represent sub-samples of the original “mother” sample.

Figure 10.

Fontainebleau sandstone.


Cracks may reduce the acoustic velocities significantly. Sub-resolution cracks are not incorporated in the processed mCT images and the corresponding pore-scale derived velocities are therefore overestimated for materials containing such cracks.

Cracks are closed during loading. Acoustic measurements performed at elevated stress levels are consequently expected to approach pore-scale derived velocities. Derzhi and Kalam (2011) compared acoustic measurements at different stress levels with pore-scale derived velocities. Their results are shown in Figure 11. Note that this assumes that the pore space and rock framework deforms without large micro-structural changes such as pore collapse, grain rotation and grain crushing.


P-wave velocity as a function of porosity at different stress levels: a) at ambient conditions, b) at 7 MPa, c) at 30 MPa and d) at 40 MPa effective hydrostatic stress. Pore-scale derived values are denoted DRP and shown as open orange triangles. Laboratory measurements are shown as filled symbols; from Derzhi et al. (2011).

Figure 11.

P-wave velocity and stress for carbonate samples.


Pore-scale derived elastic properties represent properties in the low frequency limit (f → 0). Measurements in the laboratory are generally performed in the ultrasonic range (1 Hz to 4 kHz). Acoustic velocities are independent of frequency for dry materials, while an increase with increasing frequency has been observed for fluid saturated rocks. Pore-scale derived properties are therefore expected to be comparable to laboratory measurements for dry rocks – and lower for fluid saturated rocks. NMR

NMR is simulated as a diffusion process using a random walk algorithm to solve the diffusion equations (Øren, Antonsen, Rueslåtten and Bakke, 2002). The T2 responses at Sw = 1.0 were simulated on one sample from the field A. The results are shown as T2 decay and T2 distribution curves in figures 14 and 15. The surface relaxation strength was kept at 1.65x10-5 m/sec for all minerals. The inter echo time was 200µsec at a background magnetic field gradient of 0.2 G/m, the bulk water T2 was 0.3 sec, and the diffusion constant was 2x10-9 m2/sec. The decay curves are transformed into T2 distributions using 50 exponential functions (assuming the same has been done for the lab data).

3.5.2. Multi-phase flow simulations

The reconstructed rock models were simplified into pore network models. Crucial geometrical and topological properties were retained, while the data volume was reduced to allow timely computation (Øren and Bakke, 2003). In pore network modelling, local capillary equilibrium and the Young–Laplace equation are used to determine multiphase fluid configurations for any pressure difference between phases for pores of different shape and with different fluid/solid contact angles. The pressure in one of the phases is allowed to increase and a succession of equilibrium fluid configurations are computed in the network. Then, empirical expressions for the hydraulic conductance of each phase in each pore and throat are used to define the flow of each phase in terms of pressure differences between pores. Conservation of mass is invoked to find the pressure throughout the network, assuming that all the fluid interfaces are fixed in place. From this the relationship between flow rate and pressure gradient can be found and hence macroscopic properties, such as absolute and relative permeabilities, can be determined (Øren and Bakke, 2002).

The following oil-water displacements were simulated:

Primary drainage

Imbibition at a given wettability preference

For each displacement process, capillary pressure and relative permeability curves were calculated. Resistivity index with the corresponding n-exponent was calculated after primary drainage.

Each saturation and relative permeability value corresponds to a capillary equilibrium state. In all the simulations, it is assumed that capillary forces dominate. This is a good approximation for capillary numbers Nca < 1e-6. This, however, does not necessarily mean that it is the most efficient displacement possible. In certain cases, viscous and/or gravity forces can dominate and result in higher (or lower) displacement or sweep efficiency. Relative permeability

Simulated relative permeability data are fitted to the empirical LET expression (Lomeland, Ebeltoft and Thomas, 2005). For water/oil displacements the LET equations become:


where kro and krw are the relative permeability of oil and water, respectively. The Li’s, Ei’s, and Ti’s are the LET fitting parameters, where i is either oil (o) or water (w). kro(Swi), krw(Sor), Swi and Sor are determined from the computed results and the optimised values of the fitting parameters were determined using a simulated annealing algorithm. Capillary pressure

A detailed account of the methods used to calculate capillary pressure as a function of Sw during primary drainage and waterflooding invasion sequences is given in Øren et al. (1998). Fluid injection is simulated from one side of the model (usually x-direction). Thus, the entry pressure is a function of the pore sizes present in the inlet. In a mercury injection capillary pressure simulation, fluid is allowed to enter the model from all sides. In that case, entry pressures are much lower.

The calculated capillary pressures can be expressed in terms of the dimensionless Leverett J-function:


where k and ϕ are the permeability and total porosity, respectively. Pc,ow is the oil-water capillary pressure and σow the oil-water interfacial tension. The Leverett J-function results have been fitted to the empirical Skjæveland correlation (Skjæveland et al, 2000):


where c1, c2, a1 and a2 are curve fitting parameters. Sw is the water saturation, Swi initial water saturation and Sor residual oil saturation. Sw, Swi and Sor is determined by the simulations. Here, results are reported both as J-function and capillary pressure. Water saturation

All reported Sw is total Sw, i.e. including water in microporosity and isolated pores. It should be noted that Swi is strongly dependent on capillary pressure. Thus, any comparison with laboratory data should be done at the same capillary pressure. Any isolated pore volume and any µ-porosity cannot be invaded and, thus, contributes to Swi. Therefore, the irreducible water saturation is given by:


where Φ is porosity (with suffixes denoting isolated, total and microporosity), i is the number of connected pores in the network, n the number of corners per pore (3 or 4), σ interfacial tension, Pc capillary pressure, θ contact angle and β the corner half angle. Note that initial water saturation for waterflooding depends on Pc and can be given as an input to the simulation if needed. Resistivity Index and saturation exponent, n

The resistivity index is calculated from capillary dominated two-phase flow simulations on the pore network representation of the 3D rock image. The basis for simulating capillary dominated displacements is the correct distribution of the fluids in the pore space. For two-phase flow, the equilibrium fluid distribution is governed by wettability and capillary pressure and can be found by applying the Young-Laplace equation for any imposed pressure difference between the phases. A clear and comprehensive discussion of all the mathematical details, including the effects of wettability, involved in the simulations can be found in Øren et al., 1998, and Øren and Bakke, 2003).

The current I between two connecting nodes i and j in the network is given by Ohm’s law


where Lij is the spacing between the node centres. The effective conductance gij is the harmonic mean of the conductances of the throat and the two connecting nodes


where the subscript t denotes the pore throat and the conductance gt is evaluated at the throat constriction. The effective lengths li, lj, and lt govern the potential drop associated with the nodes and the throat, respectively. By letting lt = αLij and li = lj = 0.5(1-α)Lij, the effective lengths can be calculated from α as


The conductance of a pore element k (pore body or throat) is given by gk = σwAw, where Aw is the area of the pore element filled with water. Expressions for Aw for different fluid configurations, contact angles, and pore shapes are given in Øren et al., 1998. We impose current conservation at each pore body, which means that

where j runs over all the pore throats connected to node i. This gives rise to a set of linear equations for the pore body potentials. The formation resistivity factor of the network is computed by imposing a constant potential gradient across the network and let the system relax using a conjugate gradient method to determine the node potentials. From the potential distribution one may calculate the total current and thus the formation resistivity factor F = σ0w, where σ0 is the conductivity computed at Sw = 1.

The resistivity index is computed similarly. At various stages of the displacement (i.e. different Sw values), we compute the current and the resistivity index defined as


The n-exponent is determined from a linear regression of the RI(Sw) vs. Sw. curve.

3.5.3. Upscaling of single phase and effective flow properties

Effective properties of the core samples are determined using steady state scale up methods. The CT scan of the core plug is gridded according to the observed geometrical distribution of the different rock types or porosity contributors. Each grid cell is then populated with properties calculated on the pore scale images of the individual rock types. The following properties are assigned to each grid cell; porosity, absolute permeability tensor (kxx, kyy, kzz), directionally dependent m-exponents, capillary pressure curve, relative permeability curve, and n-exponent.

Single phase up-scaling is done by assuming steady state linear flow across the model. The single phase pressure equations are set up assuming material balance and Darcy’s law


The pressure equation is solved using a finite difference formulation. From the solution one can calculate the average velocity and the effective permeability using Darcy’s law. By performing the calculations in the three orthogonal directions, we can compute the effective or up-scaled permeability tensor for the core sample. The effective formation resistivity factor is computed in a similar manner by replacing pressure with voltage, flow with current, and permeability with electrical conductivity. The up-scaled m-exponent is determined from the effective formation resistivity factor and the sample porosity.

Effective two-phase properties (i.e. capillary pressure, relative permeability, and n-exponent) are calculated using two-phase steady state up-scaling methods. We assume that the fluids inside the sample have come to capillary equilibrium. This is a reasonable assumptions for small samples (<30cm) when the flow rate is slow (<1m/day). The main steps in the two-phase up-scaling algorithm are:

  1. Select a capillary pressure (Pc) level

  2. Using the Pc (Sw) relationship, calculate Sw in each grid cell

  3. Calculate the average water saturation using pore volume weighting

  4. From the kr(Sw) curves, calculate krw and kro in each grid cell, and then the phase mobilities kw and ko by multiplying the relative permeabilities with the absolute permeability for the grid cell

  5. Perform two separate single phase steady state simulations, one for the oil and one for the water, to calculate the effective phase permeability

  6. Divide the phase permeability with the effective absolute permeability to obtain the effective relative permeability

  7. Repeat these steps with different Pc levels to construct the effective relative permeability curves

The resistivity index curve is generated in a similar manner by replacing phase permeability in the water phase calculations with electrical conductivity. The effective n-exponent is determined from a linear regression of the effective RI(Sw) vs. Sw. curve.

3.6. Uncertainty

The overall uncertainty in up-scaled properties varies from sample to sample. Uncertainties may be introduced in the following steps:

  • Generation of digitized core models

How representative are identified rock types and their corresponding distribution? In other words; how representative are the digitized model of the core samples?

  • Generation of digitized micro-plugs

How representative are μCT models of the rock micro-structure? The main uncertainty is related to size (REV) and image segmentation where the spatial distribution of pore-space and rock minerals is set.

  • Calculation of absolute permeability and formation resistivity factor

An uncertainty of ± 2% related to the accuracy of pressure solvers.

  • Calculation of elastic properties

Uncertainty related to whether all relevant physics are included in the calculations or not.

  • Simulation of two- and three-phase flow

Wettability is an input parameter to the simulations. The uncertainty in wettability is the main source of uncertainty for both two- and three-phase flow simulations.

4. DRP applications in multi-phase flow

In this section, we present some novel results of validation of multi-phase flow SCAL results using Digital Rock Physics. The reservoir cores comprise complex carbonates from giant producing reservoirs in Middle East. Figure 12 show the comparisons of water-oil capillary pressure (Pc) measured in a SCAL laboratory at reservoir temperature and net over burden pressure using a Porous Plate and MICP trims from the same cores corrected to the reservoir conditions. The DRP data were acquired from the cores after the tests were completed on the Porous Plate and core thoroughly cleaned for final SCAL reference measurements. Both limestone and dolomite samples show excellent similarity of DRP derived data with the laboratory evaluations. Figure 13 confirms the validity of such measurements on different sets of core samples comprising the same reservoir rock type (RT), provided the rock typing is valid and captures the key formation properties of rock and fluids.


Figure 12.

Water-oil Pc (Porous Plate): DRP vs lab on same core sample

Figures 14 and 15 show for the first time in industry that laboratory NMR T2 and MICP measurements done on carbonate rock types can also be captured using DRP based simulations on the same cores with distinctly different pore geometries. The robustness of DRP in capturing NMR T2 based pore bodies and MICP based pore throat distributions have far reaching consequences. This shows that DRP in essence can be used confidently to quantify pore body and pore throat distributions, and therefore the 3D pore geometry is representative of the specific core sample and pore network topology. In using DRP effectively, it is recommended that one compares and validates measured NMR T2 and MICP prior to detailed simulations to quantify various two-phase and three-phase flow properties through such reservoir rocks.

Figures 16 and 17 demonstrate example DRP based validations with respect to water-oil relative permeabilities conducted at full reservoir conditions (reservoir temperature, reservoir pressure and live fluids) on other complex carbonates, including highly permeable vuggy samples. The imbibition displacements were conducted under steady state conditions at SCAL laboratories and QC’ed thoroughly with respect to production, pressure profiles and insitu saturation data, and the corresponding numerically simulated measrements. The DRP data were acquired on cores comprising each of the composites tested. It is interesting to note that when plug DRP data are compared with composite laboratory measurements there is some scatter and divergence for each reservoir rock type. However, the divergences are significantly minimized when the DRP plugs used are digitally butted to represent the composite used in the laboratory tests. DRP captures the full reservoir condition multi-phase flow data very well, and in some cases even show the experimental artefacts of the SCAL measurements. The validity of tehse tests were confirmed on 14 different reservoir rock types comprising different formations.


Figure 13.

Water-oil Pc (Porous Plate): DRP vs lab in different core samples, but same RRT


Figure 14.

NMR T2 distribution and MICP pore throat distribution, DRP vs Lab – vuggy core


Figure 15.

NMR T2 distribution and MICP pore throat distribution, DRP vs Lab – tight core


Figure 16.

Validating water-oil kr of low permeability composite samples: RRT 6 (10-25 mD)


Figure 17.

Validating water-oil kr of high permeability composite samples: RRT 8 (350-560 mD)


ADCO and ADNOC Management are acknowledged for their permission to publish these novel Digital Rock Physics based SCAL results.

Numerical Rocks (Norway) is acknowledged for providing the detailed drafts relating to the procedures adopted in example DRP computations and the robust measurements presented in this chapter.

Ingrain Inc (Houston and Abu Dhabi) are acknowledged for introducing the author to the various challenges ahead, and the uncertainties in current DRP based dvelopments.

Digital Core (Australia) is gratefully remembered in first introducing the concept of DRP to Middle East, and involving ADCO in one of the first Joint Industy Projects offered to industry.

iRocks (Beijing and London) are acknowledged for many stimulating discussions relating to the state-of-the-art.

Finally, one must remember the ADCO DRP team for the interest generated and the numerous insights to imaging, segmentation, data evolution and their impact on different Petrophysical, SCAL and elastic properties. I thank my son, AbdAllah, for helping me in getting this draft ready despite the very busy schedules of August 2012.


1 - Bakke, S., Øren, P.-E., “3 -D Pore-Scale Modelling of Sandstones and Flow Simulations in the Pore Networks,” SPE Journal1997
2 - Blunt M., Jackson M.D., Piri M., Valvatne P.H., “Detailed physics, predictive capabilities and macroscopic consequences for pore network models of multiphase flow,” Advances in Water Resources2002
3 - Ghous, A., Knackstedt, M.A., Arns, C.H., Sheppard, A.P., Kumar, R.M., Senden, T.J., Latham, S., Jones, A.C., Averdunk, H. and Pinczewski, W.V., “3-D imaging of reservoir core at multiple scales: Correlations to petrophysical properties and pore-scale fluid distributions,” presented at International Petroleum Technology Conference, Kuala Lumpur, Malaysia, 2008, 10 p.
4 - Gomari, K. A. R., Berg, C. F., Mock, A., Øren, P.-E., Petersen, E. B. Jr., Rustad, A. B., Lopez, O., 2011, Electrical and petrophysical properties of siliciclastic reservoir rocks from pore-scale modeling, paper SCA2011-20 presented at the 2011 SCA International Symposium, Austin, Texas.
5 - A. Grader, M. Z. Kalam, J. Toelke, Y. Mu, N. Derzhi, C. Baldwin, M. Armbruster, Dayyani. T. Al, A. Clark, Yafei. G. B. Al, Stenger. B. And, 2010A comparative study of DRP and laboratory SCAL evaluations of carbonate cores, paper SCA201024presented at the 2010 SCA International Symposium, Halifax, Canada.
6 - M. Z. Kalam, Dayyani. T. Al, A. Grader, C. Sisk, 2011Digital rock physics analysis in complex carbonates’, World Oil, May 2011.
7 - Kalam M.Z., Serag S., Bhatti Z., Mock A., Oren P.E., Ravlo V. and Lopez O., SCA2012-03, “Relative Permeability Assessment in a Giant Carbonate Reservoir Using Digital Rock Physics,” SCA 2012 International Symposium, Aberdeen, United Kingdom.
8 - Kalam, M.Z., Al-Hammadi, K., Wilson, O.B., Dernaika, M., and Samosir, H., “Importance of Porous Plate Measurements on Carbonates at Pseudo Reservoir Conditions,” SCA2006-28, presented at the 2006 SCA International Symposium, Trondheim, Norway.
9 - Kalam M.Z., El Mahdi A., Negahban S., Bahamaish J.N.B., Wilson O.B., and Spearing M.C., “A Case Study to Demonstrate the Use of SCAL Data in Field Development Planning of a Middle East Carbonate Reservoir,”SCA200618presented at the 2006SCA International Symposium, Trondheim, Norway.
10 - Kalam, M. Z., Al Dayyani, T., Clark, A., Roth, S., Nardi, C., Lopez, O. and Øren, P. E., “Case study in validating capillary pressure, relative permeability and resistivity index of carbonates from X-Ray micro-tomography images”, SCA2010-02 presented at the 2010 SCA International Symposium, Halifax, Canada.
11 - Knackstedt, M.A., Arns, C.H., Limaye, A., Sakellariou, A., Senden, T.J., Sheppard, A.P., Sok, R.M., Pinczewski, W.V. and Bunn G.F., “Digital core laboratory: Reservoir-core properties derived from 3D images," Journal of Petroleum Technology, (2004) 56, 66-68.
12 - Lopez, O., Mock, A., Øren, P. E., Long, H., Kalam, M. Z., Vahrenkemp, V., Gibrata, M., Serag, S., Chacko, S., Al Hosni, H., Al Hammadi M. I., and Vizamora, A., “Validation of fundamental carbonate reservoir core properties using Digital Rock Physics”, SCA2012-19, SCA 2012 International Symposium Aberdeen, United Kingdom.
13 - O. Lopez, A. Mock, J. Skretting, E. B. Petersen, P. E. Jr Øren, A. B. Rustad, Investigation into the reliability of predictive pore-scale modeling for siliciclastic reservoir rocks, SCA201041presented at the 2010SCA International Symposium, Halifax, Canada.
14 - Mu, Y., Fang, O., Toelke. J., Grader, A., Dernaika, M., Kalam, M.Z., ‘Drainage and imbibition capillary pressure curves of carbonate reservoir rocks by digital rock physics’, SCA 2012- Paper A069, Aberdeen, United Kingdom.
15 - Øren, P.E. and Bakke, S., Process Based Reconstruction of Sandstones and Prediction of Transport Properties, Transport in Porous Media,2002200246311343
16 - P. E. Øren, F. Antonsen, H. G. Rueslåtten, S. Bakke, 2002Numerical simulations of NMR responses for improved interpretation of NMR measurements in rocks, SPE paper 77398, presented at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas.
17 - P. E. Øren, S. Bakke, O. J. Arntzen, 1998Extending predictive capabilities to network models, SPE J., 3324336
18 - Øren, P.E. and Bakke, S., “Process Based Reconstruction of Sandstones and Prediction of Transport Properties,” Transport in Porous Media,2006
19 - Øren, P.E., Antonsen, F., Rueslåtten, H.G., and Bakke, S., “Numerical simulations of NMR responses for improved interpretation of NMR measurements in rocks,” SPE paper 77398, presented at the2002SPE Annual Technical Conference and Exhibition, San Antonio, Texas.
20 - Øren, P. E., Bakke, S. and Arntzen, O. J., “Extending predictive capabilities to network models,” SPE Journal, (1998) 3, 324-336.
21 - T. Ramstad, P. E. Øren, S. Bakke, 2010Simulations of two phase flow in reservoir rocks using a Lattice Boltzmann method, SPE J., SPE 124617.
22 - Youssef, S., Bauer, D., Bekri, S., Rosenberg, E. and Vizika, O., “Towards a better understanding of multiphase flow in porous media: 3D in-situ fluid distribution imaging at the pore scale,” SCA2009-17, presented at the 2009 SCA International Symposium, Noordwijk, The Netherlands.
23 - K. Wu, Z. Jiang, G. D. Couples, M. I. J. Van Dijke, K. S. Sorbie, 2007Reconstruction of multi-scale heterogeneous porous media and their flow prediction,” SCA200716presented at the 2007 SCA International Symposium, Calgary, Canada.
24 - Wu, K., Ryazanov, A., van Dijke, M.I.J., Jiang, Z., Ma, J., Couples, G.D., and Sorbie, K.S., “Validation of methods for multi-scale pore space reconstruction and their use in prediction of flow properties of carbonate,” SCA-2008-34 presented at the 2008 SCA International Symposium in Abu Dhabi, UAE.