## Abstract

Fractures are crucial for unconventional hydrocarbon exploitation, but it is difficult to accurately observe the 3D spatial distribution characteristics of fractures. Microtomography (micro-CT) technology makes it possible to observe the 3D structures of fractures at micro-scale. In this study, micron-CT scanning is conducted on multiple mud-shale samples of source rocks in the Permian Lucaogou Formation, Junggar Basin. The Avizo® software is applied to process and segment the micron-CT images, so as to obtain the 3D fracture structure model inside rock core.

#### Authors:

#### Chao Qi^{a}, Xiaoqi Wang^{b}, Wei Wang^{a}, Jie Liu^{a,c}, Jincai Tuo^{d}, Keyu Liu^{b,e}

^{a}School of Earth Sciences and Engineering, Sun Yat-sen University, Guangdong, 510275, China. ^{b}Research Institute of Petroleum Exploration and Development, PetroChina, Beijing, 100083, China. ^{c}Guangdong Provincial Key Laboratory of Mineral Resources & Geological Processes, Guangdong, 510275, China. ^{d}Gansu Provincial Key Laboratory of Petroleum Resources / Key Laboratory of Petroleum Resources Research, Institute of Geology and Geophysics, Chinese. ^{e}Academy of Sciences, Gansu, 730000, China. School of Geosciences, China University of Petroleum (East China), Shandong, 266580, China

Received 16 April 2018 Accepted 6 August 2018

Therefore, the independently-developed CTSTA program is adopted to quantitatively describe the micro-fractures inside rock core, including fracture dimension, extension direction and extension scale. Meanwhile, this study summarizes the classification characteristics of fractures and their anisotropy. On this basis, the fractal dimensions of fractures can also be extracted.

Previous studies show that the geometric features of fractures have self-similarity at large and small scales, which can be described by exponential laws; and the fractal dimension is a typical exponent. Through the quantitative description or characterization of 3D fractures at micro-scale, the distribution characteristics of fractures in large scales could be known.

## Introduction

The exploration and development of unconventional resources is becoming increasingly attractive world-wide (Schmoker, 2002; Law, 2002; Jia et al., 2012; Zou et al., 2015; Wang, 2016; Zhao et al., 2016) since the success of shale gas/oil production in the USA. Unconventional oil and gas are formed and stored in micro-to nano-scale pores of shale, which has ultra-low porosity and permeability (Li et al., 2016; Liu et al., 2016; Jia et al., 2014). In recent years, the characteristics of pore structures of tight reservoirs at micro-scales have been studied by using different methods.

For example, a number of researchers (Clarkson et al., 2012; Yang et al., 2013; Yuan et al., 2016; Cao et al., 2016; Zhang et al., 2017) study the pore structures of shales qualitatively and quantitatively by means of electron microscopy, low-temperature nitrogen adsorption and high pressure mercury injection experiments. Dai et al. (2016) propose a new effective porosity of movable fluid based on the NMR (nuclear magnetic resonance) technology and quantitatively calibrated the effective porosity occupied by movable fluid. Bai et al. (2013) carry out multi-scale (nano-to micro-) three-dimensional (3D) CT scanning and characterize the distribution and texture of micro-scale pore throats in tight reservoirs.

Some scholars (Ambrose et al., 2010; Ma et al., 2014) obtain nano-scale pore images using focused ion beam (FIB) technology and characterize the reservoir rocks. In the development of low permeable reservoirs, hydraulic fracturing is often used to create artificial fractures to increase production capacity of tight reservoirs by interconnecting natural fractures (or cracks thereafter used in this paper) and enhancing fluid movement (Zou et al., 2015; Guo et al., 2014). More accurate characterization of cracks is thus of great significance in the development of unconventional oil and gas.

Currently more studies focus on the characterization of pore-structures of porous media and less on the characterization of fractures. Some related work of fracture characterization includes: (1) Lee et al. (2013) propose a method for extracting the width, length and direction of cracks automatically from digital images of cracks on the surface of concrete; (2) Arena et al. (2014) develop a new semi-automatic method to isolate single cracks by identifying branch points, and then carried out the statistical analysis of the geometrical characteristics such as length, average width, and aspect ratio of micro-cracks;

(3) Wang et al. (2012) introduce a method of recognizing cracks based on the characteristics of local grid of images and measuring the direction, length and width of cracks in tunnel lining; (4) Song (2011) use the concept of fitness of two opposite surfaces to describe the morphology of fractures based on 3D laser scanning data; (5) Guo et al. (2014) describe the propagation of hydraulic fractures by combining rock mechanical experiments, 3D acoustic emission (AE) recording and localization, and CT scanning technique;

(6) Zhou et al. (2004) calculate the fractal dimensions of fractures of core samples and analyze the correlation between fractal dimensions and fracture densities. However, the above studies mainly achieve two-dimensional (2D) characterization of cracks, even though the observation techniques used in some of these studies are 3D. In view of this, this paper attempts to identify individual cracks in a 3D network to achieve the characterization of cracks in 3D space.

The characterization of fractures requires 3D observations with high resolution. The development of microtomography (e.g. micro-CT) imaging technology over the past two decades makes it possible to obtain high-resolution 3D images of cracks at micro-scales. Previous studies of 2D fractures demonstrate that some characteristics of fractures can be described by exponential laws with fractures showing self-similarity at different scales (Bažant, 1997; Bour and Davy, 1997; Borodich, 1999; Davy et al., 2013). Using self-similarity, one can predict the characteristics of fractures at the scales that cannot be observed from the characteristics at the observed scales.

In this study, the 3D crack-structures of tight reservoirs rock (shale) from microtomographic images are analyzed. The position, size, shape and direction of individual cracks are characterized and the statistics of the characteristics of the cracks are given. Furthermore, the fractal dimensions of 3D cracks are also obtained.

## Sample and data

Samples are taken from a shale outcrop in the Jimusaer sag, Junggar Basin, the Permian Lucaogou Formation. Petrophysical experiments show that (1) dark shale contains average 19.96% total organic carbon (TOC); (2) the organic matter is immature and is dominated by Type I kerogen. In this study, 6 core plug samples with similar lengths are drilled from the same sample. Hydrocarbon generation and expulsion simulation experiments are performed on 5 of the 6 samples under high-temperature and pressure and cracks are generated.

The isothermal experimental temperature is fixed to be 350 °C and the pore fluid pressures are set to be 10 MPa, 20 MPa, 30 MPa, 40 MPa, 60 MPa and maintain for 48 h for different samples. The corresponding samples are named YP10, YP20, YP30, YP40 and YP60, respectively. The original sample that is not subject to any experiment is named YS. Micro-CT images are obtained by Xradia 510 Versa X-ray microscopy (XRM) scanner. Each dataset contains 2008 slices (images) and each image contains 2004 × 2048 pixels. The resolutions of the 6 samples range from 6.18 to 9.12 μm.

## Image processing and characterization

### Workflow

These samples are processed in the following steps as shown in Fig. 1: image loading, cropping, segmentation, island removal and crack isolation.

**Fig. 1.** Workflow for core sample processing.

(1) The first step is to use the Avizo® software package Fire version to read 2008 grayscale images of the sample. Typical 2D greyscale images and segmented sub-zone are shown in Fig. 2.

**Fig. 2.** Typical 2D greyscale images and segmented sub-zone. (a) Sample YP10, (b) Sample YP20, (c) Sample YP30, (d) Sample YP40, (e) Sample YP60, (f) Sample YS.

(2) Fig. 2 shows that the samples only partly occupy the 3D image space, thus, we crop the maximum volume containing only the sample information, as indicated in the rectangles in Fig. 2.

(3) To identify cracks from the greyscale images is a prerequisite of the analysis of cracks. A suitable threshold value should be carefully selected by using the threshold tool in Avizo® and then the binary label-data are exported for each sample. Voids (cracks or pores) are labelled as 1 and solid matrix as 0 in these binary data.

(4) Voids in these binary data may be real pores and cracks, low-density minerals, or noise in the images. We focus on the analysis of cracks in this study, thus it is necessary to remove the isolated small pores by using the functionality of “removing islands” in Avizo®.

(5) We check the intersection of cracks in 3D space. The values of voxels at the intersecting positions are modified to be 0 (solid matrix) manually in Avizo® slice by slice. Then the intersected cracks are separated, thus the 3D quantitative description of the individual cracks can be achieved.

Fig. 3 shows the 3D volume rendering of separated cracks of three samples. Some separate cracks are shown in the same color because of the limitation of color-map. Comparing Fig. 3(b) and (c), it can be seen that the function of removing islands is to filter out some very small non-cracked structures.

**Fig. 3.** 3D volume rendering of the fractures of three typical core samples. (a) Sample YP10, (b) Sample YP30 (islands removed), (c) Sample YP30 (before islands removal), (d) Sample YP40.

### Methods for quantitative analysis

Using the processed binary data, we perform the characterization of the micro-cracks by running the independently-developed CTSTA code (Liu et al., 2013), which is available to the public (Liu et al., 2016) now. The characterization is mainly based on the percolation theory and the concept of cluster is used. Parameters that can be extracted from the analyzed volumes include: porosity, specific surface area, the connectivity of the volume, pore size distribution, and the position, size, shape, and anisotropy of each cluster. The definitions of parameters used in this paper are as follows:

(1) The specific surface area (SSA) s_{p}:

It is a parameter closely related to hydration property, weeping property and the properties of mechanics (Wilcox, 1990; Qiu et al., 1999).

(2) Cluster: A cluster is a group of voxels of the same material in the binary dataset connected by common planes. A cluster can be considered as a pore-structure that is not connected to other pores, or it can be a crack that does not intersect with other cracks.

(3) Percolation (or the connectivity of the volume): We define a direction as percolating when one boundary of the parallelepiped in this direction has the same cluster label as its opposite boundary. This means two boundaries of the direction are connected by the same cluster and the model is permeable. Fluid can only flow when the model is permeable. Thus the percolation is essential for fluid migration. The isolated pores and cracks in the model have no effect on the migration of the fluid. Knowing the percolating cluster(s), the effective porosity of the model can be calculated.

(4) Shape and orientation of clusters: A vector a_{i}=(a_{xi},a_{yi},a_{zi})^{T} is defined from the center of the cluster to the position of every voxel. Then the shape and orientation of a cluster can be described approximately by a tensor T define as:

The matrix T has 3 eigenvalues τ_{1}< τ_{2} <τ_{3} and corresponding eigenvector u_{1}, u_{2}, u_{3}. The matrix can be visually represented by an ellipsoid. The half lengths of the three principal axes of the ellipsoid are corresponding to the square roots of three eigenvalues of the clusters. In this paper, we describe cracks using clusters, and an equivalent ellipsoid of a crack is shown in Fig. 4.

**Fig. 4.** The equivalent ellipsoid for fractures.

(5) The anisotropy of a cluster is described by using a parameter named isotropy index which is defined as:

It is the ratio of the minimum and maximum eigenvalues of the cluster, and its value is between 0 and 1. The isotropy index of 1 indicates that the structure is spherical and isotropic. For the cracks studied in this paper, the isotropy index is the ratio of the thickness (eigenvalue (τ_{1})^{1/2}) to the length (eigenvalue (τ_{3})^{ 1/2}) of crack. The isotropy index of ideal Griffith crack should be close to 0 (minimum semi-axis of the ellipsoid tends to 0) (Griffith, 1921), but in the realistic microtomographic images, the isotropy index always is greater than 0 because of thickness, curvature and roughness of cracks. In our analysis, clusters with I≤0.2 are classified as a typical shape of cracks, while I≥0.2 indicates the crack is curved or does not have a typical shape.

(6) Orientation of crack plane: It is expressed by the normal direction of the crack surface, which coincides with the orientation of the eigenvector µ_{1.}

Note that the orientation of the crack surface obtained in this analysis is relative to the coordinate system of the analyzed parallelepiped, we do not re-orient cracks to the original directions in the field.

(7) Crack aperture distribution: For a plane parallel to x-axis, for example, we count the cracks in different sizes in the x direction on the plane. After counting other planes parallel to x-axis, we obtain the counting summation of different apertures in the x-direction. In the counting in the y- and z-directions the operations are the same.

(8) Fractal dimension: The most commonly used method, box-counting method (Russell et al., 1980), is used to calculate the fractal dimension which is derived from

The basic implementations cover the binary image with cubic boxes of length δ (regardless of the border effect) and counting the number of boxes that contain cracks. Using different δ values and by fitting the slope of LN(N(δ)) versus LN(δ) the fractal dimension of cracks is derived.

In the above quantitative parameters, porosity, specific surface area and the size of a cluster are all based on voxels, the results are precise values which are only affected by the resolution of the image. The shape and orientation of a cluster is a cumulative expression as shown in Eq. (2), which is not an accurate value, but an approximation of statistics. The larger the clusters, the higher the credibility of the statistical results of anisotropy and direction, as shown in Fig. 10 of reference (Liu et al., 2009). Clusters which have 3^{3} voxels are considered statistically significant. The larger the cluster, the better the morphological characterization.

## Results

We analyze all the 6 samples using the procedures described above. Here we only illustrate the results of Samples YP10, YP30 and YP40. Basic parameters of these three samples are listed in Table 1. In order to show the effect of “removing islands” on the quantitative analysis, the results of YP30 without “removing islands” are also shown.

**Table 1.** Basic characteristics of core samples.

Table 1 shows that the porosities of these models are less than 3%, and the lowest porosity is only 1.49%. The specific surface areas are relatively small. Our analyses show that (1) there is no percolating crack in Sample YP10, (2) Sample YP30 percolates in the x direction, and (3) Sample YP40 percolates in the x and z directions. Our analyses also show that “removing islands” can dramatically reduce the total number of clusters in the models, as shown in Table 1.

Before removing islands, the total number of clusters of Sample YP30 is much more than those of samples after the islands are removed. Analysis also resolves that the total numbers of small clusters with less than 10 voxels are 57342. The influence of removing islands on the porosity is about 0.2% for Sample YP30 and these isolated pores do not contribute to connectivity or fluid flow.

Statistics of crack dimensions along the coordinate axes in x, y and z directions are shown in Fig. 5. Fig. 5(a) shows that in all three directions, more than 90% of the openness of cracks are smaller than 10 voxels. There is no obvious difference among the three directions for Sample YP10. Fig. 5(b) shows that more than 90% of the openness of cracks are smaller than 19, 6 and 22 voxels in the x, y and z directions, respectively.

This means that the cracks in Sample YP30 are mainly perpendicular to the y direction, which is equivalent to the normal direction of crack-surfaces parallel to the y axis. Sample YP40 in Fig. 5(c) is similar to Sample YP30.

**Fig. 5.** Statistics of fracture dimensions along the coordinate axis. (a) Sample YP10, (b) Sample YP30 (with islands removed), (c) Sample YP40.

The statistical results of the orientation of crack planes with respect to the three coordinate axes are shown in Fig. 6. Fig. 6(a) shows that the normal direction of cracks in Sample YP10 is mostly parallel to the x-axis and perpendicular to the y and z directions. Fig. 6(b) and (c) shows that the normal direction of cracks in both Samples YP30 and YP40 is mostly parallel to the y direction and perpendicular to the x and z directions. It is quite clear that the concentration of the directions of cracks in Sample YP10 is higher than those of the other two samples. Combining with Fig. 2, Fig. 3, we know that most cracks are parallel to the bedding of the shale.