Physical Transport Properties of Porous Rock with Computed Tomography

Abstract

In this chapter, three-dimensional digital rock models can be constructed by the micron X-ray computed tomography (CT). Then, lattice gas automata was applied to simulate the flow of electrical current in the saturated digital rocks to reveal the non-Archie relation of resistivity index and water saturation (I-Sw). The flow of single-phase Newtonian fluid in pore space had been studied with LBM for calculating the absolute permeability. Moreover, we have developed a model based on digital rock to simulate thermal neutrons transporting for imaging the anisotropy of pore structure. The advantages of the model over traditional methods indicate that it can simultaneously consider both the separation of matrix and pore and the distribution of mineral components. The results of numerical simulation with Monte Carlo are in good agreement with the pore distribution from X-ray CT, which can further verify the validity of the new model. In contrast to the conventional conclusion, we find that the porosity calculated with neutron data can be affected by the anisotropy. Therefore, a new formula to relate the resolution of array detectors to the quality of imaging, had been proposed to analyze the critical resolution and to optimize the number of neutrons in each simulation.

1. Introduction

Due to extremely sensitive to the hydrogen [1, 2], neutron radiography (NR) is commonly used to investigate the physical properties, which are dependent on the element in porous media, such as the distribution of water and mineral contents in ceramic, brick and concrete [3–5]. During petroleum exploration, it is proved that the transport property of thermal neutrons (TN) can be used to identify the fluid types bearing in rock and soil [3, 6]. Moreover, NR is often taken as an effective non‐destructive method in the research of pore structure. However, a few research studies have been reported to investigate the anisotropy of pore connectivity with TN in the last decade. Besides, many researchers believe that porosity computed with neutron data cannot be affected by the anisotropy. Therefore, a few research studies have been done to investigate the anisotropic transport of TN. Moreover, it is very difficult to study the anisotropy in the traditional rock experiment, because the intrinsic pore structure cannot be measured directly.

With the application of micron X‐ray CT in digital rock physics (DRP) [7], the advantages of DRP in obtaining pore structure with a high resolution have made it an effective method to research properties of rock core. Based on DPR, static elastic and electrical properties can be simulated by the numerical methods, such as finite element and finite difference methods [8–10]. The lattice Boltzmann method had been used to simulate the transport of fluid in pore space for revealing the effects of micro‐factors on permeability [11–13]. Besides, the research on the acoustic wave and the nuclear magnetic resonance can be efficiently conducted by the numerical simulation methods on the basis of digital rock model [14, 15]. However, since only matrix and pore [7, 16] can be segmented and taken as the effective contents, it is impossible to set up the distribution of chemical multi‐components in conventional models. Therefore, little attention has been paid to the numerical research on the TN properties. Actually, the reaction between TN and material contents filling in matrix and pore will play a significant role in the study of thermal neutron transport. Accordingly, a new method has been proposed to construct a digital rock model including not only the three‐dimensional pore structure, but also the material components. Thus, the transport of neutrons in the new models, considering the reactions between them and the components, can be simulated by the Monte Carlo method (MCM) for revealing the relationship between TN properties and the anisotropy of pore structure. Based on these results, the effects of the anisotropy on the evaluation of porosity can be further analysed through the neutron transport testing along different directions of the new models.

We use two categories of digital core samples obtained by micron X‐ray CT [7, 17]. One group of rock image volumes is constructed from the carbonate (pure limestone) sample with the density of 2.63 g·cm⁻³. The other group is constructed from the pure sandstone with the density of 2.42 g·cm⁻³. The physical resolution of all digital core samples is 2.2046 μm per pixel, and the utilized facility is ultra XRM‐L200. After cropping, we can get fourteen 400 × 400 × 400 carbonate cubes and a 400 × 400 × 400 sandstone cube for further property simulations. The composition of chemical elements in a real rock is so complicated that only the main elements can be considered in our research. Therefore, CaCO₃ and H₂O are set to represent the chemical components of matrix and pore‐filling fluid in carbonate, respectively (pure limestone), whereas SiO₂ and H₂O are set to be the chemical components of matrix and pore‐filling fluid in sandstone, respectively. In both cases, these chemical contents are evenly distributed in their own space respectively.

For transport of thermal neutrons, the macroscopic reaction cross‐section (MRCS) of polyatomic molecule can be written as

$\text{Σ}=\frac{\rho N_A}{A}{\displaystyle \sum }_{i = 1}^n{l}_i{\sigma }_i$ E1

where Σ is the macroscopic reaction cross‐section, cm⁻¹; ρ is the density of molecule, g·cm⁻³; N_A is the Avogadro constant, 6.02 × 10²³; σ_i is the reaction cross‐section of the ith atom for thermal neutron in the molecule, barn (10⁻²⁴ cm²); A is the molecular weight; n is the category number of atoms in the molecule and l_i is the number of the ith atom in the molecule.

Table 1 shows the reaction cross‐section of ${}_1{}^1\text{H}$ , ${}_6{}^{12}\text{C}$ , ${}_8{}^{16}\text{O}$ , ${}_{14}{}^{28}\text{S}\text{i}$ and ${}_{20}{}^{40}\text{C}\text{a}$ , where the first column is the abundance of primary elements, and the last column is the element abundance of their isotopes. Obviously, the abundance of the primary elements is the maximum and is commonly higher than 92%. Thus, it is reasonable for us to consider only these primary elements in our research. The density for the components CaCO₃, SiO₂ and H₂O can be set to the value of 2.71, 2.65 and 1.00 g·cm⁻³, respectively, when the pressure is one standard atmosphere. MRCS of CaCO₃, SiO₂ and H₂O will be: ${\Sigma }_{{\text{CaCO}}_3}$ = 0.27, ${\Sigma }_{{\text{SiO}}_2}$ = 0.20, ${\Sigma }_{{\text{H}}_2\text{O}}$ = 1.52 cm⁻¹. Due to containing more hydrogen, MRCS of pore‐filling fluid (H₂O) is much larger than that of matrix (CaCO₃ and SiO₂). It is clear that porosity will have a significant influence on the transmission of incident neutrons [18], if the contents and the composition of rock model can keep unchanged. The transmittance will decrease with the increase of porosity, which can accordingly cause the change of the intensity of detected thermal neutron.

Element	Cross‐section
	(n, total)	(n, elastic)	(n, nonelastic)	(n, gamma)	Abundance (%)
${}_1{}^1\text{H}$	2.076834e + 1	2.043633e + 1	/	3.320126e − 1	99.99
${}_6{}^{12}\text{C}$	4.742600e + 0	4.739240e + 0	3.360000e − 3	3.360000e − 3	98.93
${}_8{}^{16}\text{O}$	3.852000e + 0	3.851810e + 0	1.900000e − 4	1.900000e − 4	99.76
${}_{14}{}^{28}\text{S}\text{i}$	0.000000	0.000000	0.000000	0.000000	92.22
${}_{20}{}^{40}\text{C}\text{a}$	2.399990e − 3	0.000000	2.399990e − 3	0.000000	96.94

Table 1.

Reaction cross‐section of H, C, O, Si and Ca for thermal neutron and the corresponding natural abundance of isotopes.

Figure 1 shows the two examples of the carbonate and the sandstone core sample, respectively. From Figure 1(b) and (d), it can be observed that the pore space is homogeneous in the sandstone core in YZ, XZ and XY planes, whereas the pore space is heterogeneous in a carbonate core in the three orthogonal planes. The cracks can observed in XZ and XY planes of the carbonate core (Figure 1(b)), but not in the YZ plane, which further demonstrates that its pore space is strongly anisotropic.

Figure 1.

3D reconstructed pore structure of digital core based on micron X‐ray CT. The total size of the dataset is 400 × 400 × 400 voxels with the edge length being 2.2046 μm. (a) The pore space of carbonate core, where the blue part is the pore structure; (b) the surface of the carbonate; (c) the pore distribution of sandstone, where the blue part is the pore structure; (d) the surface of the sandstone.

Neutrons can be classified into different types according to the energy level [19, 20], such as cold neutron, thermal and epithermal neutron. The energy of thermal neutrons involved in our research is 0.0253 eV. In each simulation, these thermal neutrons will irradiate the constructed new models from a vertical direction of one surface plane. On another side of the model, we use array detectors to only receive the neutrons with the energy from 0.0250 to 0.0256 eV in order to record those unreacted thermal neutrons. Many parallel thermal neutron beams will randomly emit from the source (the green part in Figure 2) with the energy of 0.0253 eV for each neutron. The total amount of emitted neutrons is about 40 million moving from the source towards the model. For the neutrons irradiating the core, some of them will react with the core model mostly via elastic scattering, which may change their energy and moving direction during these reactions. Meanwhile, those unreacted thermal neutrons can directly pass through the sample model and reach to the array detectors; therefore, their distribution can be recorded for further analysis and imaging. On the basis of above fact, we use the neutron transmission imaging method to obtain the distribution of difference of counts in YZ, XZ and XY planes.

Figure 2.

Model of thermal neutron irradiating core samples. The cube in the middle is the core sample; the plane on the left is the 50 × 50 array detector; the neutron source (square surface source) is at the right end of the lines. The detector, the model and the source are supposed to be in a straight line, and all of them are set under normal atmosphere. Moreover, they are of the same size in sections perpendicular to the line.

2. Methods

2.1. Reconstruction of core with X‐ray CT

The method to reconstruct the digital rock model is based on micron meter X‐ray CT imaging. Wilhelm Röntgen firstly found the X‐ray in 1895, and X‐ray CT was invented by Hounsfield and Cormack in 1979. Generally, the practical application of X‐ray CT was in medical fields for scanning bone. Rapidly, it had spread into many fields [21], such as materials research and rock physics, as a non‐destructive technique [17, 22]. For those monochromatic X‐rays passing through the medium, some of them will react with the medium, which leads to the energy attenuation. Different components are of different attenuation coefficients, the transport process can be represented as below

$I=I_0{e}^{{-\int }^{\text{}}\mu {\displaystyle \left(s\right)}{d}_s}$ E2

where I₀ is the initial intensity of X‐ray; I is the detected intensity after passing through the medium and μ(s) is the local attenuation coefficient along the transport paths. The components of a substance can be identified by measuring the attenuation coefficient of X‐ray in experiments. For a cylindrical rock sample with diameter of a few millimetres or less, the resolution of X‐ray CT can reach a few micron, or even better. The reconstruction of the new model used in our research can be divided into three steps including image acquisition, image processing and three‐dimensional (3D) reconstruction. With micron X‐ray CT instrument, lots of raw radiographs (as shown in Figure 3(a)) can be obtained through scanning a real rock sample, such as carbonate or sandstone core. Then, the radiographs will experience a series of processing before reconstruction, including noise reduction, smoothing, cropping and segmentation [7]. After that, the segmented binary data (400 × 400) with matrix and pore voxels (Figure 3(b)) can be achieved for each radiograph. Based on the binary data, the three‐dimensional digital rock model can be reconstructed, which can allow us to consider both the pore structures, but also the material components, during the numerical simulation. The chemical component of pore‐filling fluid is only water (H₂O), whereas the component of matrix is only CaCO₃ for carbonate and SiO₂ for sandstone, respectively. Therefore, these different components will not anymore be considered as a uniform mixture according to the ratio of elements like that in the traditional methods. Instead, they are unevenly distributed in their own space in a form of completely separated.

Figure 3.

(a) Cross‐sectional view of carbonate image with voxel of 2.2046 μm; (b) binary image with the blue areas representing the pore space and grey representing the matrix.

The scale of the model should be optimized before simulation in order to reduce the consumption of time and computer memory. Therefore, we have utilized the representative elementary volume (REV) [23] method to implement the optimization by investigating the variety of porosity with the scale. The scale of REV can be determined when the porosity stops changing with the increase of model size. REV scale will give us the minimum elementary volume used to represent the whole sample in simulation. Actually, the specific REV is 400 × 400 × 400 voxels in our research.

2.2. Lattice gas automata

The lattice gas automata (LGA) are derived from the cellular automata by improving the evolution and introducing the collision of particles. The complete discrete LGA model, named as the Hardy, Pomeau and de Pazzis (HPP) model, was introduced in 1973 to simulate the two‐dimensional (2D) fluids flow through dividing the computational field into square lattices. For all LGA model, the discrete particles with unit mass and velocity are supposed to compose the fluid, they will move in the lattice following some local rules. At each step, the discrete particles will move forward one lattice unit to their neighbouring node along their headed directions. Then, they will undergo collisions following the rules that must conserve mass and momentum. Therefore, the local continuity and momentum conservation conditions can be satisfied for incompressible fluids in LGA. Due to the square shape of discrete lattice, the lattice tensor will be anisotropic in the HPP model, which may cause disagreement with the description of the Navier‐Stokes Equation (NSE), and lead to some unreasonable simulating results consequently. A hexagon lattice was introduced in the Frisch, Hasslacher and Pomeau (FHP) model to simulate the viscous flow, which has six neighbouring nodes at each site for the particles to move after collision. With such lattices, the FHP model has successfully solved the asymmetric problem that is inherent in the HPP model.

In LGA, the evolution of distribution function of particle density can be used to describe the collision and streaming of particles, as follows:

$f_i{\displaystyle \left(x+e_i, t+1\right)}=f_i{\displaystyle \left(x, t\right)}+\text{Ω}{\displaystyle \left(f_i{\displaystyle \left(x, t\right)}\right)}$ E3

where the f_i(x, t) represents the distribution function at node x, at time t, and moving along the direction i with the velocity of e_i. The Ω represents the local rules utilized to control the redistribution of particles.

Based on the simple rules, the NSE for flow of an incompressible fluid could be retrieved at macro scale from the evolution of particle. Later, the improvement of FHP by introducing the rest particles makes the LGA method more reasonable for simulation of fluid flow.

2.3. Lattice Boltzmann

As a powerful and flexible approach for mesoscopic numerical simulation, especially in mono‐phase and multi‐phase fluid flow in porous medium, the lattice Boltzmann method (LBM) has been developed to a series of discrete space models including the two‐dimensional model with multi‐speed, since 1988. In the two‐dimensional model (D2Q9), the square lattice with nine velocity vectors is used in discrete space of the computational domain and each node has eight neighbouring lattice nodes. At each time step, a particle at a node described by the distribution function of particle density moves to a neighbouring node, meanwhile, some particles move to this node together and undergo a collision at next time step to form a new distribution of particles.

The time evolution of LBM equation can be written as

$f_i{\displaystyle \left(x+e_i, t+1\right)}-f_i{\displaystyle \left(x, t\right)}=-\frac{1}{\tau }{\displaystyle \left(f_i{\displaystyle \left(x, t\right)}-f_i^{eq}{\displaystyle \left(x, t\right)}\right)},\text{ i}=0,1,2,3,\dots ,8$ E4

where the f_i(x, t) is the distribution function of particle density at lattice site x, time t, moving along the direction i with the velocity of e_i. The τ is relaxation time and f_i^eq(x, t) is the equilibrium distribution of particles.

It was proved that the Navier‐Stokes equation (NSE) can be derived from the time evolution of LBM in a mesoscopic scale using multi‐scale and Champman‐Enskog expansion technology with definition of density and momentum as the function of particle distribution and unit velocity vector, respectively.

$\rho {\displaystyle \left(x,t\right)}={\displaystyle \sum }_{i=0}^b{f}_i{\displaystyle \left(x,t\right)}, \rho u{\displaystyle \left(x,t\right)}={\displaystyle \sum }_{i=0}^b{f}_i{\displaystyle \left(x,t\right)}{e}_i$ E5

2.4. Monte Carlo method

The simulation for researching the transport and scatter of thermal neutrons in the model is conducted by the Monte Carlo method [24, 25] in our research. As a random or statistical method, the Monte Carlo method (MCM) is widely used in many physics fields, especially in nuclear physics. MCM has been proved that it can realistically model the experiment of nuclear physics [24, 25]. In this research, we focus on revealing the effects of elastic scattering on the transport of thermal neutron in the digital rock model. The neutrons moving in the medium will undergo several collisions to lose energy. Between every two collisions, they are supposed to transport along a straight line without loss of energy. Usually, it is composed of five key steps to determine the states during collisions, including setting up the initial state, locating the position, identifying atomic nucleus, determining collision type and measuring the direction and energy of neutrons after collision.

As described in Section 2.1, the constructed digital model, based on X‐ray CT, consists of 400 × 400 × 400 pixel cubes containing the matrix and the pore space inside. Therefore, the model can give the projected image of 400 × 400 pixels in the YZ plane. In the image, the value of each pixel is porosity calculated with the 400 cubes along the X direction of the pixel. Thus, the projected image of porosity in the YZ plane can be used for a further analysis in our research. In the same way, the projected image of porosity in XZ and XY planes can be obtained, respectively. The same colour bar has been applied in the figures of both porosity image and the counts difference to eliminate the effects of the data magnitude, as shown in Figure 4. Accordingly, we can conveniently calculate the similarity between them in these figures by comparing each other. The pictures can be taken as two‐dimensional matrices through converting them into grayscale, thus, the similarity between them can be computed by Eq. (6), as follows:

$r=\frac{{{\displaystyle \sum }}_m{{\displaystyle \sum }}_n{\displaystyle \left(A_{mn}-\overline{A}\right)}{\displaystyle \left(B_{mn}-\overline{B}\right)}}{\sqrt{{\displaystyle \left({{\displaystyle \sum }}_m{{\displaystyle \sum }}_n{{\displaystyle \left(A_{mn}-\overline{A}\right)}}^2\right)}{\displaystyle \left({{\displaystyle \sum }}_m{{\displaystyle \sum }}_n{{\displaystyle \left(B_{mn}-\overline{B}\right)}}^2\right)}}}$ E6

Figure 4.

Thermal neutron imaging of the core. The similarity of (a) and (b) is 0.8385; the similarity of (c) and (d) is 0.7322; the similarity of (e) and (f) is 0.8063. The resolution of (b), (d) and (f) is 17.6368 μm.

where r represents the calculated similarity, A_mn is the element value of matrix A, B_mn is the element value of matrix B, $\overline{A}$ is the average value of all elements in A and $\overline{B}$ is the average value of elements in B.

2.5. Neutron imaging method

The elements of a material can be efficiently identified by the neutron transmission imaging method, which is generally taken as a non‐destructive detection method in the past decades. Hydrogen is the most effective neutron moderator, which mainly distributes in the pore‐filling fluid (for example, water) as constituents of the rock. Therefore, the elastic scattering between thermal neutrons and hydrogen will mostly occur in the pore space. Thus, we can research the pore structure through detecting the neutron. The attenuation of monochromatic (single wavelength) neutron beam transporting in medium can be expressed as Eq. (7) [26, 27]

$I=I_0{e}^{-\mu \tau }$ E7

where I is the intensity of neutron beam transporting in medium, I₀ is the initial intensity of incident neutron beam, μ is the attenuation coefficient, cm⁻¹ and τ is the sample thickness.

For the method, the key steps in simulation will be implemented following this process: pure matrix model without pore space, such as pure carbonate or sandstone, was first irradiated by TN beam to obtain the reference background. Then, the reconstructed digital rock (as shown in Figure 1) was put into the same environment for conducting the numerical experiment of simulating the transport of TN beam along directions of X‐, Y‐ and Z‐, respectively. Consequently, we can obtain the count difference through subtracting the simulated data of the reconstructed model, from that of the pure matrix model. Finally, we will image the distribution of count difference to further analyse the pore structure.

3. Results

3.1. TN imaging of the model

Firstly, the carbonate model is selected for the numerical experiment. TN beam will irradiate the model along directions of X, Y and Z, respectively. The porosity images in YZ, XZ and XY planes have been, respectively, plotted in Figure 4(a), (c) and (e). Correspondingly, the count difference has been plotted in Figure 4(b), (d) and (f). As shown in Figure 4, the similarity of Figure 4(a) vs. (b), Figure 4(c) vs. (d) and Figure 4(e) vs. (f) can be calculated by Eq. (6), which can quantitatively verify that the new method is effective. The obvious count difference in Figure 4(b) can be observed in the middle and right part of the YZ plane, which is consistent with Figure 4(a). It is clear that the distribution of pore space in Figure 4(d) is uniform in the XZ plane of the model. A crack has been measured in Figure 4(f) in the XY plane. Therefore, it is intuitively observed that the distribution of porosity will obviously vary in different planes of the same model, according to Figure 4(b), (d) and (f). This observation can further demonstrate that the carbonate model has strong heterogeneity and anisotropy.

Ten cores were selected from the fourteen 400 × 400 × 400 carbonate model. For each core, we have conducted the numerical experiments and data processing along different directions, the corresponding count difference (CD) has been obtained, as shown in Figure 4(b), (d) and (f). On basis of this fact, the total count difference (TCD) and porosity can be calculated for each core. Relationship between TCD and total porosity has been constructed through linear fitting. The following equations are the relationship in different direction:

$\text{Along X}-\text{direction}:{\varphi }_t=1.62546589556\times {10}^{\mathrm{-7}}x-5.02730504255319\times {10}^{\mathrm{-4}}$ E8

$\text{Along Y}-\text{direction}:{\varphi }_t=1.62519532496\times {10}^{\mathrm{-7}}x-5.73821908133962\times {10}^{\mathrm{-4}}$ E9

$\text{Along Z}-\text{direction}:{\varphi }_t=1.62543339256\times {10}^{\mathrm{-7}}x-4.24629774524067\times {10}^{\mathrm{-4}}$ E10

where x represents the TCD obtained by the simulation along a given direction and ϕ_t is the real total porosity of rock core. We have obtained three relationships between ϕ_t and x along X, Y and Z directions, respectively, they are Eqs. (8)–(10). All the three fittings have the squared correlation coefficient R² = 0.99, which can demonstrate the reliability of them. It is clear that all the three equations are in good agreement with the simulated data. To verify the validity of the equations, they are utilized to predict porosity of the rest core samples. The TCD obtained along the X‐direction will be substituted into Eqs. (8)–(10) to compute porosity of each core. The predicted results are shown in Table 2.

Counts difference (x)	Real porosity	Relative error X (%)	Relative error Y (%)	Relative error Z (%)
33,2780	0.053464516	0.23	0.08	0.38
66,1348	0.108629469	−1.50	−1.58	−1.43
309,735	0.049879859	−0.07	−0.23	0.08
426,181	0.068628984	0.21	0.09	0.32

Table 2.

Prediction of the porosity.

In Table 2, the first column is TCD obtained along the X‐direction, and the second column is the real porosity of each core. The relative errors of predicted results with the three equations have been show in the next three columns, respectively. From the last three columns, it is observed that the relative errors of the predicted porosity will have small difference varying with the directions of neutron transport, when porosity is calculated by the three equations with the same simulated data. Therefore, it is reasonable to deduce that the linear relationship between TCD and porosity will be affected by the anisotropy, which is a contrast to the conventional conclusion.

The neutron transmission image in Figure 4 intuitively shows the distribution of pores, cracks and fractures in rock. On the basis of this fact, the heterogeneity and anisotropy of the model can be analysed with these images. In our research, the scale of array detectors will gradually vary from 4 × 4 to 400 × 400 for revealing the effects of resolution. It is the edge length of each detector in the array that can decide the resolution. Therefore, a variable range of resolution from 2.2046 to 220.46 μm can let us conveniently investigate the effects of the resolution on evaluation of anisotropy. Moreover, we repeat the numerical experiments on the carbonate core as shown in Figure 1 along the X‐direction and record the results with the detector array of different resolution. Thus, a series of images with different resolutions can be obtained through processing the results, as shown in Figure 5. By comparing Figure 4(a) with Figure 5, it is clear that a higher resolution does not always lead to a better image, which means that there is a critical resolution for the neutron imaging. Besides, the similarity (correlation coefficient) between the original model and the detected image has been computed by Eq. (6) under different resolutions. Figure 6(a) shows the calculated results, where the upper points and line represent the carbonate, the lower points and line represent the sandstone, the abscissa is the resolution and the ordinate is the similarity. According to Figure 6(a), it is observed that the similarity will be gradually improved with the increase of the resolution, then reach at the maximum resolution and rapidly decrease after that. Accordingly, we have developed a new equation to relate the resolution to the similarity as follows:

$y=A{x}^a{e}^{b{x}^a}-B$ E11

Figure 5.

Detection results of carbonate core with different resolution. (a) is a image with resolution 110.23 μm; (b) resolution 17.6368 μm; (c) resolution 2.2046 μm.

Figure 6.

The resolution and the number of thermal neutrons emitted in each simulation.

where x represents the resolution of array detector and y represents the similarity and A, B, a, b are constant parameters with related to the porous medium.

Specifically, the relation between the similarity and the resolution in Figure 6(a) can be written as

$\text{Carbonate core}:y=14.67x^{-0.2554}{e}^{-1.951x^{-0.2554}}-1.923, R^2=0.971$ E12

$\text{Sandstone core}:y=0.1285x^{0.9649}{e}^{-0.0877x^{0.9649}}+0.3025, R^2=0.984$ E13

Obviously, Eq. (11) is suitable for both carbonate and sandstone, although they have different material compositions and pore distributions.

For the carbonate core, the simulation has been repeated eight times by modifying the amount of the emitted TN. Figure 6(b) shows the effect of the resolution on the similarity under different intensities of emitted neutron. From Figure 6(b), it is observed that the effect of the resolution on the similarity is similar with that in Figure 6(a), although the amount of emitted TN is different in each repeat. In order to make it clear, the rectangular zone is enlarged and is shown in the lower right of Figure 6(b). We can see that the peak of similarity will vary with the number of emitted TN. Generally, the peak will increase with increasing the TN number. However, the trend will stop when the TN number is more than 40 million. Moreover, it is observed that the peak only appears in a range of the resolution from 10 to 20 μm. Currently, the maximum physical resolution for array detector is about 15 μm [20]. However, based on our results under the condition of 40 million TN, the optimum resolution of array detectors should be set to 17.6368 μm for carbonate in our research.

3.2. Simulations for electrical transport properties

In order to investigate the electrical transport properties, LGA has been applied to simulate the current flow in digital rock samples virtually saturated with oil and water, for revealing the relationship between resistivity index and water (oil) saturation (I‐Sw) at a pore scale. Based on the digital core, we use LGA to research the variation of the resistivity with the saturation for investigating the origin of the non‐Archie I‐Sw relation. The simulated results indicate that I‐Sw relation is a non‐linear function in a log‐log plot, which is a contrast to the Archie’s formula. With the decreasing of water saturation and porosity, the I‐Sw relation will gradually turn towards the abscissa (water saturation axis), which means a non‐Archie phenomenon appears in a low saturation zone. In this case, Archie’s formula is valid only when the water saturation is in a high zone. Therefore, a new equation has been developed to describe the non‐Archie I‐Sw relation based on this study, which can be used to improve the accuracy of saturation calculation in reservoir evaluation. The calculated results with the new equation have shown that it may fit the data of laboratory measurements very well.

In Figure 7, the horizontal axis is the water saturation. The vertical axis is the resistivity index. The range of the water saturation applied in the numerical simulation is set to span from 0.1 to 1 to further investigate the widely changing of the I‐Sw relation. Figure 7 shows I‐Sw relation obtained by the virtual experiments of LGA on the digital rock models, with decreasing water saturation to below 0.3. In this figure, the non‐linear I‐Sw relation can be obviously observed, especially when the water saturation is less than 0.3. In a log‐log plot, I‐Sw relation will gradually drop off towards the abscissa with decreasing Sw. Therefore, it is reasonable to deduce that the Archie exponent n is a function of Sw, instead of a constant in Archie’s theory.

Figure 7.

The non‐linear relation of I‐Sw for water saturation from 0.1 to 1.

3.3. Simulation of fluid flow by LBM

The numerical experiments involved in this study are modelling the flow of Newtonian fluid in a three‐dimensional channel with rigid boundaries at the top and bottom. For a stationary flow driven by a constant horizontal pressure gradient, if the velocity is sufficiently low to neglect the inertia terms, the NS equation can be simplified to Darcy’s law, i.e.

$Q=K\nabla P/\mu$ E14

where Q is the volume flow rate in unit of m/s, K is the permeability of the medium in unit of m² and μ is the dynamic viscosity of the fluid.

Single‐phase flow in pore space has been simulated by LBM with a digital rock of 400 × 400 × 400 voxels. The distribution of fluid velocity is plotted in Figure 8. In order to clearly show the distribution inside the rock, only slices had been chosen to draw in the figure.

Figure 8.

Velocity distribution profile of 3D digital rock model.

Based on the simulations, the absolute permeability can be calculated with Eq. (3). After comparing the results with experimental data, it is proved the feasibility and validity of LBM in the research of simulating fluid transport properties at the pore scale. Moreover, the phenomenon of oil‐water two‐phase flow can be simulated by introducing the Shan‐Chen potential model.

4. Discussion

On the basis of the digital core, LGA was utilized to simulate the variation of resistivity with the saturation for revealing the origin of the non‐Archie I‐Sw relation. The simulated results further verify that I‐Sw relation is a non‐linear function in a log‐log plot, which is a contrast to Archie’s formula. Furthermore, we have proved that Archie’s formula is valid only when the water saturation is in a high zone, and the exponent n is a function of Sw, instead of a constant in Archie’s theory. Moreover, the flow of single‐phase Newtonian fluid in a three‐dimensional digital rock had been implemented with LBM for calculating the absolute permeability. The good agreement between the calculated results and the experimental data has further proved the validity of LBM in research of fluid flow in pore space.

A new method has been developed to reconstruct a digital rock model based on the data from scanning a real rock with micron X‐ray CT, which can simultaneously take into account both the pore structure and the material components. MCM is used to simulate the transport of TN in medium at a pore scale. After analysing the results, we have further investigated the effect of anisotropy on the porosity calculation. According to the results, it is proved that the method of TN imaging can intuitively show the distribution of pore, cracks and fractures inside a rock. Moreover, we have found that the linear relationship between TCD and total porosity can be slightly affected by the anisotropy. Therefore, porosity may have a little difference if the data of different directions are utilized, when we use the neutron method to evaluate porosity of porous media. In the new model, the distribution of chemical elements of different components is supposed to be uniform in their space, such as in matrix or in pore space. Hydrogen, as the most effective neutron moderator, mainly distributes in pore‐filling fluids, such as water, oil and natural gas. Thus, the content and distribution of these substances in rock can be quantitatively investigated with the neutron method. The results of our research may offer an evident to support the anisotropic analysis of pore structure with the neutron method. Furthermore, the distribution of substances containing hydrogen‐4 is of great interest in many research fields, which is mainly detected by the neutron radiography at present. Accordingly, our new method may be a significant reference for them due to its flexible ability in handling complex pore structure and components distribution. Besides, we have found that a higher resolution of array detector cannot always lead to a better neutron image, when we research the pore connectivity and distribution with TN transmission. A critical resolution for the optimum image can be determined with our method, under the condition of a suitable TN number emitted in each simulation.

Acknowledgments

The authors thank financial supports from the National Natural Science Foundation of China (Grant Nos. 41374143 and 41074103).