## Abstract

A statistical technique for the pore-scale analyses of heterogeneity and representative elemental volume (REV) in unconventional shale rocks is hereby presented. First, core samples were obtained from shale formations. The images were scanned using microcomputed tomography (micro-CT) machine at 6.7 μm resolution with voxels of 990 × 990 × 1000. These were then processed, digitised, thresholded, segmented and features captured using numerical algorithms. This allows the segmentation of each sample into four distinct morphological entities consisting of pores, organic matter, shale grains and minerals.

#### Author

##### James O. Adeleye, Lateef T. Akanji

Petroleum Engineering Division, School of Engineering, University of Aberdeen, Aberdeen, UK.

Received: 2 February 2017 / Accepted: 23 July 2017

© The Author(s) 2017.

In order to analyse the degree of heterogeneity, Eagle Ford parallel sample was further cropped into 96 subsamples. Descriptive statistical approach was then used to evaluate the existence of heterogeneity within the subsamples. Furthermore, the Eagle Ford parallel and perpendicular samples were analysed for volumetric entities representative of the petrophysical variable, porosity, using corner point cropping technique. The results of porosity REV for Eagle ford parallel and perpendicular indicated sample representation at 300 μm voxel edge. Both pore volume distribution and descriptive statistical analyses suggested that a wide variation of heterogeneity exists at this scale of investigation. Furthermore, this experiment allows for adequate extraction of necessary information and structural parameters for pore-scale modelling and simulation. Additional studies focusing on re-evaluation at higher resolution are recommended.

## Introduction

Pore geometry, tortuosity, grains size and shape are properties that are important to describe and characterise fluid flow in shale rock. However, getting a single material point measurement at which this shale rock property is determined depends on the ability to accurately extract information from its structure.

Determination of the nature and extent of heterogeneities at pore scale can enhance fluid characterisation in porous media. In shale rock, variability is extreme, because shales are detrital sediments formed as a result of alteration of mud or clay deposits. They occur as fine-grained rich in illite and fragments with particle sizes generally less than 0.062 mm (Leith 2016). Hence, its petrophysical properties (lithology, porosity, permeability water saturation) become difficult to evaluate.

Another major challenge of shale petrophysical evaluation is its heterogeneity as this makes quantification of porosity and other properties difficult. These parameters depend on knowing the details of how shape or size grains (Milner et al. 2010), pore throats (Curtis et al. 2012) and tortuosity (Bai et al. 2013; Katsube et al. 1991) are distributed. Several studies have reported the complex microstructure of the shale rocks (Curtis et al. 2012), heterogeneity of shale rocks (Chen et al. 2013) as well as microcracks (Bai et al. 2013) without necessarily quantifying the degree of heterogeneity. Extensive research into the effects of particle shape and polydispersity on flow through porous media using synthetic models has shed more light on how pore geometry and complexity of different porous media affect porosity and permeability (Mota et al. 2005; Garcia et al. 2009). However, little has been performed to fully characterise same in unconventional shale geometries.

Microcomputed tomography (micro-CT) is one of the most powerful non-destructive imaging techniques that can be used to characterise the microstructure of porous materials at the microscopic scale (Rozenbaum and du Roscoat 2014). Over the years, this technique has contributed immensely to the understanding and characterisation of fluid flow through pore geometries of shale rocks (Watson and Mudra 1994; Boruah and Ganapathi 2015). Ma et al. (2016) recently reported that they were able to scan sample size of 100 µm and obtained spatial pixel size of 130 nm. This, however, brings about the question of the ability of the acquired images to adequately represent the microstructural characteristics of the data sets. In order to simulate fluid flow processes, an average of seven grains diameter is required for a fully developed representative elemental volume (REV) (Garcia et al. 2009).

Replacing a heterogeneous property of the rock with an equivalent homogeneous one through a continuum description could best be done through REV method, which represents appreciable property to capture the heterogeneity (Bear 1972). This analysis is carried out by plotting different sample volumes with their corresponding measured property. As shown in Fig. 1, the measured property (n) varies intensely with small changes in the sample length (L) and begins to reduce until measured property is relatively insensitive to small changes in volume or location.

**Fig. 1** Graphical representation of how representative elementary volume (REV) is determined for a specific property (Bear 1972).

Then, representative amount of porosity can be confidently determined (Bear 1972). According to Sahimi (1995), since the continuum approach breaks down if the correlation lengths in the system approach the size of the system, it is important to evaluate the effect of various length scales at which the system may be considered homogeneous.

In this paper, we developed a workflow for the characterisation of heterogeneity of unconventional shale rock samples. The degree of heterogeneity was evaluated using 1 microcomputed tomography (micro-CT) image from Eagle Ford parallel. The sample was divided into 96 for descriptive statistical analysis. Meanwhile, REV for porosity was evaluated using a set of two microcomputed tomography (micro-CT) images from Eagle Ford parallel and perpendicular. Finally, corner point reference techniques where each of the sample was cropped into 7 were used to achieve REV for porosity. The images were processed, digitised, thresholded, segmented and each feature captured using marching cube algorithm. This allows the segmentation of each sample into four distinct entities, consisting of pores, organic matter, shale grains and minerals.

## Image acquisition and processing

The 3D image of the samples used for the analysis was obtained using an industrial micro-CT device phoenix v|tome|x s. It has a 180-kV, micro- and nano-focus X-ray tube and digital detector array (1000 × 1000 pixels). The flow chart for this study is shown in Fig. 2. Full sample (Fig. 2a) was placed in the sample holder, and 6.7 mm size was scanned by rotating it 360°, and the signal-to-noise ratio used for the acquisition of these data was 37.8%. The resulting projections were converted into a 3D image stack using the PhoenixImaging resolution 3D reconstruction software developed from filtered backprojection Feldkamp algorithm (Feldkamp et al. 1984). The detailed procedure for the sample imaging can be found in Singhal et al. (2013).

**Fig. 2** a Eagle Ford (EGF) shale rock sample obtained parallel to the bedding plane. b Micro-CT image of the shale sample at resolution of 6.7 μm showing 3D volume raw data with artefacts. c The 3D volume after artefacts have been removed. d Graphical annotation design of six (6) REV image windows of 50, 100, 200, 300, 400, 500 and 600 μm pixel length. e Graphical annotation design of ninety-six (96) 1 mm × 1 mm × 1 mm subsamples for REV and heterogeneity analyses, respectively. f Sample description used for heterogeneity analysis. Note: The length of each subsample is denoted by lx where x∈{1,2,3,4,⋯,n} and total length L=∑^{n}_{x=1}lx.

In order to accurately estimate petrophysical properties of the samples, proper segmentation of solid(s) and void phases takes a priority. Segmentation involves the segregation of the grey-level voxels of the 3D image into distinct phases. The presence of artefacts such as streaks, brightness, non-uniformity, and phase-contrast fringes at edges and/or noise would reduce the accuracy of segmentation of these 3D images (Ketcham and Carlson 2001) and subsequently lead to misidentification of shale components (Fig. 2b). Due to complex intrinsic shale properties, advanced image processing such as artefact removal and multiband thresholding is necessary (Iassonov et al. 2009). In addition to acquiring high-quality images, the images were cropped to remove edges artefacts present (Fig. 2b). This will also reduce the processing power requirements of the computing machine.

## Polygon surface representation: marching cubes algorithm

The marching cube algorithm is applied to reconstruct a surface CT scan volumetric data sets. It generates inter-slice connectivity, through linear interpolation of the scan line to calculate triangle vertices. The algorithm works in two steps: initially the surface that matches to the user-specified value is located, and subsequently, the normal to the surface of each vertex of each triangle is calculated.

The algorithm basically operates by determining how the surface intersects one cube and then proceeds (marches) to the subsequent cube. Full details of the implementation of the marching cube algorithm can be found in Lorensen and Cline (1987). In summary, the marching cubes algorithm for creating a polygonal surface representation of an isosurface is given as follows:

- Four slices are read into the memory;
- Two slices are scanned to form a cube from neighbouring slices, and four others from the next slice;
- By contrasting between the density values of the eight vertices while keeping the surface constant, the cube index can be calculated;
- From the index value obtained, using a pre-calculated table, a list of edges is determined;
- Applying linear interpolation to find the surface edge intersection using the densities at each vertex edge;
- The unit normal at each cube vertex is calculated using central difference approach and interpolating the normal to each triangle vertex by using the equations:

where D(i,j,k) is the density at pixel (i, j) in slice k and Δx,Δy,Δz are the lengths of the cube edges. When the gradient is divided by its length, it produces the unit normal at the vertex required for rendering.

7. Output the vertex normal and triangle vertices.

### Representative elementary volume (REV) approach

In pore-scale modelling, smaller subsamples can be treated as REV in order to be used in the evaluation of Darcy fluxes. REV is assumed to be obtained when the computed variable (in this case porosity) plotted against increasing sample size (Fig. 2d) does not change significantly at a plateau value (e.g. Figure 1). This is particularly useful when determining what sample to be used in modelling and simulating fluid flow through them.

In order to evaluate deterministic REV, each sample was further incrementally cropped into seven subsamples using corner point reference technique (see Fig. 2d). From the origin (at x, y, z = 0), we extract seven (7) three-dimensional subsample image windows of 50, 100, 200, 300, 400, 500 and 600 μm pixel lengths. The length of each subsample is denoted by l_{x} where x∈{1,2,3,4,…,n} and totallength,L=∑^{n}_{x=1}l_{x} (Fig. 2d). Similar REV data to the conceptual diagram shown in Fig. 1 were obtained for each sample, for porosity variable. REV was determined using the approaches of Costanza-Robinson et al. (2011) and was taken as the minimum window length scale (l) at which the absolute value of the relative gradient error (ε_{g}) in the porosity (∅) remained below 0.003:

where l is the window increment identity/number and δ is the magnitude of l_{x}. In physical terms, the relative gradient error REV criterion requires changes in the measured variable over a given length-scale increment to be relatively small proportional to the increment size.

The choice of the ε criterion (e.g. <0.002) depends on region II analysis given by (Bear 1972). Costanza-Robinson (2011) states that the semi-quantitative ε approach used here is advantageous because it makes REV estimation to be automated and reproducible across numerous images. This allows quantitative relationships between REV and other variables to be readily evaluated; however, it is not applicable to shale rocks which are highly heterogeneous.

### Image processing and heterogeneity analysis technique

In order to investigate the degree of sample heterogeneity, each image sample was subdivided into 1 mm × 1 mm × 1 mm as shown in Fig. 2e. Hence, x, y, z represent row (R), layer (L), and column number (N), respectively. Pore volume was measured for each of the subcropped samples after processing (Fig. 3).

**Fig. 3** Illustration of calculation of absolute value of the relative gradient error.

Each of the subcropped images was morphologically processed. A morphological filtering algorithm opens up holes and gaps in a mask to get rid of small and potentially spurious features. This is to enhance component segmentation (foreground from background as the case may be, Fig. 4a, b).

**Fig. 4** Micro-CT image of Eagle Ford shale rock sample obtained parallel to the bedding plane at resolution of 6.7 µm and physical extent is 1 mm × 1 mm. a 2D view of typical subsample before morphological processing. b 2D view of typical subsample after morphological processing.

In order to estimate physical properties and identify pore space as well as other components of each of the Eagle Ford parallel shale rock samples, manual greyscale thresholding was used to segment the images and to distinguish all its entities. Thresholding is a common technique used in segmentation; it assigns all the pixels that belong to the object based on their grey colour or brightness (represented by grey values) into individual group (Young et al. 1998). This option selects a window of greyscale values and is useful where segmentation can be achieved based on greyscale intensities.

Only pixels that have a greyscale value within the lower value and upper value were included in segmentation created by the thresholding tool. Starting with the darkest sections (i.e. the lowest density material, typically pores) to the brightest sections (i.e. highest density material, typically shale grains, organic matter and mineral components), eight (8) components were segmented. These eight (8) entities were identified with the names pores, organic matter, Shale grain 1 (SG1), Shale grain 2 (SG2), Shale grain 3 (SG3), Shale grain 4 (SG4), Shale grain 5 (SG5) and Minerals. The range of greyscale values used in the thresholding of the main sample is shown in histogram (Fig. 5).

**Fig. 5** A typical histogram showing the thresholding greyscale values of identified components.

However, the major mineral identified is pyrite because of its high density compared to other possible constituent minerals in a typical shale rock sample. Its identification is also in agreement with the previous research (Drillskill et al. 2013; Joe et al. 2007). Each entity surface is captured using marching cube algorithm (see “Polygon surface representation: Marching cubes algorithm” section for further description of the marching cube algorithm).

Marching cube is used to capture an isosurface by a divide-and-conquer approach of a region of space into 3D. The cube is created logically from eight pixels/vertices of a cell. Each of the vertices is assigned a value. The value at each vertices of each box is compared to the designated minimum value, and any value less than or equal to the minimum is inside the surface (Lorensen and Cline 1987). Figure 6 shows 3D and 2D of pores, organic matter and minerals for subsample L_{1}R_{1}N_{1}.

**Fig. 6** Micro-CT image of Eagle Ford shale rock subsample L_{1}R_{1}N_{1} at resolution of 6.7 µm and physical extent of 1 mm × 1 mm × 1 mm. Subsample with a segmented three-dimensional CT images showing pores (red), organic matter(turquoise) and minerals (green), e.g. pyrite. b Segmented two-dimensional CT images showing pores (red), organic matter (turquoise) and minerals (green), e.g. pyrite.

## Statistical tools for the analysis of pore geometry

Descriptive analyses were employed in order to estimate statistical measure of central tendency and dispersion necessary to organise and summarise pore information.

### Mean and standard deviation of pore volume

The mean of pores is computed using Eq. 5. Pores volume was labelled in each sample as x_{1},x_{2},………,x_{N},

where x_{1} is the first pore volume, x_{2} is the second and so on until the last pore which is x_{N}, where N is the population size. Hence, the population mean µ can be calculated as (Keller 2014):

The standard deviation σ is thus calculated to measure variability of the pore volume in the sample:

### Measure of dispersion: population skewness and kurtosis

For further characterisation of the pore volume distribution data, the population skewness (b_{1}) is thus calculated:

and kurtosis G_{2} is calculated as:

where k_{2} is the unbiased estimate of the second cumulant (i.e. sample variance) (Joanes and Gill 1998).

Unlike skewness which was used to measure lack of symmetry of the frequency curve of the pore volume distribution, we used kurtosis to evaluate the degree of flatness of pore volume distribution near its centre. This will be the extent to which the sample pore volume tails, lighter or heavier, in relation to a normal distribution. Sample pore volume data may show heavy tail (long tails) or light tail (short tails) depending on the presence or lack of microcracks. Hence, kurtosis could be classified as platykurtic, mesokurtic and leptokurtic when the value is less than 3, equal to 3 and greater than 3, respectively.

## Results and discussion

### Pore volume characterisation

In this section, we compare two algorithms for the computation of pore volume in order to determine the accuracy of the applied methods to be adopted in the computation of the properties on all the samples. The two methods available for pore volume computation are voxel method (VOX) and object-oriented bounding box (OBB). Voxel method is computed using the equation:

where Vxel is pore volume, f is the frequency and η is the voxel size, while the OBB method is computed using marching cube algorithm presented in “Statistical tools for the analysis of pore geometry” section.

In order to further evaluate the degree of accuracy of the two numerical algorithms, pore volume computation using simulation-based OBB and simulation-based VOX was compared with analytically calculated method (ANA). The computation of the pore volumes using ANA is similar to the VOX method, but, unlike the VOX method where the voxels were selected for computation numerically, the voxels count in ANA method was identified manually.

Figures 7 and 8 show the pore volume distribution for subsamples L_{2}R_{2}N_{3}, L_{3}R_{2}N_{3}, L_{5}R_{2}N_{3} and L_{6}R_{2}N_{3} using VOX and OBB methods. It can be seen that small pore volumes have higher frequency while big pore volumes have lower frequency. This observation is consistent in both methods. However, comparing the two methods, it is evident that OBB method over-predicted the number of pores with reference to a specific pore volume. For instance, for the class of pores with volume of 8000–8500 µm^{3} as shown in both Figs. 7 and 8, OBB over-predicted by up to a factor of 5 compared to the VOX method. This overestimation can be attributed to the result of marching cube algorithm over-circumscribing the pore volume especially as the voxel sizes increase.

**Fig. 7** Pore volume distribution of subsamples L2R2N3 (green), L3R2N3 (pink), L5R2N3 (black) and L6R2N3 (cyan) obtained from voxel method.

**Fig. 8** Pore volume distribution of subsamples L2R2N3 (green), L3R2N3 (pink), L5R2N3 (black) and L6R2N3 (cyan) obtained from object-oriented bounding box method.

Similar analysis was conducted to evaluate the degree of accuracy for the computation of three (3) entities: pore volume fraction (Fig. 9), organic matter volume fraction (Fig. 10) and minerals volume fraction (Fig. 11). For this purpose, the total volume of voxel was plotted against voxel number for each of the three (3) entities for the subsample L1R1N1. Figure 9 shows that the results of pore volume computation were the same and consistent for the simulation-based pore volume OBB, simulation-based pore volume VOX and the analytically calculated method ANA only when the voxel number is not more than 3.

**Fig. 9** Comparison of methods of computing pore volume fraction using ANA, VOX and OBB methods.

However, appreciable deviation of simulation-based pore volume (from OBB method) was observed from voxel number 3 and above. For instance, when voxel number is 95, simulation-based pore volume result was about seven (7) times more than the other two methods. Similar variation is observed in Figs. 10 and 11; however, the difference was observed much earlier starting from voxel number 3 (Fig. 10).

**Fig. 10** Comparison of methods of computing organic matter volume fraction using ANA, VOX and OBB methods.

**Fig. 11** Comparison of methods of computing minerals volume fraction using ANA, VOX and OBB methods.

In general, Table 1 summarises the comparison of the three methods explained above for the computation of volume fraction of pores (pore volume), organic matter and minerals using the total volume of voxel. It shows that object-oriented bounding box method always over-predicts bigger entity. From this investigation, it is therefore evident that the two numerical algorithms diverge in cases where the sample is much more in mineral components. To this end, the descriptive statistics will be done with simulation-based pore volume voxel (VOX) method.

**Table 1** Comparison of methods of computing volume fraction of pores, organic matter and minerals using voxel and object-oriented bounding box methods.

Histogram of the pore volume distribution for subsamples L1R4N4, L2R1N1, L2R1N2, L2R2N2, L4R1N3, L4R2N3 and L6R1N1 is shown in Figs. 12, 13 and 14. They all have exponential decay trend similar to pore-size distribution presented by (Curtis et al. 2012; Chen et al. 2013). The results of subsample analysis for Eagle Ford parallel indicated varying range of pore volume distribution.