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

available literature data.

**3.1. Introduction** 

**3. Typical DRP Workflow** 

thin sections and micro-plugs.

**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

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

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

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

computed from Nano and micro scale to plug scale to a whole core scale.

input parameters were set according to expected wettability conditions.

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

available experimental data.

experimental data.

1.0

1.5

2.0

**Cementation Exp. (m)**

2.5

3.0

**Figure 4.** Porosity-Cementation exponent correlation

**Figure 3.** Porosity-FRF correlation

**Lab DRP**

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

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

0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40

**Porosity [Frac.]**

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

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


Network based flow relations


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

2-phase flow

206 New Technologies in the Oil and Gas Industry

Project workflow and calculated properties.

Grid based reservoir core properties

Absolute permeability kabs

wave velocity Vp and shear-wave velocity Vs)

Total porosity

**Figure 7.** Digital rock analyses on carbonate samples

analysed for each micro-plug, generating 3D digital rock models.

on network representations of the pore space (network based calculations):

Formation Resistivity Factor FRF and corresponding cementation exponent m

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

Micro-plug properties have been calculated directly on the digital rock (grid calculations) or

 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-


#### 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


Flow properties


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*

208 New Technologies in the Oil and Gas Industry

sub-plugs.

**3.4. Segmentation** 

software).

**3.3. Thin section preparation, sub-sampling** 

on micrite, pore sizes and cementation features.

*3.4.1. Step 1: Cropping and filtering* 

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

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

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

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

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

the corresponding micro-plug images. A noise reduction filter is applied.

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:

$$\phi\_{\text{vis}} = \sum\_{n=o}^{GV\_1 - 1} \mathbf{x}\_n \tag{1}$$

$$\phi\_{\mu} = \sum\_{n=GV\_1; m=1}^{GV\_2} \propto\_2 \left( 1 - \frac{m}{GV\_2 - GV\_1} \right) \tag{2}$$

$$
\phi\_{res} = \phi\_{vis} + \phi\_{\mu} \tag{3}
$$

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:

$$\sum\_{n=0}^{GV\_p} x\_n = f\_{\text{res}} \tag{4}$$

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

$$\phi\_{tot} = \phi \text{res} + \sum\_{n=GV\_p+1; m=GV\_2-GV\_p+1}^{GV\_2} \mathbf{x}\_n \left( 1 - \frac{m}{GV\_2 - GV\_1} \right) \tag{5}$$

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:


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*

#### *3.5.1.1. Porosity*

210 New Technologies in the Oil and Gas Industry

For this purpose, the grey value (GVp) is found where:

 

**Figure 9.** Practical pore threshold.

1.E+04

1.E+05

1.E+06

**Count**

1.E+07

1.E+08

3. representing a coarser image where

1. representing

2. representing

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.

*n res*

<sup>2</sup> 1; 1 2 1 1

(4)

res and matrix are represented as one phase and





**Count**

50000

100000

150000

200000

0

*m*

*GV GV*

(5)

*x f*

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.

The grey values GVp and GV2 give the segmentation thresholds for the micro-plug images.

0 50 100 150 200 250 **Grey Value**

**Porosity = 1 Porosity = 0**

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,

res and matrix (including remaining micro-porosity)

Three images are prepared according to this segmentation method:

res, micro-porosity and matrix

micro-porosity is segmented into 6 different porosity ranges

0

2

*GV*

*res x*

*tot n n GV m GV GV*

**GV0 GV1 GVp GV2**

*p p*

*GVp*

*n*

Three different porosities are reported: total porosity ( tot), micro-porosity ( <sup>μ</sup>) 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). <sup>μ</sup> 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 <sup>μ</sup> 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:

$$\begin{aligned} \phi\_{\text{tot,corplug}} &= \phi\_{\text{res,nugs}} + f\_{\mu-\text{plug},A}\phi\_{\text{tot,A}} + f\_{\mu-\text{plug},B}\phi\_{\text{tot,B}} + f\_{\mu-\text{plug},C}\phi\_{\text{tot,C}} \\ \phi\_{\text{tot,\mu-\text{plug}}} &= \phi\_{\mu} + \phi\_{\text{res}} \\ \phi\_{\text{res,\mu-\text{plug}}} &= \phi\_{\text{perc}} + \phi\_{\text{iso}} \end{aligned} \tag{6}$$

where *f* is fraction.

#### *3.5.1.2. 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).

#### *3.5.1.3. Formation resistivity factor*

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

$$\begin{aligned} \nabla \bullet \mathbf{J} &= 0 \\ \mathbf{J} &= \sigma\_w \nabla \Phi \end{aligned} \tag{7}$$

subject to the boundary condition n=0 on 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* = σ*i* /σ*w*. 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 reported 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 *<sup>µ</sup>* and *m* are the micrite porosity and the cementation exponent of the micrite phase.

*3.5.1.4. 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:

$$V\_p = \sqrt{\frac{K + \frac{4}{3}\mu}{\rho}}\tag{8}$$

$$V\_s = \sqrt{\frac{\mu}{\rho}}\tag{9}$$

Inputs to the calculations are:

212 New Technologies in the Oil and Gas Industry

difference method (*Øren and Bakke*, 2002).

**Pore-scale modeling of elastic properties** 

throughout the pore space (Arns *et al.*, 2002).

moduli and the effective density according to:

*-m*.

where

phase.

*Numerical code* 

*3.5.1.4. Elastic moduli* 

=0 = *<sup>w</sup>* 

subject to the boundary condition n=0 on 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

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* = σ*i* /σ*w*. 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* =

Formation factor and cementation exponent reported 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* (

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

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

*K*

*p*

*Vs*

*V*

4 3

(8)

(9)

*<sup>µ</sup>* and *m* are the micrite porosity and the cementation exponent of the micrite

 **J**

**<sup>J</sup>** (7)

σ*µ*)*m*,


#### **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.

#### *Stress*

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.

#### *Frequency*

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.

#### *3.5.1.5. 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*

214 New Technologies in the Oil and Gas Industry

collapse, grain rotation and grain crushing.

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

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

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

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

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 =

triangles. Laboratory measurements are shown as filled symbols; from Derzhi *et al*. (2011).

**Figure 11.** P-wave velocity and stress for carbonate samples.

rocks – and lower for fluid saturated rocks.

velocities are therefore overestimated for materials containing such cracks.

*Stress* 

*Frequency* 

*3.5.1.5. NMR* 

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.

#### *3.5.2.1. 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:

$$k\_{ro} = k\_{ro} \left( S\_{wi} \right) \frac{\left( S\_{on} \right)^{L\_o}}{\left( S\_{on} \right)^{Lo} + E\_o \left( 1 - S\_{on} \right)^{T\_o}} \tag{10}$$

$$k\_{rw} = k\_{rw} \left(\mathcal{S}\_{ov}\right) \frac{\left(\mathcal{S}\_{wu}\right)^{L\_w}}{\left(\mathcal{S}\_{wu}\right)^{L\_w} + E\_w \left(1 - \mathcal{S}\_{wu}\right)^{T\_o}}\tag{11}$$

$$\mathcal{S}\_{\rm uvn} = \frac{\mathcal{S}\_{\rm uv} - \mathcal{S}\_{\rm uvi}}{1 - \mathcal{S}\_{\rm uvi} - \mathcal{S}\_{\rm oy}}, \quad \mathcal{S}\_{\rm ov} = 1 - \mathcal{S}\_{\rm uvn} \tag{12}$$

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.

#### *3.5.2.2. 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:

$$J(S\_w) = \frac{P\_{c,ow}}{\sigma\_{ow}} \sqrt{\frac{k}{\nu}} \tag{13}$$

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):

$$J\_{\left(\mathcal{S}\_w\right)} = \frac{c\_1}{\left[\frac{\left(\mathcal{S}\_{av} - \mathcal{S}\_{wi}\right)}{\mathbf{1} - \mathcal{S}\_{wi}}\right]} + \frac{c\_2}{\left[\frac{\left(\mathbf{1} - \mathcal{S}\_{av} - \mathcal{S}\_{or}\right)}{\left(\mathbf{1} - \mathcal{S}\_{or}\right)}\right]}\tag{14}$$

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.

#### *3.5.2.3. 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:

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

$$\mathbf{S}\_{wi} = \frac{\phi\_{\text{isolated}} + \phi\_{\mu}}{\phi\_{\text{tot}}} + \mathbf{S}\_{w, corner} \tag{15}$$

$$S\_{w,corner} = \sum\_{1}^{i} \sum\_{1}^{n} \left(\frac{\sigma\_{ow}}{P c\_{ow, max}}\right)^{2} \left(\cos\ \theta\_{ow}^{\prime} \frac{\cos\left(\theta\_{ow}^{\prime} + \beta\right)}{\sin\beta} - \frac{\pi}{2} + \beta + \theta\_{ow}^{\prime}\right) \tag{16}$$

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.

#### *3.5.2.4. Resistivity Index and saturation exponent, n*

216 New Technologies in the Oil and Gas Industry

*3.5.2.2. Capillary pressure*

function:

where *k* and

entry pressures are much lower.

*wS*

Here, results are reported both as J-function and capillary pressure.

capillary pressure and

*3.5.2.3. Water saturation*

irreducible water saturation is given by:

*k kS*

parameters were determined using a simulated annealing algorithm.

*rw rw or L T*

, 1 <sup>1</sup> *w wi wn on wn wi oy S S S SS S S*

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

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,

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

 *c ow*, *w*

*c c <sup>J</sup> SS SS*

1 1

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.

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

1 2 1

*w wi w or wi or*

*S S* 

have been fitted to the empirical Skjæveland correlation (Skjæveland et al, 2000):

*ow <sup>P</sup> <sup>k</sup> J S* 

are the permeability and total porosity, respectively. *Pc,ow* is the oil-water

*ow* the oil-water interfacial tension. The Leverett *J*-function results

(13)

(14)

1

*S ES*

*S*

*wn w wn*

*w w o*

(12)

(11)

*L wn*

> 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 twophase 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

$$I\_{ij} = \frac{\mathcal{S}\_{ij}}{L\_{ij}} \left(\boldsymbol{\Phi}\_i - \boldsymbol{\Phi}\_j\right) \tag{17}$$

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

$$\frac{L\_{ij}}{\mathcal{S}\_{ij}} = \frac{l\_i}{\mathcal{S}\_i} + \frac{l\_t}{\mathcal{S}\_t} + \frac{l\_j}{\mathcal{S}\_j} \tag{18}$$

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

$$\alpha = \frac{2g\_t \left(\frac{1}{\mathcal{G}\_{ij}} - \frac{1}{\mathcal{G}\_i} - \frac{1}{\mathcal{G}\_j}\right)}{\mathcal{g}\_{ij} \left(1 - \frac{\mathcal{G}\_t}{\mathcal{G}\_i} - \frac{\mathcal{G}\_t}{\mathcal{G}\_j}\right)}\tag{19}$$

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

$$\sum\_{j} I\_{ij} = 0 \tag{20}$$

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* = *σ0/σw*, 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

$$RI\left(\mathcal{S}\_w\right) = \frac{\sigma\_0}{\sigma\left(\mathcal{S}\_w\right)} = \mathcal{S}\_w^{-n} \tag{21}$$

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

$$\begin{aligned} \nabla \bullet \begin{pmatrix} k \nabla P \end{pmatrix} &= O \quad \text{with boundary conditions} \\ p = P\_1 \text{ at } \mathbf{x} = o, p = P\_0 \text{ at } \mathbf{x} = L \\ \mathbf{v} \bullet \mathbf{n} &= \mathbf{0} \text{ at other faces} \end{aligned} \tag{22}$$

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:


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**

218 New Technologies in the Oil and Gas Industry

current conservation at each pore body, which means that

factor *F* = *σ0/σw*, where *σ0* is the conductivity computed at *Sw* = 1.

*3.5.3. Upscaling of single phase and effective flow properties* 

0

1 0 ,

*p P at x o p P at x L v n at other faces*

and *n*-exponent.

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

0 *ij*

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

The resistivity index is computed similarly. At various stages of the displacement (i.e.

0 *n w w w*

*RI S S S* 

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,

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

*k P O with boundary conditions*

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

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.

*<sup>I</sup>* (20)

(21)

(22)

*j*

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

