New pore space characterization method of shale matrix formation by considering organic and inorganic pores

Accepted Manuscript New Pore Space Characterization Method of Shale Matrix Formation by Considering Organic and Inorganic Pores Yongfei Yang, Jun Yao, Chenchen Wang, Ying Gao, Qi Zhang, Senyou An, Wenhui Song PII:

S1875-5100(15)30086-X

DOI:

10.1016/j.jngse.2015.08.017

Reference:

JNGSE 925

To appear in:

Journal of Natural Gas Science and Engineering

Received Date: 30 May 2015 Revised Date:

5 August 2015

Accepted Date: 6 August 2015

Please cite this article as: Yang, Y., Yao, J., Wang, C., Gao, Y., Zhang, Q., An, S., Song, W., New Pore Space Characterization Method of Shale Matrix Formation by Considering Organic and Inorganic Pores, Journal of Natural Gas Science & Engineering (2015), doi: 10.1016/j.jngse.2015.08.017. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

New Pore Space Characterization Method of Shale Matrix Formation by Considering Organic and Inorganic Pores Yongfei YANG a, Jun YAO a,*, Chenchen WANGb, Ying GAO a, Qi ZHANG a, Senyou AN a, Wenhui Song a

School of Petroleum Engineering, China University of Petroleum (East China), Qingdao, 266580, China b iRock Technologies Co., LTD., Beijing 10020, China

RI PT

a

Corresponding author. Tel.: +86 186 5321 6342, E-mail address: [email protected]

ABSTRACT

TE D

M AN U

SC

A shale matrix is too tight to be described using conventional methods, and digital core technology is becoming an alternative method. Because both organic and inorganic pores exist in the shale matrix, a digital core and a pore network model that could describe these two types of pores at the same time are constructed in this paper. Firstly, the inorganic pore digital core is constructed based on the multiple-point statistics method, and the organic pore digital core is constructed based on the Markov chain Monte Carlo method. The two types of digital cores are superposed together according to a superposition algorithm, which includes information about the shale organic and inorganic pores. The pore network models of different constructed digital cores are extracted using the pore space medial axis method. Finally, based on these platforms, the geometry and topology structure properties, the pore size distribution and the coordination number of a shale sample are analyzed. The results show that the pore size distribution of the shale sample generally ranges from 2 nm to 100 nm, mainly distributing from 5 nm to 20 nm. The coordination number is almost always in the range of 2 to 3. The digital core results match well with the experimental results to some extent for our study case.

Key words：shale; digital core; multiple-point statistics (MPS) method; Markov chain Monte Carlo (MCMC) method; superposition algorithm

AC C

EP

HIGHLIGHTS: A superposition algorithm is used to combine the inorganic pore digital core and the organic pore digital core. The multiple-point statistics method and the Markov chain Monte Carlo method are used to construct the inorganic pore digital core and the organic pore digital core, respectively. The pore space medial axis method is used to extract the pore network model from the digital core. The geometry and topology structure properties of the shale matrix sample are obtained from the pore network model and match well with experimental results.

1. Introduction

A complex shale pore structure includes nano organic pores in organic matter, micron-nano inorganic pores in inorganic mineral and micron-millimeter natural fractures (Wang and Reed, 2009). The shale pore is too small to be studied using conventional methods. At present, the methods commonly used to study the shale pore structure include a scanning electron microscope (SEM) (Clarkson et al., 2013), a focused ion beam scanning electron microscope (FIB-SEM) (Bai et al., 2013), atomic force microscopy (AFM) (Javadpour, 2009), nano CT, and transmission electron microscope scanning methods. Combined with an energy spectrum (ESD) or a backscatter 1

ACCEPTED MANUSCRIPT

M AN U

SC

RI PT

image (BEI), the three-dimensional distribution image of different mineral compositions can also be obtained. Based on the above equipment, the digital core images of shale can be captured. Digital core construction methods are divided into the following two classes: (1) the physical construction method in which focused ion beam scanning electron microscopy or nano CT can directly construct real 3D digital cores; and (2) the random reconstruct algorithm in which 2D scanning electron microscopy, atomic force microscopy, and a transmission electron microscope can only obtain two-dimensional image data, and digital cores are reconstructed using random reconstruction algorithms. The reconstruction methods commonly used include a Gauss simulation (Joshi, 1974), simulated annealing (Hazlett, 1997), the process-based simulation method (Bryant and Blunt, 1992), multiple point statistics (Okabe and Blunt, 2004), the sequential indicator simulation (Keehm, 2003), the Markov stochastic reconstruction method (Wu et al., 2006) and various hybrid methods (Hidajat et al., 2002; Liu et al., 2009; Okabe and Blunt, 2007; Yao et al., 2013). However, for either the physical construction method or the random reconstruction algorithm using 2D image data, only single-scale or single-component digital cores can be constructed, which are not accurate enough for the strong heterogeneity of shale (Yang et al., 2015). It is also very difficult to consider organic and inorganic pores simultaneously. Therefore, constructing digital cores containing both an organic pore and an inorganic pore is the focus of this paper. First, the inorganic pore digital core was constructed based on the multiple-point statistics method, and the organic pore digital core will be constructed based on the Markov chain Monte Carlo method. The two digital cores were combined together according to a superposition method. This superimposed digital core contains both inorganic pore information and organic pore information. Second, the pore network models were extracted from the superimposed digital core according to the medial axis thinning algorithm. Finally, based on the pore network models, the geometry and topology properties of the shale samples were analyzed and compared with the experimental results.

TE D

2. Superposition method to construct the digital core of shale

AC C

EP

A high resolution 2D image is often available using optical or electronic equipment. The 3D pore space can be reconstructed from 2D images using statistical methods by considering the porosity information, lineal path function, two-point autocorrelation function, etc. Traditional porosity and two-point statistical techniques usually cannot reproduce the long-range connectivity of the pore space for using low-order information, especially for tight porous media (Adler et al., 1992; Biswal et al., 1999; Hazlett, 1997; Levitz, 1998; Quiblier, 1984). To establish a complicated pore space with submicron structures, the multiple-point statistics method is used to capture the large-scale inorganic pore space, and the Markov chain Monte Carlo (MCMC) method is used to capture the small-scale organic pore space.

2.1. Multiple-point statistics (MPS) method for the inorganic pore space Multiple-point statistics comes from a geostatistical method to reconstruct large-scale patterns using pixel-based representations (Caers, 2001; Mariethoz and Caers, 2014). This method is a good way to reconstruct higher resolution 3D images from readily available 2D microscopic images. MPS can characterize the correlation between many points compared with traditional two-point geostatistics. The main difference between MPS and two-point geostatistics is that they determine the conditional probability differently. MPS uses training images instead of a variogram function to reproduce the geological spatial structure. We use the following steps to construct the pore space: (1) Use a template scanning two-dimensional training images and build a search tree. 2

ACCEPTED MANUSCRIPT

Suppose the data template is τ and the set template center is u , and the other positions on this template are

uα = u + hα ; α = 1, …, n , where hα is the vector from the center position pointing to the other positions. Therefore, the data template τ is composed of n locations uα and a central location u . Take Fig. 1 as an example. The 2D template size is 9×9, and uα is determined by the center element u and the 80 other

RI PT

elements.

TE D

M AN U

u

SC

uα

Fig. 1. A 9×9 template to capture multiple-point statistics. The template is used to scan the training image, and each data event of any possible pattern of the pore space and skeleton is recorded

{sk ; k = 1,…, m} . A data event

EP

Suppose an attribute S has m possible states

d (u ) composed of states at

different locations, shown in Fig. 2(c), is defined by

{

}

AC C

d (u ) = i (u )； i (uα )=i (u + hα ) = skα ; α = 1, … , n

where i (u + hα ) indicates the state value at position uα ;

(1)

d (u) indicates the state value of

i (u1 ),…, i (u n ) at n positions, which are sk1 ,…, skn . Fig. 2 shows the process of obtaining a data event using the template to the scan training image. The purpose of scanning the training image using the data template is to get statistical probability of a data event d (u ) .

{

}

Prob {d (u )} = Prob i (u )；i (uα ) = skα ; α = 1, … , n When a data event in the training image is the same as a data event

3

(2)

d (u) , we call it a replicate. The ratio

ACCEPTED MANUSCRIPT between the replicate number of data events n ( d (u) ) and the size of the effective training image (or eroded image) N n is the occurrence probability of this data event. N n is also the number of different center locations of the template over the training image.

{

}

n ( d (u) )

(3)

Nn

RI PT

Prob i(u)；i(uα ) = skα ;α = 1,…, n =

For any simulation point u , the conditional probability distribution function (CPDF) of attribute

i(u)

takes one value of m possible states and should be determined at the conditions that give n condition data values

2010)

{

}

Prob i (u) = sk d ( u ) =

}

Prob i (u) = sk and i (uα ) = skα ; α = 1,…, n

{

}

Prob i (uα ) = skα ; α = 1,…, n

}

M AN U

{

{

SC

i (uα ) . According to the Bayesian conditional probability formula, the CPDF can be expressed as (Zhang et al.,

where the denominator Prob i (uα ) = skα ; α = 1, … , n

(4)

is the occurrence probability of the data event,

and it can inferred by counting the number n ( d (u) ) of replicates of the conditioning data event

{

}

d (u ) = i (u )； i (uα )=i (u + hα ) = skα ; α = 1, … , n

{

in the associated effective training image (or eroded

}

TE D

training image); the numerator Prob i (u ) = sk and i (uα ) = skα ; α = 1, … , n

is the probability of both the

data event and the simulation point u taking some state value existing together, and it can be obtained by counting the number nk ( d (u ) ) of replicates, among the n ( d (u ) ) previous ones, associated with a central

EP

value d (u) equal to sk . The CPDF can also be expressed as

{

}

AC C

Prob i (u) = sk i (uα ) = skα ; α = 1,…, n =

nk ( d (u) ) n ( d (u) )

(5)

where nk ( d (u ) ) is a replicate number of state i (u) = sk in the existing replicate n ( d (u ) ) . For the porous medium, the simulation point can only take two types of states: pore or skeleton. The simulation point state (pore or skeleton) can be determined by the nearest n points’ state and the CPDF obtained from training images scanning. If the simulation of every node needs to scan the training image to get the CPDF, it will seriously affect the speed of the simulation. This paper uses a search tree data structure to accelerate the reconstruction process, and this search tree can be obtained by scanning the training image once. The simulated images can be generated directly from the CPDF mode "feature library," which stores all possible patterns in the search tree (Strebelle, 2002; Zhang et al., 2015). (2) Define a random path to visit all unsampled voxels. For every new voxel, obtain the data event using the 4

ACCEPTED MANUSCRIPT

M AN U

SC

RI PT

same template in step (1), and obtain the CPDF of this voxel from the search tree. Take the state value stochastically, and put it into conditional events for the future voxel simulation. Simulate the following voxels on the random path until a new 2D image is produced. (3) Take the new 2D image produced in step (2) as an example, repeat the above process, and another new image will be produced. Producing images one by one, N images will be obtained. Match every pixel in the 2D image with the 3D space voxels, combine the newly produced N images and the original training image, and obtain an N+1 layer digital core. The detailed flow diagram of the multiple-point statistics reconstruction of porous media is given in the literature (Okabe and Blunt, 2005).

(a) 5×5 template;

(b) 17×17 training image;

(c) a data event (or a pattern)

Fig. 2. A data event obtained from the scanning training image. The frequency of a data event (or a pattern) is obtained by

TE D

scanning the 5×5 template over the 17×17 training image.

2.2. Markov chain Monte Carlo (MCMC) method for the organic pore space

AC C

EP

The Markov chain Monte Carlo (MCMC) method includes a Markov chain and the Monte Carlo method. The Markov chain (MC) was proposed by Markov in 1906, and this theory is extensively used in engineering technology, economic management, natural science and other fields. The Monte Carlo method was first applied by Stanislaw Ulam et al. in the 1940s (Metropolis and Ulam, 1949), and it has been widely applied in mathematics, physics, statistics, etc. The basic process of Monte Carlo is to construct or describe the probability of a problem, sample from a known probability distribution, and estimate various statistics parameters. The MCMC method is developed from Markov random fields (MRFs), which are widely used for image processing research (Geman and Geman, 1984). The MCMC method was first used to reconstruct 2D thin section images of soil samples by Wu et al., and a 3D reconstruction of reservoir rocks was extended to generate pore space of real heterogeneous porous media (Wu et al., 2004; Wu et al., 2006). A Markov chain can be used for describing systems that follow a chain of linked events, where what happens next depends only on the current state of the system. Suppose X(t) is a known random process describing the state of a process at time t0. If the state of random process X(t) in the future time t(t>t0) depends only on the state at time t0 and has nothing to do with the state before time t0, this random process X(t) is called a Markov process. During a Markov process, the relationship between “PAST” and “FUTURE” can be described as “FUTURE” linked to the “PAST” through “NOW”, and if “NOW” is confirmed, “FUTURE” is irrelevant to the “PAST”. Suppose a system is composed of n voxels ( i = 1, … , n ), and X = X 1 , …, X n indicates the state of the voxels. Therefore, X i 5

ACCEPTED MANUSCRIPT means that site i has a state xi . Specifically, for a certain voxel s, Λ − s represents all of the other points except s. The neighbor of voxel s, N s , must exist and can be described as (Wu et al., 2006):

p( xs x(Λ− s )) ≈ p( xs x( Ns ))

(6)

RI PT

The probability of voxel s has a certain state, so the states of all of the other points except s can be approximated by the conditional probability equation (6), which only considers the neighbor of voxel s. Generally, the distribution in equation (6) is difficult to simulate, and it takes a long time to reach convergence with existing iterative methods. Thus, its practical application is limited, especially for high dimension applications. For the two-dimensional case, a new single-pass method is used to simulate the conditional probability (Qian and Titterington, 1991). We extend the two-dimensional method to a three-dimensional situation in this paper.

layers of a rectangular grid filled with cube voxels.

( i, j , k )

SC

Assume VLMN = {(l , m, n) : 0 < l ≤ L,0 < m ≤ M , 0 < n ≤ N } represents L rows, M columns, and N represents the intersected voxel of row i, column j

M AN U

and layer k, and its associated state is expressed as X ijk . V ijk represents a rectangular parallelepiped array of voxels, and its associated state is expressed as vector X (Vijk ) . Any (i, j , k ) ∈ VLMN has the joint probability function: i

j

k

p( x(Vijk )) = ∏∏∏ p( xlmn xl −1,m,n , xl , m−1,n , xl ,m,n−1 )

(7)

TE D

l = 2 m=2 n=2

The conditional probability of each voxel for its random Markov field is:

p( xijk {xlmn : (l , m, n) ≠ (i, j, k )}) = p( xijk {xlmn : (l, m, n) ∈ N (ijk )})

(8)

where N (ijk ) represents the neighbor of (i , j , k ) . For instance, for a 19 neighborhood model, the nearest

EP

neighbor of the target voxel is given by:

AC C

 (i − 1, j, k )  (i + 1, j, k )   (i, j − 1, k + 1) N19 (ijk ) =   (i, j + 1, k − 1) (i, j + 1, k + 1)   (i, j − 1, k − 1)

(i, j − 1, k ) (i, j + 1, k ) (i − 1, j , k + 1) (i + 1, j, k − 1) (i + 1, j, k + 1) (i − 1, j, k − 1)

(i, j , k − 1)  (i, j , k + 1)  (i − 1, j + 1, k )   (i + 1, j − 1, k )  (i + 1, j + 1, k )   (i − 1, j − 1, k ) 

(9)

In the absence of three-dimensional information, three mutually perpendicular independent two-dimensional images can be used to reconstruct a three-dimensional chain. The core slice data from three mutually perpendicular planes xy, yz and xz are extracted, and the corresponding binary images are obtained (state only 0 or 1). The three-dimensional Markov model is reconstructed using the above core slice data, namely, the three-dimensional digital core. To get more efficient computation, we choose the double voxel method, which simultaneously produces two new voxels, i.e., voxel (i, j, k) and (i, j +1, k), during the modeling process.

6

RI PT

ACCEPTED MANUSCRIPT

Fig. 3. Sketch map of a heterogeneous porous media simulation based on three perpendicular thin sections.

M AN U

SC

Based on the Markov chain Monte Carlo method, the reconstruction steps for the three-dimensional digital core are as follows: Firstly, the porosity of the horizontal core slices is used as the conditional probability. Whether the first voxel in the one-dimensional chain is porous or not can be determined by the porosity (Fig. 4a). Secondly, the voxels on the first line of the first layer are simulated along the y-direction. The second voxel of the first line is simulated using the 2-neighborhood model, and after that, the 3-neighborhood model is used at the beginning of the simulation of the third voxel. Its conditional probability is calculated by the two-dimensional core slice on the xy plane (Fig. 4b). Thirdly, the first layer was simulated along the x direction line by line. For double elements (i, j ) and (i , j +1) , 3 and 4 neighborhoods are used when simulating the edge voxels, and 5 and 6 neighborhoods are

TE D

used when simulating the internal voxels. The conditional probability is derived by the two-dimensional core slice on the xy plane (Fig. 4c). Lastly, each layer of voxels was simulated along the Z-direction, and the three-dimensional model was reconstructed. For double elements (i, j ) and (i , j +1) , the construction method of the first row on the second

EP

layer is similar to that of the second row on the first layer, i.e., 3 and 4 neighborhoods were used when simulating the edge voxels, and 5 and 6 neighborhoods were used when simulating the internal voxels. The only difference in the conditional probability is from the two-dimensional slices in the other direction. From the beginning of the second row on the second layer, 9 and 10 neighborhoods were used when simulating the edge voxels, and 14 and 15 neighborhoods were used when simulating the internal voxels (Fig. 4d). The 15 3D neighborhoods can be expressed as:

AC C

 (i − 2, j , k − 1) (i − 2, j + 1, k − 1) (i − 1, j , k − 1)  (i − 1, j + 1, k − 1) (i, j − 1, k − 1) (i, j, k − 1)   N15 (ijk ; i, j + 1, k ) =  (i, j + 1, k − 1) (i − 2, j, k ) (i − 2, j + 1, k )    (i − 1, j , k ) (i − 1, j + 1, k )   (i − 1, j − 1, k )  (i, j − 1, k ) 

7

(10)

SC

RI PT

ACCEPTED MANUSCRIPT

Fig. 4. The 3D digital core reconstruction process based on the MCMC method (white elements with lines are the

M AN U

"NOW" simulating voxel; white elements with dotted lines are the "FUTURE" simulating voxel; dark elements are the voxels already produced, while light elements are the neighbors of the voxel currently being simulated, and they both belong to "PAST" voxels)

TE D

Markov chain Monte Carlo (MCMC) is a special type of Monte Carlo method, and a Markov chain is applied in a random process to realize the Monte Carlo dynamic simulation. The MCMC method is a stochastic simulation method to achieve the simulation target distribution by an ergodic constraint based on the dynamic structure of the Markov chain. The MCMC method constructs the digital core from three directions, and the spatial distribution characteristics of the digital core are similar to the real core. The MCMC method is a random statistics method with very fast modeling speed, and the established digital core has good pore connectivity.

2.3. Superposition method for the shale total pore space

AC C

EP

The 3D digital core data are represented in binary 0 and 1, where 0 indicates the rock skeleton and 1 indicates the pore space. We use the following steps to reconstruct the shale digital core. First, scanning electron microscopy (SEM) is used to obtain two-dimensional images with different resolutions based on different types of cores, and image processing is used to get the two value images to describe the characteristics of the organic and inorganic shale pores. Second, the MPS method and the MCMC method are used to establish the inorganic and organic digital cores, respectively. Lastly, the superposition algorithm is proposed to reconstruct the shale gas reservoir digital rock considering the organic and inorganic pores. The constructed shale digital rock, which retains the inorganic and organic pore geometry and topology information, can describe the pore distribution characteristics of shale rocks. The superposition algorithm has been successfully applied to characterize a heterogeneous carbonate reservoir (Yao et al., 2013). We will extend this method to construct the digital core for shale samples. The inorganic pore digital core and organic pore digital core have the same physical size and different resolutions (i is the resolution ratio between the low resolution pore digital rock and the high resolution pore digital rock). The low resolution pore digital rock is refined into i*i*i voxels (Fig. 5).

8

i

ACCEPTED MANUSCRIPT

i

i i

RI PT

i

i

Fig. 5. Sketch map of voxel refinement for the low resolution digital rock (i=5, a black box is one pore space element, a

SC

white box is one rock matrix element, and they are refined into 125 (5*5*5) elements.)

After the refinement, the detailed steps of the superposition operations between two digital rocks are as follows:

M AN U

Ω = Ω1 + Ω2

Ω skeleton = Ω1skeleton + Ω 2 skeleton  organic pore = Ω1skeleton + Ω 2 organic pore Ω  inorganic pore = Ω1 pore + Ω 2 skeleton Ω Ωinorganic pore = Ω pore + Ω pore  1 2

0 = 0+0 1 = 0 +1 2 = 1+ 0 2 = 1+1

AC C

EP

TE D

where Ω indicates the shale matrix digital rock, Ω1 denotes the inorganic pore digital rock and Ω2 denotes the organic pore digital rock. For the inorganic pore digital core, 0 means skeleton and 1 means inorganic pore; for the organic pore digital core, 0 means skeleton and 1 means organic pore; and for the superposed shale matrix digital core, 0 represents the skeleton, 1 represents the organic pore and 2 represents the inorganic pore to show the difference. When the organic pore meets the skeleton, the superposed results become porous to incorporate the organic pore space in the shale matrix, so the corresponding data operations are: 0+0=0, 0+1=1, 1+0=2, and 1+1=2 (Wang et al., 2014). Taking 1+0=2 as an example, the 1 means an inorganic pore in an inorganic pore digital rock, the 0 means a skeleton in an organic pore digital rock, and the 2 means an inorganic pore in a shale matrix digital rock.

3. Case study and parameter analysis 3.1. Image segmentation

Based on scanning electron microscopy (SEM), a low resolution inorganic pore scanning image and a high resolution organic pore scanning image are obtained. For the further superposing digital cores, the same physical size images are selected (Fig. 6). The pixel size of the inorganic pore image and the organic pore image are both 411×326, the inorganic pore image resolution is 0.022 µm per pixel, and the organic pore image resolution is 0.0044 µm per pixel. Therefore, the image resolution ratio of the inorganic and organic pore images is 5:1. The corresponding segmentation binary images can be obtained using the maximum class separation distance method proposed by Otsu (Otsu, 1979). For the segmentation binary images, white is the rock skeleton and black is the 9

ACCEPTED MANUSCRIPT

SC

RI PT

rock pore space (Fig. 7). Taking pixels (pixel) as the basic unit, the pixel value is set to 0 when the pixel is located in the shale skeleton and the pixel value is set to 1 when the pixel is located in the shale pore space. Therefore, the corresponding shale core samples can be characterized by a series of 0 and 1 data for future superposition.

(a) Inorganic pore SEM image

(b) Organic pore SEM image

EP

TE D

M AN U

Fig. 6. Typical inorganic and organic pore SEM images in shale rock

(a) Inorganic pore binary image

(b) Organic pore binary image

AC C

Fig. 7. Binary image of inorganic and organic pore in shale rocks (white is the rock matrix, and black is the pore space)

3.2. Digital core construction and pore network model extraction Based on the MPS method and the segmented inorganic pore SEM image in Fig. 7(a), the shale inorganic pore digital core is constructed and shown in Fig. 8(a). Based on the Markov chain Monte Carlo (MCMC) method and the segmented organic pore SEM image in Fig. 7(b), the shale organic pore digital core is constructed and shown in Fig. 8(b). The superposition digital core is shown in Fig. 8 (c) after being superimposed using the described algorithm in Section 2.3.

10

ACCEPTED MANUSCRIPT

RI PT

(a) Inorganic digital core (b) Organic digital core (c) Superposed digital core Fig. 8. Inorganic digital core, organic digital core and superposed digital core of shale rock

EP

TE D

M AN U

SC

It is usually difficult to analyze the geometric characteristics from digital cores with irregularly shaped pore space. Therefore, the corresponding pore network models are extracted using a medial axis thinning algorithm based on the reconstructed digital cores (Lee et al., 1994). The medial axis of the digital core is composed of curves formed by connected pixels in the center position of the pore space (Lindquist et al., 1996). The main steps are as follows. First, based on image analysis theory, the isolated pore and rock skeleton from image segmentation errors are removed. Second, based on the medial axis thinning algorithm, the pore space’s medial axis system is established, and the redundant structure caused by the rough surface of the rock wall is deleted to ensure the accuracy of the topological structure of the pore system. Third, the pore center is positioned, and the pore throat is determined. Finally, the geometric parameters of the pores and throats, including the volume, the radius of the inscribed circle, the length, and the shape factor, are measured. A pore network model with real topology and geometry is established. Using the above steps, the inorganic pore network model, the organic pore network model and the superposed pore network model extracted from the corresponding digital cores in Fig. 8 are shown in Fig. 9.

AC C

(a) Inorganic pore network model (b) Organic pore network model (c) Superposed pore network model Fig. 9. Inorganic pore network model, organic pore network model and superposed pore network model of shale rock (green is pore, red is throat)

3.3. Parameter analysis and validation Based on the two research platforms, the digital core and the pore network model, the geometrical and topological features of the shale pore space can be evaluated. The pore radius size of the organic shale pore space is 2 to 22 nm, and the pore radius size of the inorganic shale pore space is 12 to 112 nm. The superposed digital rock and the pore network model can have a similar pore size distribution, generally ranging from 2 nm to 100 nm and mainly distributed from 5 nm to 20 nm (Fig. 10). In the superposed digital core, the inorganic pore numbers take a very small percentage of total inorganic and organic pores. Therefore, the organic pore characteristics have a greater influence in the pore radius distribution curve. Fig. 11 shows that the coordination number is mostly in the range of 2 to 3 and that the consistency for three digital cores is good. 11

ACCEPTED MANUSCRIPT 0.4 0.35 Inorganic digital core

0.3

superposition digital core

0.2

RI PT

Frequency

Organic digital core 0.25

0.15 0.1 0.05 0 0.02

0.04

0.06 Radius (micron)

0.08

0.1

0.12

SC

0

M AN U

Fig. 10. Comparison of the pore size distribution obtained from the inorganic, organic and superposed digital cores of shale rocks

0.45 0.4

Inorganic digital core

0.35

Organic digital core Superposition digital core

0.25

TE D

Frequency

0.3

0.2 0.15

0.1

EP

0.05 0

1

2

3

4

5 6 7 8 Coordination Number

9

10

11

12

AC C

Fig. 11. Comparison of the coordination number obtained from inorganic, organic and superposed digital cores of shale rocks

12

ACCEPTED MANUSCRIPT

RI PT

The nitrogen adsorption method is currently used to characterize the shale micro pore structure. The basic theory is to calculate the pore structure parameters through the shale adsorption and desorption isotherms obtained by nitrogen adsorption experiments (Barrett et al., 1951). We compared the pore size distribution between the digital core analysis results and the nitrogen adsorption experimental results, and there is a reasonable match, as shown in Fig. 12. For the future, more cases should be studied to test the applicability, and new methods are needed for validation. 0.4 0.35

Superposition digital core

0.25

SC

experiment

0.2 0.15 0.1 0.05 0 0

0.02

M AN U

Frequency

0.3

0.04

0.06 0.08 Radius (micron)

0.1

0.12

Fig. 12. Comparison of the pore size distribution between the superposed digital core and the experiment

4. Conclusions

AC C

EP

TE D

(1) Based on the multiple point statistics method (MPS), the inorganic pore digital core was established. The MPS method can predict the long-range connectivity, so the large-scale characterizations of inorganic pore can be studied. Additionally, based on the Markov chain Monte Carlo (MCMC) method, the organic pore digital core was established. The MCMC method has high calculation speed, so it can save lots of reconstruction time for a small-scale pore space. (2) Using the proposed superposition coupling algorithm, the superposed digital core was established considering both the inorganic pore space and the organic pore space. They can take full advantages of the two methods and consider scale at the same time. (3) The corresponding pore network models were extracted based on the pore space medial axis method. The digital core and pore network models of the shale sample were both established, and the geometry and topology of the sample were analyzed. (4) The pore size distribution of the shale sample generally ranges from 2 nm to 100 nm and is mainly distributed between 5 nm to 20 nm. The coordination number is mostly in the range of 2 to 3, and the consistency for the inorganic digital core, the organic digital core and the superposed digital core is good. Comparing the pore radius distribution obtained from the pore network model with experiments shows a very good match, which verifies the feasibility of the construction methods to some extent.

13

ACCEPTED MANUSCRIPT

RI PT

Acknowledgements: We would like to express appreciation for the following financial support: the National Natural Science Foundation of China (No. 51304232, 51490654, 51234007), the National Basic Research Program of China (No. 2014CB239103), the Major Programs of the Ministry of Education of China (No. 311009), the Specialized Research Fund for the Doctoral Program of Higher Education (20120133120017), the Fundamental Research Funds for the Central Universities (No. 14CX05026A, 14CX06090A), Introducing Talents of Discipline to Universities (B08028), and the Program for Changjiang Scholars and Innovative Research Team in University (IRT1294).

References

Adler, P.M., Jacquin, C.G. and Thovert, J.-F., 1992. The formation factor of reconstructed porous media.

SC

Water Resources Research, 28(6): 1571-1576.

Bai, B., Elgmati, M., Zhang, H. and Wei, M., 2013. Rock characterization of Fayetteville shale gas plays. Fuel, 105(1): 645-652.

M AN U

Barrett, E.P., Joyner, L.G. and Halenda, P.P., 1951. The Determination of Pore Volume and Area Distributions in Porous Substances. I. Computations from Nitrogen Isotherms. Journal of the American Chemical Society, 73(1): 373-380.

Biswal, B., Manwart, C., Hilfer, R., Bakke, S. and Øren, P.E., 1999. Quantitative analysis of experimental and synthetic microstructures for sedimentary rock. Physica A: Statistical Mechanics and its Applications, 273(3-4): 452-475.

Bryant, S. and Blunt, M., 1992. Prediction of relative permeability in simple porous media. Physical Review A, 46(4): 2004-2011.

TE D

Caers, J., 2001. Geostatistical reservoir modelling using statistical pattern recognition. Journal of Petroleum Science and Engineering, 29(3-4): 177-188. Clarkson, C.R. et al., 2013. Pore structure characterization of North American shale gas reservoirs using USANS/SANS, gas adsorption, and mercury intrusion. Fuel, 103(1): 606-616. Geman, S. and Geman, D., 1984. Stochastic relaxation, gibbs distributions, and the bayesian restoration

EP

of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, PAMI-6(6): 721-741.

Hazlett, R.D., 1997. Satistical characterization and stochastic modeling of pore networks in relation to

AC C

fluid flow. Mathematical Geology, 29(6): 801-822.

Hidajat, I., Rastogi, A., Singh, M. and Mohanty, K.K., 2002. Transport properties of porous media reconstructed from thin-sections. SPE Journal, 7(1): 40-48.

Javadpour, F., 2009. Nanopores and Apparent Permeability of Gas Flow in Mudrocks (Shales and Siltstone). Journal of Canadian Petroleum Technology, 48(8): 16-21.

Joshi, M., 1974. A class of stochastic models for porous media, University of Kansas, Kansas, 163 pp. Keehm, Y., 2003. Computational rock physics: Transport properties in porous media and applications. 3085197 Thesis, Stanford University, United States -- California, 135 pp. Lee, T.C., Kashyap, R.L. and Chu, C.N., 1994. Building Skeleton Models via 3-D Medial Surface Axis Thinning Algorithms. CVGIP: Graphical Models and Image Processing, 56(6): 462-478. Levitz, P., 1998. Off-lattice reconstruction of porous media: Critical evaluation, geometrical confinement and molecular transport. Advances in Colloid and Interface Science, 76-77(1): 14

ACCEPTED MANUSCRIPT 71-106. Lindquist, W.B., Lee, S.M., Coker, D.A., Jones, K.W. and Spanne, P., 1996. Medial axis analysis of void structure in three-dimensional tomographic images of porous media. Journal of Geophysical Research, 101(B4): 8297-8310. Liu, X., Sun, J. and Wang, H., 2009. Reconstruction of 3-D digital cores using a hybrid method. Applied Geophysics, 6(2): 105-112. Images. Wiley-Blackwell.

RI PT

Mariethoz, G. and Caers, J., 2014. Multiple-point Geostatistics: Stochastic Modeling with Training Metropolis, N. and Ulam, S., 1949. The Monte Carlo Method. Journal of the American Statistical Association, 44(247): 335-341.

Okabe, H. and Blunt, M.J., 2004. Prediction of permeability for porous media reconstructed using multiple-point statistics. Physical Review E, 70(6): 066135. Petroleum Science and Engineering, 46(1-2): 121-137.

SC

Okabe, H. and Blunt, M.J., 2005. Pore space reconstruction using multiple-point statistics. Journal of Okabe, H. and Blunt, M.J., 2007. Pore space reconstruction of vuggy carbonates using microtomography and multiple-point statistics. Water Resources Research, 43(12): W12S02.

M AN U

Otsu, N., 1979. A threshold selection method from gray-level histograms. IEEE Transactions on Systems, Man and Cybernetics, 9(1): 62-66.

Qian, W. and Titterington, D.M., 1991. Pixel labelling for three-dimensional scenes based on Markov mesh models. Signal Processing, 22(3): 313-328.

Quiblier, J.A., 1984. A new three-dimensional modeling technique for studying porous media. Journal of Colloid and Interface Science, 98(1): 84-102.

Strebelle, S., 2002. Conditional simulation of complex geological structures using multiple-point

TE D

statistics. Mathematical Geology, 34(1): 1-21.

Wang, C. et al., 2014. Organic and inorganic pore structure analysis in shale matrix with superposition method, Unconventional Resources Technology Conference. Society of Petroleum Engineers, Denver, Colorado, USA.

Wang, F.P. and Reed, R.M., 2009. Pore Networks and Fluid Flow in Gas Shales, SPE Annual Technical

EP

Conference and Exhibition. Society of Petroleum Engineers, New Orleans, Louisiana. Wu, K.J., Nunan, N., Crawford, J.W., Young, I.M. and Ritz, K., 2004. An efficient Markov chain model for the simulation of heterogeneous soil structure. Soil Science Society of America Journal,

AC C

68(2): 346-351.

Wu, K.J. et al., 2006. 3D stochastic modelling of heterogeneous porous media - Applications to reservoir rocks. Transport In Porous Media, 65(3): 443-467.

Yang, Y., Wang, C., Yao, J. and Gao, Y., 2015. A new voxel upscaling method based on digital rock. International Journal for Multiscale Computational Engineering, 13(4).

Yao, J., Wang, C., Yang, Y., Hu, R. and Wang, X., 2013. The construction of carbonate digital rock with hybrid superposition method. Journal of Petroleum Science and Engineering, 110(1): 263-267. Zhang, T., Du, Y., Huang, T. and Li, X., 2015. Reconstruction of porous media using multiple-point statistics with data conditioning. Stochastic Environmental Research and Risk Assessment, 29(3): 727-738. Zhang, T., Li, D., Lu, D. and Yang, J., 2010. Research on the reconstruction method of porous media using multiple-point geostatistics. Science China Physics, Mechanics and Astronomy, 53(1): 122-134. 15