Digital-Rock Construction of Shale Oil Reservoir and Microscopic Flow Behavior Characterization

: In shale oil reservoirs, nano-scale pores and micro-scale fractures serve as the primary fluid storage and migration space, while the associated flow mechanism remains vague and is hard to understand. In this research, a three-dimensional (3D) reconstruction of the shale core and micro-pore structure description technique is established; digital core technology for shale reservoirs was developed using X-ray computed tomography (X-CT), scanning electron microscope (SEM) and a focused ion beam scanning electron microscope (FIB-SEM). Microscopic oil–water two-phase flow is mimicked using the lattice Boltzmann method (LBM), a well-acknowledged approach to exploring nanoconfined fluid dynamics. In addition, coupled with digital cores, the flow characteristics of shale reservoirs are characterized. The total porosities of bedding fractures in shale and lamellar shale are 2.042% and 1.085%, respectively. The single-phase oil flow inside bedding fractures follows Darcy’s linear flow principle. This work can deepen the understanding of the microscopic flow characteristics of continental shale reservoirs and provide a reference for similar problems that may be encountered.


Introduction
An apparent difference between shale and conventional reservoirs is that the pores in the shale matrix are mainly nano-pores with very low permeability.There are many types of shale pores with complex structures and obvious multi-scale [1].Because of the micro-nano scale characteristics of shale pores, using traditional methods [2] to study the rock properties generally lacks quantitative description and cannot shed light on the fluid flow behavior and mechanism in the core microstructure.In the early 1990s, Dunsmuir [3] advanced X-ray computed tomography (X-CT) scanning technology and introduced it to the petroleum industry for the first time.Scholars have investigated the flow mechanism in the microstructure porous media in the petroleum engineering industry using digital core technology [4][5][6].Based on the computed tomography (CT) scanning images of the realistic core sample collected from the field [7], digital core technology and a 3D model of the physical structure are established using reverse reconstruction software.Rosenberg et al. [8] introduced X-CT scanning imaging technology into the field of the digital core and established the 3D digital core of Fontainebleau sandstone.Gao [9] verified the accuracy and reliability of digital cores constructed by CT and suggested that a stricter constraint shall be applied to the alteration principle of average throat radius.Based on CT images, Liu [10] used the pore-scale reconstruction model to simulate water and CO2 displacement, revealing the influence of pore characteristics on water displacement efficiency and the two-phase flow mechanisms.
Compared to traditional numerical methods, the lattice Boltzmann method (LBM) [11,12] is an efficient method widely used to simulate fluid percolation characteristics in a digital core.The LBM contains different discrete models for mathematical and physical problems one may face in fluid flow modeling [13,14], and its complete parallelism dramatically improves calculation efficiency.In the 1990s, Martys [15] used the LBM model of two-phase immiscible fluid to estimate the relative permeability of 3D digital cores, and its small model size shall take responsibility for the undesirable reliability of calculation results.Walsh [16] established a new partial boundary regression algorithm, which has mass conservation in heterogeneous media and derived the analytical expression for calculating the permeability of different models.Leclaire [17] performed 3D experiments on water absorption and drainage in digital Berea sandstone, and three primary states of fluid invasion are reproduced [18].Based on the LBM and the multi-component Shan-Chen method, Lautenschlaeger [19] proposed a homogenization method to simulate the fluid-fluid and solid-fluid interactions in the porous media.To solve the problem of inaccurate throat simulation with the help of machine learning, Rabbani [20] combined the advantages of pore network modeling (PNM) and LBM for the first time, and LBM is utilized to calculate throat permeability.Zhao [21] established an improved PNM to simulate single-phase flow in porous media, using LBM to simulate the fluid flow in realistic throats.Suo et al. [22] combined LBM and PNM, PNM modeled the coarse-scale multiphase transport in porous media; meanwhile, LBM solved the multi-component fluid dynamics.
To sum up, the digital core technology based on CT can genuinely reproduce the internal microstructure inside the core, and its visual accuracy advantage is irreplaceable by other methods.LBM has developed into a widely used and efficient numerical simulation method, which will play an increasingly important role in stimulating the flow characteristics of rocks.In this paper, the digital core constructed by the CT scanning method realizes the realistic microscopic pore structure of the shale core, presented in Section 1.After that, in Section 2, scanning electron microscope (SEM) images are processed to obtain the characteristic parameters of the pore structure in the shale matrix.In Section 3, the PNM of different matrix cores was obtained.Subsequently, LBM is used to quantitatively describe fluid flow characteristics in the microscale pore network in Section 4. At last, conclusions are drawn in Section 5.

3D Reconstruction of Digital Shale Core by CT Scanning
The digital core technology has been widely used to advance oilfield development [23,24], including micro-flow mechanism research, oil displacement mechanism research, and reservoir production performance simulation and prediction.Meanwhile, it is realized that CT scanning is a favorable method for reconstructing the digital core, as it can reproduce the micro-pore structure of the actual core with high accuracy.

Digital Core Reconstruction Method
The core target selected in this paper is the typical shale reservoir.There are six cores divided into pure shale and Lamellar shale, and the basic parameters of the above six cores are summarized in Table 1.CT detection and digital core reconstruction using the above core mainly include the following six steps: a. Drill 25 mm diameter core on the full-diameter core with wire cutters.b.Use a CT scanner to conduct 3D scanning of a 25 mm diameter core with a scanning accuracy of 13.6 μm.c.Processing the measurement results using Avizo software, identifying identifiable bedding cracks by threshold segmentation, and calculating the width, spacing, and porosity of bedding cracks.d.Drill a core with a diameter of 2 mm at the center of the core with a diameter of 25 mm.e. Use a CT scanner to scan the core with a diameter of 2 mm in three dimensions and a scanning accuracy of 1 μm.f.Processing and calculating the measurement results by using Avizo software.

Digital Core Reconstruction Method
Figure 1 shows the 3D reconstruction results of 25 mm and 2 mm diameter shale cores; the dark flake and strip areas are shale bedding fractures, and the dark spots are matrix pores; and Table 2 shows the statistical table of bedding fracture characteristic parameters of shale.The porosity of fractures in shale cores with a width greater than 13.6μm is about 1.153%, and the average fracture spacing is about 5.15 mm.The porosity of fractures in lamellar shale cores with a width greater than 13.6 μm is about 0.571%.

Scanning Images of SEM
Taking core fragments, the shale samples polished by argon ion were scanned by scanning electron microscope in the direction parallel to the bedding plane and perpendicular to the bedding plane of the core, and five viewing zones were taken in each direction for scanning.Figure 2 shows the approximate position of the viewing areas in different directions.The characteristic parameters of matrix pore structure can be obtained through image processing and statistics.In this chapter, six cores in the first Section are scanned by SEM.Figures 3-5 shows the scanning results of CZ-Y 47-1.Although the resolution of Figure 3d is 20 μm, its HFW value increases accordingly, so these images can be reasonably used in the process of digital core reconstruction.

Pore Structure Analysis
After processing the SEM images, the parameters such as specific surface area, porethroat radius, tortuosity, shape factor and pore--throat coordination number were obtained.In this Section, CZ-Y47-1, one of the six shale samples is taken as an example, after binarizing the SEM images; as shown in Figure 6, the black pixels are matrix pores, and the matrix surface porosity results in different directions are shown in the title.The plane porosity is the ratio of the pore area to the total image area.Because of the anisotropy of shale, the core porosity is about the average porosity of each plane porosity.Firstly, the processed images are statistically analyzed.According to Equation (1), the plane porosity of three different observation directions can be calculated: 1.8%, 3.6%, and 3.4%, as shown in Figure 6: where Φf is the plane porosity; Ap is the total pore area, μm 2 ; Aall is the total area of view, μm 2 ; Φ is core porosity; Φfv1 is the plane porosity in #1 vertical bedding direction; Φfv2 is the plane porosity in #2 vertical bedding direction; Φ fp is the plane porosity in parallel bedding direction.Then, the average core porosity of the CZ-Y47-1 core is calculated as 2.9% according to Equation (1).According to the above statistical results and Equation ( 2), the pore diameter distribution is the ratio of the number of pores with different diameters to the total number of pores; it can be estimated that the average pore diameter of the shale matrix is close to 0.11 μm: where Dpi is the diameter of a single pore, μm; Api is the area of a single pore, μm 2 ; p D is the average pore diameter, μm; N is the total number of pores.
The throat diameter distribution refers to the ratio of the number of throats with different diameters to the total number: where Dti is the diameter of a single throat, μm; Ati is the area of a single throat, μm 2 ; t D is the average throat diameter, μm; M is the total throat number.The tortuosity is the ratio of the actual length of the flow channel to the apparent length through the flow medium, which indicates the complexity of the fluid flow channel: where τti is the tortuosity of a single throat; Si is the actual length of a single throat between two pores, μm; Li is the linear length of a single throat between two pores, μm; t τ is the average tortuosity.The coordination number refers to the number of connected throats of a single pore, indicating the pore network's complexity: where CNpi is the coordination number of a single pore.The shape factor is the ratio of the cross-sectional area of the pore and throats to the square of the perimeter, which can quantitatively show the shape characteristics of the pore-throat unit: where p F is the average shape factor.By calculation according to Equations ( 3)-( 6), the average throat diameter of the shale matrix approaches 0.034 μm; the above results are shown in Figures 7 and 8.The average tortuosity of the shale matrix is 6.53, and the average coordination number of the shale matrix is 1.77, and using the above formulas and image statistics, the parameters are summarized in Table 3.

Construction of Pore Network Model Based on FIB-SEM Tests
PNM, pore network modeling, can reflect the topological structure of the core pore space and has geometric parameters equivalent to the real core pore space [25,26].Shale cores were scanned by focused ion beam scanning electron microscope (FIB-SEM) technology, and the slice spacing was close to 10 nm.Then, the PNM of different matrix cores is obtained through image processing and reconstruction.

Construction of PNM
The pixel size is 5.42 nm, and the effective measurement range is 9100 × 6476 × 3990 nm.The digital core of CZ-Y47-1 was obtained by cutting the image using Avizo software.Image segmentation of organic matter, pore, and matrix is performed to extract the structure of different components, as shown in Figure 11.Reconstructing the binarized images to obtain the segmented 3D photos, as shown in Figure 12, and pore network modeling is generated, as presented in Figure 13.

Structural Characteristics Analysis of PNM
The classification of reservoir space type depends on the structural characteristics of the core pore network, which can be displayed intuitively and clearly by extracting the digital core pore network model.The microstructure of pores can be evaluated in detail by statistical analysis of structural parameters [27].
According to the 3D digital core of CZ-Y47-1, volume percentage statistics of the pore, organic matter, and core matrix are 1.85%, 4.12%, and 95.03%, respectively.The fractal dimension of the generated PNM is 1.95, and the tortuosity is 8.72.The statistical distribution of core structure characteristics is shown in Figure 14.The PNM parameters are shown in Table 4, and it can be observed that the porosity of micropores is relatively large; the increase of micropores enhances the connectivity between some isolated macropores.

Microscopic Flow Characteristics
Based on the reconstructed digital core, LBM can be used to simulate the displacement process and quantitatively characterize the flow characteristics of the fluid in the microporous network [28].The underlying influence of wettability, interfacial tension, and fluid viscosity on the single-phase oil flow characteristics are investigated, expecting to guide shale oil production prediction and development scheme optimization.

Establishment of the Flow Model
In 1992, Qian [29] developed the DdQm model that laid a foundation for the extensive application of the LBM.The LBM can also be deduced from the continuous Boltzmann Equation: where f is the particle number density distribution function, eis the particle velocity, a is the particle motion acceleration, and ( ) f Ω is the collision operator.The collision operator is a complex nonlinear integral form.A linear collision operator is proposed based on Maxwell equilibrium distribution, and the Boltzmann-BGK Equation is obtained: BGK is a continuous equation, and the discretization of velocity and space obtains the discretized lattice Boltzmann Equation.From the macroscopic statistical theory, it can be considered that all molecules are moving in specified directions, so the rate of particle movement is simplified to several rates along the specified direction {e0, e1,…, eN}, and N is the direction number of the rate.The distribution function f of particle number density is also discretized into {f0, f1, …, fN}, which defines fα = fα (r, e, α, t), α = 0, 1, ..., N. The Boltzmann Equation for velocity dispersion can be obtained as follows: The equation is discretized in time and space, and the fully discretized LB Equation is obtained: The last term on the right side of the equal sign is the external force term, and the above formula is the Lattice-BGK Equation considering external force.The LBM comprises lattice type, evolution equation of particle number density distribution, equilibrium distribution function , and boundary condition.Using the D3Q19 lattice, the evolution equation is: The equilibrium distribution function is: e u e u (12) where, s c is the lattice sound speed, , and α ω is the weight coefficient: Consider the boundary conditions of flow: Based on the above equations, a single-phase flow LBM model is established considering the interfacial tension and liquid-solid wettability.The viscosity of the oil phase can influence the relaxation time: The oil phase adsorption layer on the pore wall affects the oil phase flow in the pore.The oil-solid interfacial tension (γ) and contact angle (θ) jointly determine the adhesive force.The contact angle is only used to describe the strength of fluid-solid force and does not indicate the existence of multiphase fluid in porous media; additionally, with the increase of the contact angle, the wettability of the wall gradually tends to water.

Percolation Characteristics of Crude Oil in Matrix Pores
Based on the established shale matrix digital cores with average pore-throat diameters of 130 nm and 30 nm, the permeability of the digital cores is calculated, and the matrix permeability of the shale core is measured by a matrix permeability meter.The simulation parameters are shown in Table 5.Compared with the two models, the error between the model verification measurement results and the simulation results is only 6.2%, clarifying the accuracy and reliability of the proposed model.In this paper, the length scale  equals 5 × 10 −8 m, the mass scale  equals 1.51 × 10    As illustrated in Figure 16, with the increase in the proportion of hydrophilic minerals, the adequate area oil flow increases, and the pressure decreases in large-size pores.The adsorption force the water-wet surface imposes on oil molecules is less than that with the oil-wet surface; as a result, oil-flowing resistance in water-wet pores is reduced compared to oil-wet pores, leading to relatively high flow efficiency.It can be seen from Figure 17 that the starting pressure gradient of single-phase oil flow in the shale matrix decreases with the increase of wetting angle.As can be seen from Figure 18, with the increase of wetting angle, the area where oil can flow increases, and the pressure in large-size pores decreases.Figure 19 shows that the starting pressure gradient of single-phase oil flow in the shale matrix increases with the increase of interfacial tension.Figure 21 shows that as the viscosity of the oil phase increases, the starting pressure gradient of single-phase oil flow in the shale matrix increases.It can be seen from Figure 22 that when the viscosity of the oil phase increases, the oil-flow area decreases, and the pressure in large-size pores decreases.The increasing oil viscosity makes it more difficult to flow in the pores.According to the single-phase starting pressure gradient curve, the average porethroat diameters are 30 nm and 130 nm, and the starting pressure gradient of flow in the matrix is 0.103-0.285and 0.045-0.110MPa/mm, respectively.

Flow Characteristics of the Oil Phase in Bedding Fractures
Based on this digital core, the percolation characteristics of single-phase oil in the bedding fractures are simulated.As seen from Figure 23, there is no starting pressure gradient in the bedding fracture flow, while the effect of stress is not considered, and the velocity has a linear relationship with the pressure gradient.Figure 24 shows the relationship between the single-phase flow velocity of crude oil and pressure gradient under different influential factors.The results show that the singlephase flow follows Darcy's principle, and there is no threshold pressure gradient.The more significant the proportion of hydrophilic minerals is, the greater the wetting angle of the oil phase is, the smaller the interfacial tension is, and the higher the flow velocity of the oil phase is.In addition, the higher the viscosity of the oil phase, the lower the flow rate of the oil phase.

Conclusions
In this paper, the CT scans of cores with diameters of 25 mm and 2 mm were carried out, six 3D digital cores were reconstructed, and the characteristics of shales and lamellar shales were obtained.SEM and FIB-SEM were used to rebuild the digital core, obtain the characteristic parameters of the pore structure of the shale matrix, and reconstruct the 3D pore network model of the shale matrix.Considering the viscosity, interfacial tension, and wettability of fluid, the LBM model of shale matrix/bedding fracture flow is established, and the microscopic flow characteristics of shale oil under different conditions are calculated.The following conclusions were drawn: (1) Total porosities of bedding fractures in shale and lamellar shale are 2.042% and 1.085%, respectively; the average fracture spacing is 0.5 and 0.9 mm, and the average fracture widths are 6.8 and 5.5 μm.
(2) The average porosities of the shale and lamellar shale matrix cores obtained by SEM are 3% and 2.23%, the average pore diameters are 113 and 154 nm, the average throat diameters are 32.6 and 33.3 nm, the average tortuosities are 5.65 and 5.06, and the average shape factors are 0.035 and 0.054, respectively.(3) The average pore throat diameter is 30 nm, and the starting pressure gradients in the matrix are 0.103~0.285and 0.045~0.110MPa/mm, respectively.The mechanism of single-phase flow in shale bedding fractures is the same as Darcy's linear flow.(4) The accuracy of the digital core model constructed by the three modeling methods needs further improvement.Although it can reflect the structural characteristics of the reservoir, the single structure method and mixed numerical reconstruction method still need to be improved to simulate the pore structure and improve the application results accurately.

Figure 1 .
Figure 1.3D CT reconstruction results of shale cores.

Figure 2 .
Figure 2. Schematic diagram of SEM observation position of the rock sample.

Figure 3 .
Figure 3. SEM image of CZ-Y47-1 in the direction parallel to the bedding plane.

Figure 9 .Figure 10 .
Figure 9. Reconstructed digital core of shale with an average pore-throat diameter of 130 nm.

Figure 11 .
Figure 11.Image segmentation and structure extraction with different compositions.
(a) Reconstruction of pore structure (b) Reconstruction of feldspar minerals (c) Reconstruction of matrix minerals

Figure 12 .
Figure 12.Reconstruction of different pores and minerals.
(a) Pore diameter distribution (b) Coordination number distribution (c) Throat radius distribution (d) Throat length distribution

Figure 15 .
Figure 15.Oil pressure diagram of hydrophilic minerals with different proportions.

Figure 16 .
Figure 16.Simulation results of hydrophilic minerals with different proportions.

Figure 17 .
Figure 17.Fluid pressure diagram under different wetting angles of oil.

Figure 18 .
Figure 18.Simulation results under different wetting angles of the oil phase.

Figure 19 .Figure 20 .
Figure 19.Fluid pressure diagram under different interfacial tension.It can be seen from Figure20that with the increase of interfacial tension, the oil-flow area decreases, and the pressure in large-size pores increases.The adsorption force of oil on the pore wall increases, and the flow characteristics become more complicated.

Figure 21 .
Figure 21.Fluid pressure diagram under different oil phase viscosities.

Figure 22 .
Figure 22.Simulation results under different oil phase viscosities.

Figure 23 .
Figure 23.Simulation results of oil in bedding fractures under different pressure gradient conditions.
(a) Proportion of hydrophilic minerals (b) Oil phase contact angle (c) Interfacial tension (d) Oil phase viscosity

Figure 24 .
Figure 24.Relationship curve between oil phase flow rate and pressure gradient under different conditions.

Table 1 .
Basic parameter of collected shale cores.

Table 2 .
Characteristic parameters of shale bedding fractures.

Table 3 .
Summary of pore and throat structure parameters.

Type No. Sample No. Porosity/% Average Pore Diameter/nm Average Throat Diameter/nm Tortuosity Coordination Number Shape Factor
Figures 9 and 10 present digital cores of shale matrix with average pore-throat diameters of 130 and 30 nm, respectively.They are reconstructed by the random growth method, which provides a model for subsequent micro-flow simulation with less than 5% error.