Characterization of Flow Parameters in Shale Nano-Porous Media Using Pore Network Model: A Field Example from Shale Oil Reservoir in Songliao Basin, China

1. Introduction

Shale oil reservoirs have a wide distribution of pore throats (ranging from nanoscale to micrometer scale), with complex pore geometries and pore throat structures [1,2]. The micro pore-throat structure is an important factor affecting the macroscopic reservoir properties and fluid flow in shale oil reservoirs [3]. Therefore, quantitative evaluation and characterization of the micro pore throat structure of shale reservoirs [4] (pore geometry, size distribution, connectivity, etc.) is of great significance for understanding the fluid flow patterns of shale reservoirs. The commonly used methods for characterizing the complex pore-throst structure of porous media are digital cores [5,6,7,8,9,10] and pore network models [11,12,13,14,15,16,17]. Direct flow simulation methods based on digital cores (such as lattice Boltzmann method [18,19,20,21], computational fluid dynamics method [22,23] and Monte Carlo simulation [24,25]) take a long time to calculate, require a large amount of memory, and are powerless for large-scale flow simulation and parameter sensitivity analysis. Considering the extremely small pore-throat and strong micro heterogeneity of shale oil reservoirs, which make it difficult to obtain characterization units or have larger sizes, direct flow simulation methods based on digital cores of shale reservoirs will be very challenging. Based on previous studies, it has been shown that pore network simulation is basically consistent with core scale experiments in obtaining macroscopic physical parameters (capillary pressure curves, relative permeability curves, etc.) [26,27,28,29,30,31,32,33,34]. Moreover, pore network simulation computation time is relatively short, which is very advantageous for large-scale flow calculations. Therefore, the pore network model is selected as the method for pore level flow simulation in shale reservoirs in this paper.

Although porous flow parameters such as flow rate-differential pressure curve, permeability, capillary pressure curve and relative permeability curve can be obtained through core experiments [35,36], micro-nano pore throat development of shale rock core and fluid flow in it is very slow, it will be time-consuming and expensive to obtain macro physical parameters such as phase permeability through experiments, and the flow rate is small, resulting in large measurement errors. At present, direct flow simulation methods based on digital cores are time-consuming and occupy a large amount of memory, while pore network simulation can quickly and accurately predict macroscopic physical parameters [37,38,39,40,41,42,43], making it a powerful method for predicting macroscopic physical parameters of shale reservoirs. The pore network model uses regular geometry to replace the characteristics of complex pore throat structure, and uses form factor to characterize the irregularity of pore throat structure. Based on the percolation theory, the pore level flow characteristic of porous media fluid is studied [44,45,46]. In 1957, Broadbent and Hannersley [47] first proposed the concept of percolation and pointed out its prospects in the application of porous flow. Afterwards, many scholars conducted extensive research on the simulation of pore network flow. In the 1970s, Dullien [48] applied the percolation theory to the pore network model and explored the flow rules of fluids in the network model. Lenormand [49] proposed the invasion percolation model and applied it to the simulation of the displacement process simulation of the displaced fluid in porous media. Oak et al. [50] used the pore network model to study the oil-water two-phase and oil gas water three-phase relative permeability curves of water wet Berea sandstone. In recent years, Blunt, Ø ren, and van Dijke et al. [51,52,53,54] have conducted extensive research on the effects of reservoir wettability on pore level flow, three-phase flow, and prediction of capillary pressure and relative permeability curves using pore network models. Capillary pressure is the main driving force determining saturation changes in quasi-static network models. For the capillary pressure at the inlet of each network model, the equilibrium position of the interface between the fluid and the fluid is determined according to the fluid displacement model. Based on the criteria of each displacement mode, displacement and fluid location in the network model occur gradually in order. In addition to being used for theoretical research, many scholars have proposed network models which can be used to predict macroscopic physical property parameters according to research needs. Vogel [55,56] predicted the soil relative permeability using a network model generated by a series of thin slice techniques that can characterize the complex pore structure of soil. In 2002, Blunt et al. [57] pointed out that by combining pore-scale physical phenomena with geometrically equivalent pore-throat cross-section shapes, a representative network model could be used to accurately predict capillary pressure curves and relative permeability curves. In 2005, Piri and Blunt [58] established a pore network flow model that included all the important characteristics of pore scale immiscible fluid flow (such as water film, oil layer, wetting hysteresis and wettability change). The corresponding model was used to calculate oil-water relative permeability, invasion path and capillary pressure curve in different displacement processes. In the water wetting system, the oil-water relative permeability curve predicted by the model was in good agreement with the experimental results. With the deepening of research, the capillary pressure curves and relative permeability curves during different displacement processes under mixed wetting conditions have gradually been successfully predicted. The quasi static pore network model for the three-phase flow of oil, gas, and water has gradually been developed.

In view of the above problems, this paper proposes a characterization method of flow parameters in shale nano-porous media using pore network model. The structure of this article is as follows. In section 2, the FIB-SEM experiment and the calculation method of flow parameters in porous media are introduced. The section of results and discussion is introduced in Section 3. Finally, the conclusion is drawn in section 4.

2. Methodology

2.1. Establishment of pore network model based on FIB-SEM

The core sample of this experiment is from the Gulong shale oil reservoir in Songliao Basin, China. The experiment used a focused ion beam scanning electron microscope (FIB-SEM) with a resolution of 2nm, Helios 5 CX DualBeam model, to perform high-precision scanning of the shale sample. The electron beam current range is 0.8pA-176nA, and the acceleration voltage range is 200V-30kV. The Ion beam current range is1pA-100nA, and the acceleration voltage range is 500V-30kV.

For the images generated by FIB-SEM, the angle between the ion beam and the electron beam is 52 degrees. Due to the angle between the milled surface and the imaging axis being not 90 degrees, the original image presents geometric artifacts compared to the real profile. Therefore, it is necessary to cut and stretch the image stack to achieve image alignment. In addition, when FIB cutting on the surface of rock samples, the difference in hardness on the surface of the rock sample can cause a difference in cutting speed, resulting in the formation of marks parallel to the cutting direction on the surface of the rock sample. In the FIB-SEM image, it appears as a vertical stripe noise with varying brightness, resembling a curtain, which is called the "curtain effect". Therefore, the preprocessing of FIB-SEM images mainly involves image alignment and denoising filtering. This article uses Avizo software to process FIB-SEM images, providing the FIB Stack Wizard and image filters that can align and preprocess FIB-SEM image stacks. Firstly, the FIB-SEM image stack is aligned and shadow corrected using the FIB stack wizard. Then, the image is denoised and filtered using non-local mean filtering. Then, the striped artifacts caused by the "curtain effect" in the image are eliminated through fast Fourier transform filtering. Finally, the filtered image is edge enhanced through non-sharpening masking to complete image preprocessing. The final FIB-SEM image is shown in Figure 1.

After preprocessing, the signal-to-noise ratio (SNR) of the image can be effectively improved, and then the preprocessed image needs to be segmented, mainly to gain a deeper understanding of the pore structure of the rock and identify the pore phase. Therefore, the Avizo software is used to combine two segmentation methods, interactive threshold segmentation and interactive top hat segmentation, in order to accurately segment the pore structure. Then, the pore network model with true topological relationships is obtained through the pore centered axis method and the maximum sphere method. Figure 2 shows the three-dimensional pore structure of the Gulong shale core after FIB-SEM segmentation and the extracted pore network model.

2.2. Quasi-static single-phase flow equations considering diffusion

The flow rate into and out of a certain pore in the pore network model satisfies the law of mass conservation (Figure 3(a)) [59,60]:

$${\displaystyle \sum _{j=1}^{N_i}{q}_{ij}}=0$$

(1)

Where, N_i is the number of pores connected to pore i; q_ij is flow into or out of pore i; j is the pore number connected to pore i. The flow from one pore into or out of its adjacent pore q_ij is satisfied as follows:

$$q_{ij}=g_{ij}\frac{p_i-p_j}{L_{ij}}$$

(2)

Where, g_ij is conductivity; p_i, p_j are pore pressure of pore i, j; L_ij is the distance between two adjacent pores. The conductivity g_ij can be obtained by harmonic average of the conductivity of adjacent pores and the middle throat (Figure 3(b)):

$$g_{ij}=\frac{L_{ij}}{\frac{L_i}{g_i}+\frac{L_t}{g_t}+\frac{L_j}{g_j}}$$

(3)

In addition to pressure difference, diffusion also contributes to flow rate in shale nano-pores. Therefore, based on Knudsen number, the mass transfer modes of Gulong shale oil at different scales are divided. In addition to differential pressure mass transfer, Knudsen diffusion at small pore size, molecular diffusion at large pore size and transition diffusion at medium pore size are respectively established to characterize diffusion coefficients at different scales. Knudsen number K_n is defined as follows:

$$K_{\mathrm{n}}=\frac{\lambda }{d}=\frac{k_{\mathrm{B}}T}{\sqrt{2}\pi {\delta }^{2}pd}$$

(4)

where, λ is the mean free path; d is the equivalent pore-throat diameter; k_B is the Boltzmann constant, 1.3806488×10^-23 J/K; δ is the collision diameter of molecule, m; T is the temperature, K；p is the total pressure, Pa.

Gas diffusion mechanisms in porous media mainly includes Fick diffusion, transition diffusion and Knudsen diffusion. When K_n < 0.1, the flow regime is primarily continuous and slip flow, mainly occurring in large pores. At this time, gas molecules collide with each other and ignore collisions with the pore wall. Diffusion takes the form of Fick diffusion. When K_n > 10, the flow regime is primarily free molecular flow, mainly occurring in small pores. At this time, only a single molecule can pass through the pore, and molecules collide with the pore wall. Diffusion takes the form of Knudsen diffusion. When 0.1 < K_n < 10, the flow regime is primarily transitional flow. At this time, gas molecules and the wall both experience collisions, and diffusion takes the form of transitional diffusion, as shown in Figure 4.

When calculating gas flow velocity caused by diffusion, one of the most important variables is the diffusion coefficient. When Knudsen diffusion is dominant, the diffusion coefficient is the Knudsen diffusion coefficient D_K, which can be calculated using the following formula:

$$D_K=\frac{d}{3}\sqrt{\frac{8RT}{\pi M}}$$

(5)

When Fick diffusion is dominant, the diffusion coefficient is the Fick diffusion coefficient D_F, which can be calculated using the following formula:

$$D_F=\frac{k_{B}T}{3\sqrt{2}\pi {\delta }^{2}p}\sqrt{\frac{8RT}{\pi M}}$$

(6)

Where, M is molecular molar mass, g/mol; R is the ideal gas constant, 8.314462 J/(mol·K). When transitional flow is dominant, the diffusion coefficient is the transition diffusion coefficient D_T, which is a combination of the Fick diffusion coefficient and the Knudsen diffusion coefficient. The calculation formula is as follows:

$$\frac{1}{D_{\mathrm{T}}}=\frac{1}{D_{\mathrm{K}}}+\frac{1}{D_{\mathrm{F}}}$$

(7)

Therefore, after considering diffusion, the flow equation in a single pore-throat is as follows:

$$Q=Q_{\mathrm{p}}+Q_{\mathrm{d}}=\left(\frac{\pi d^4}{128\mu L}\frac{\rho }{M}+\frac{\pi d^2}{4L}\frac{D_t}{RT}\right)\Delta p$$

(8)

The pressure at each pore node of the network model can be obtained by combining (1), (2), (3) and the single-phase conductivity coefficients of any shape section. Inlet and outlet nodes were used to calculate the inlet and outlet flow, and Darcy's law was used to calculate the absolute permeability K of the network model:

$$k_{\mathrm{abs}}=\frac{Q\mu L_{\mathrm{m}}}{A\Delta p}$$

(9)

Where, Q is the outlet or inlet flow of the network model, m³/s; L_m is length of pore network model, m; A is the cross-sectional area of the network model, m²; △p is the differential pressure at both ends of the network model, Pa.

2.3. Quasi-static two-phase flow equations

The relative permeability of oil and water in the initial displacement process can be calculated at any point in the displacement process. In the two-phase flow process, if there is at least one connected cluster in the network model (the connected cluster in the two-phase flow is defined as: a fluid of a phase in the network model is connected to the inlet and outlet through pores and throats), then the relative permeability of the phase is not zero. On the connected cluster, for each phase of fluid, the flow into and out of a certain pore satisfies the law of conservation of mass [59,60]:

$${\displaystyle \sum _j{q}_{p,ij}}=0$$

(10)

The flow rate of phase p between pore i and its connected pore j is satisfied as follows:

$$q_{p,ij}=\frac{g_{p,ij}}{L_{ij}}\left(P_{p,i}-P_{p,j}\right)$$

(11)

The conduction coefficient g_p,ij between two adjacent pores can be used as the harmonic average of the conduction coefficient between two adjacent pores and the connecting throat:

$$\frac{L_{ij}}{g_{p,ij}}=\frac{L_i^n}{g_{p,i}^n}+\frac{L_{ij}^b}{g_{p,ij}^b}+\frac{L_j^n}{g_{p,j}^n}$$

(12)

Where, g_p,ij is the effective conduction coefficient between pore i and connected pore j; L_ij is the distance between pore i and adjacent pore j; P_p,i, P_p,j is the pore i and p phase fluid pressure in j.

In the initial displacement process, oil is the non-wetting phase and water is the wetting phase. At the throat section with corner, the water film is located at the corner, and the oil phase is located at the center of the pore. In the process of two-phase flow, when two phases of oil and water coexist in the same pore or throat, the oil phase conduction coefficient is the bulk phase conduction coefficient, and the water phase conduction coefficient is the sum of all the existing water film conduction coefficients. The expression is:

$$g_o=g_{o,b}，g_w={\displaystyle \sum _{k=1}^n{g}_{w,c}^k}$$

(13)

By combining the above formulas, the pressure of p-phase fluid at each pore node in the network model can be obtained, and then the effective permeability of p-phase fluid can be obtained by using formula (13):

$$k_p=\frac{{\mu }_p{Q}_p{L}_T}{A_T\Delta p}$$

(14)

The calculation formula of relative permeability of p phase is shown in Equation (14):

$$k_{\mathrm{r}p}=\frac{k_p}{k_{\mathrm{abs}}}$$

(15)

In the initial displacement process, water saturation of the pore network model is:

$$S_{\mathrm{w}}=\frac{{\displaystyle \sum _{i=1}^M\frac{{\displaystyle \sum _{j=1}^n{A}_{c,j}}}{A_i}{V}_i+{\displaystyle \sum _{i=1}^N{V}_i}}}{{\displaystyle \sum _{i=1}^{M+N}{V}_i}}$$

(16)

Where, S_w is the water saturation of the pore network model; V is volume of pore network model, m³; A is the cross-sectional area of the pore network model, m².

Table 1. Parameters used in pore network flow simulation.

Simulation parameter	Unit	Value
Pore number	/	8717
Throat number	/	18494
Oil density	g/cm3	0.7
Water density	g/cm3	1
Oil viscosity	mPa·s	0.5
Water viscosity	mPa·s	1
Oil-water contact angle	°	60
Oil water interfacial tension	mN/m	20
Initial temperature	℃	137
Initial pressure	MPa	37.5

Table 1. Parameters used in pore network flow simulation.

Simulation parameter	Unit	Value
Pore number	/	8717
Throat number	/	18494
Oil density	g/cm3	0.7
Water density	g/cm3	1
Oil viscosity	mPa·s	0.5
Water viscosity	mPa·s	1
Oil-water contact angle	°	60
Oil water interfacial tension	mN/m	20
Initial temperature	℃	137
Initial pressure	MPa	37.5

3. Results and discussion

3.1. Pore-throat parameters characterization

The three-dimensional slices of Gulong shale core were obtained by FIB-SEM, and the three-dimensional digital core was constructed. The total porosity is about 9.6%, but the effective connected porosity is only 4.6%. Figure 5(a) and (b) respectively show the radius distribution of pore and throat. The minimum pore radius and throat radius are respectively 2nm and 1nm, the maximum pore radius and throat radius are respectively 438nm and 334nm, the mean pore radius is 87nm, and the mean throat radius is 41nm. Figure 5(c) and (d) respectively show the coordinate number and shape factor of Gulong shale core, and the mean coordination number is 3.97. Generally speaking, pore connectivity of shale is not poor, and the difficulty of development of shale oil reservoir is due to the extremely small pore throat radius.

3.2. Apparent permeability characterization

Figure 6 shows the flow rate-differential pressure curve calculated with and without diffusion. The permeability predicted by the combined action of differential pressure displacement and diffusion in micro and nano pores is greater than that under pure differential pressure mass transfer condition, and the permeability calculated with diffusion is 0.000124mD. The permeability calculated without considering diffusion is 0.0000531mD.

3.3. Relative permeability curves characterization

Figure 7 shows the distribution of oil and water in the pore network model under different water saturation conditions. Injected water preferentially flows along large pore channels, and only with the increase of water saturation can injected water gradually enter small pores.

The relative permeability curve of oil and water in Gulong shale reservoir was calculated through the simulation of oil-water two-phase pore network. The irreducible water saturation of Gulong shale reservoir was calculated to be about 0.26, and the residual oil saturation was about 0.73. The relative permeability curve is shown in Figure 8.

4. Conclusion

In this paper, a characterization method of flow parameters in shale nano porous media using pore network model is proposed to characterize the key flow parameters of the Gulong shale reservoir in the Songliao Basin, China. The conclusions are as follows:

(1) The minimum pore radius and throat radius of Gulong shale are 2nm and 1nm respectively, the maximum pore radius and throat radius are 438nm and 334nm respectively, the mean pore radius is 87nm, the mean throat radius is 41nm, and the mean coordination number is 3.97.

(2) The permeability predicted by considering the combined effects of medium pressure differential displacement and diffusion in micro nano pores is greater than that under pressure differential mass transfer conditions. The calculated permeability considering diffusion is 0.000124mD, while the calculated permeability without diffusion is 0.0000531mD.

(3) The oil-water relative permeability of Gulong shale oil reservoir is calculated through oil-water two-phase pore network simulation. The irreducible water saturation of Gulong shale oil reservoir is about 0.26, and the residual oil saturation is about 0.73.