1 Introduction

Low-permeability (\(< 50{\text{ mD}}\)) oil and gas reservoirs are being widely considered by the oil and gas exploration and development in recent years (Sun et al. 2017). Compared to conventional reservoirs, pores and throats in low-permeability reservoirs are narrow (Sun et al. 2015; Zhang et al. 2018a; Huang et al. 2018) and have a large specific surface (Xiong et al. 2009). When fluid flows in a low-permeable porous media, the non-Darcy flow will occur, and there will be a pseudo-threshold pressure gradient (Lei et al. 2008; Zeng et al. 2011). Huang (1998) proposed that the existence of boundary fluid can explain the physical nonlinear phenomenon. In a porous media, liquid molecules concentrated on the surface of pores result in a boundary region with higher fluid viscosity and form a stagnant liquid layer that is called adsorption boundary layer (ABL) (Fig. 1), which is one of the main reasons for the occurrence of nonlinear flow (Huang 1998; Song et al. 2010; Huang et al. 2013; Yang et al. 2020c). The mechanism of pseudo-threshold pressure gradient can be interpreted by the molecular interaction theory and ABL theory (Xiong et al. 2009; Li et al. 2015b; Shen et al. 2019). In porous media, fluid viscosity is influenced by the ABL (Song et al. 2016). The physicochemical properties of crude oil affect the viscosity and thickness of this layer (Li et al. 2011; Tian et al. 2015; Guo et al. 2015). Li et al. (2018) built a capillary permeability model and studied the influence of ABL on tight reservoirs. Yang et al. (2018) studied the ABL in the plate space using nuclear magnetic resonance. Song et al. (2019a) proposed a new single-channel flow model for describing nonlinear flow in low-permeability reservoirs based on ABL theory. Most of previous studies (Li et al. 2015a; Tian et al. 2016; Chen et al. 2018) investigated the effect of ABL on fluid flow in a single pipe, whereas real porous media are more complicated. Therefore, it is of great significance to study the influence of ABL by considering more complex and realistic porous media. The inner structure of rocks, in micro- and nanoscales, can be revealed through high-resolution imaging techniques such as focused ion beam scanning electron microcopy (FIB-SEM) (Desbois et al. 2008, 2011; Li et al. 2015c) and computed tomography (CT) (Coles et al. 1998; Schembre and Kovscek 2003; Shabaninejad et al. 2018; Yang et al. 2020a, b).

Fig. 1
figure 1

Diagrammatic sketch of adsorption boundary layer (ABL)

Full size image

In this study, we investigate the effects of thickness and viscosity of ABL on fluid flow in porous media by constructing digital rock models using the CT scanning method. Pore geometries are extracted from real rock samples, and Delaunay triangulation algorithm is utilized to generate finite element grids, which is explained in Sect. 2 (Lee and Schachter 1980; Shewchuk 2002). The fluid flow simulation considering a stagnant ABL for digital rock is performed using Navier–Stokes (N–S) equations. Because it is difficult to derive formulations to calculate the ABL thickness theoretically, and even more challenging to define the existence of ABL directly in the numerical simulations, an indirect approach is employed to study the absorption boundary layer. In this work, we calculate an average viscosity of fluid by considering viscosity and thickness of ABL. The flow behavior of digital rock models under the effect of ABL is discussed for different rock samples.

2 Reconstruction of geometric model

According to different physical characteristics, four core samples are selected to study the effect of ABL. Berea sandstone, which has a high degree of homogeneity, is used as standard stone to compare the calculation results. The sandstone samples 1, 2 and 3 are from low-permeability formations in Shengli oil field. Samples 1–3 have been selected because they have different pore structures, and thus, we could gain a better understating of the effect of ABL on fluid flow with respect to different pore configurations.

The CT scanning technique provides a nondestructive way to image internal structures of a rock sample (Shabaninejad et al. 2018; Song et al. 2019b). In this study, the instrument used for X-ray CT scanning is Zeiss Versa MicroXCT-400. Its minimum spatial resolution is 0.35 μm, the view pixel is 20482, and the focus size is 5 μm. The resolution of samples 1–3 from Shengli oil field is 3.7594 μm, and the Berea sandstone downloaded from (Dong and Blunt 2009) is 5.345 μm.

The CT scanning yields a three-dimensional (3D) gray image for each sample. Prior to segmentation, for image enhancement, we apply a non-local means filter to enhance the signal-to-noise ratio (Buades et al. 2005). Thereafter, we use the watershed segmentation method to the filtered image to convert grayscale image into two phases (pore and matrix) so that the resulted pore structures can be used for geometric modeling (Saarinen 1994). The pore structures of digital core samples are shown in Fig. 2. The size of samples is 150 × 150 × 150 pixel3.

Fig. 2
figure 2

3D pore structures of digital core samples

Full size image

The main parameters of these core samples are presented in Table 1. The porosities of samples are calculated based on reconstructed digital rock with Avizo software, and permeabilities are estimated through flow simulation using Lattice Boltzmann method (Qian et al. 1992; Ren et al. 2015; Zhang et al. 2019). By adopting the pore network model, we obtain the average radius (Yang et al. 2015, 2019; An et al. 2016).

Table 1 Parameters of sandstone core samples
Full size table

For low-permeability sandstones, pore structure characteristics play an important role in flow mechanism. Therefore, the pore structure of the digital rock sample is characterized by extracting a pore network model (Dong and Blunt 2009; Wang et al. 2020). From Fig. 3a, the pore radius of samples 1, 2 and 3 have a range mainly within 20 μm. From Fig. 3b, we note that the throat radius of Berea sandstone is generally larger than those of samples 1, 2 and 3. In addition, the throat radius of samples 1, 2 and 3 range less than 10 μm. The peak values of pore and throat radius distribution curve of Berea sandstone are larger than that of three low-permeability samples. This indicates that low-permeability sandstone samples are mainly composed of pores and throats with smaller size.

Fig. 3
figure 3

Distribution of a pore and b throat radius of the samples obtained from pore network model

Full size image

We discretize the geometrically complex flow domain of a binary model using unstructured mesh and simulate fluid flow using the finite element method (FEM) (Madadi and Saadatfar 2017). A high mesh density is assigned to the target areas to capture the effects of the adsorption layer at narrower pore throats. Figure 2 shows pore structures of Berea sandstone sample and samples 1–3 with isolated pores. Since isolated pores have no contribution to the fluid flow, before meshing the geometry, they are removed to reduce the computational load. We use Delaunay triangulation approach, a typical surface rendering method, to mesh the geometry (Shewchuk 2002; Fabri and Pion 2009). And, we put a 3D image array of 0 and 1 without isolated pores into Iso2Mesh generator to form finite element grids (Fang and Boas 2009). In the process of modeling, we optimize the quality of the target model by adjusting the model parameters, which is easy to maintain the authenticity of model to a maximum extent. The results of pore structure meshing are shown in Table 2 and Fig. 4. The resulted grid file is imported into COMSOL Multiphysics software for flow simulation in the z direction.

Table 2 Elements of meshing results about pore geometry
Full size table
Fig. 4
figure 4

Meshing results of pore structure without isolated pores. Isolated pores are removed to reduce computational load. Connected pore structure after meshing is used to simulate the flow

Full size image

3 Numerical simulation of flow considering the ABL effects

3.1 Flow equations and model establishment

Fluid flow in unconventional low-permeable rocks is influenced by confinement effect, and the applicability of flow equations should be discussed (Bahukudumbi and Beskok 2003; Barber and Emerson 2006; Zhang et al. 2012; Yao et al. 2013). Nanoscale liquid flow behavior was investigated by molecular dynamic simulations and nanofluidic experiments, which demonstrated that the continuity assumption is valid at nanoscale when the size of flow channel exceeds several nanometers (Sparreboom et al. 2010; Wang et al. 2016; Wu et al. 2017; Zhang et al. 2018b). In this paper, we only consider liquid phase, and the pore and throat sizes are in the micrometer scale (as shown in Fig. 3); therefore, the continuity assumption is reasonable and N–S equations are applicable.

We assume that the fluid flow in porous media is continuous and incompressible. When the gravity is neglected, the N–S equations can be expressed as follows:

$$\frac{{{\text{d}}\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}}{u} }}{{{\text{d}}t}} = - \frac{1}{\rho }\nabla p + \frac{\mu }{\rho }\nabla^{2} \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}}{u}$$
(1)

where \(\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}}{u}\) is flow velocity and \(\rho\) is density. According to Huang’s conclusion (Huang 1998; Huang et al. 2013), we expect that the fluids in both low-permeable samples (samples 1, 2 and 3) and Berea sandstone sample follow the laminar flow. Therefore, we adopt laminar flow pattern to simulate the fluid flow at a pressure difference of 7.5 Pa and a fluid density of 998.3 kg/m3.

3.2 Simulation results of models without ABL

In order to study the effect of ABL in different samples, flow simulation is performed before and after considering the boundary layer. The velocity fields of models without existence of ABL are shown in Fig. 5, from which we observe the fluid flows mainly in a relatively short path between the inlets and outlets. In addition, we note that the pores with smaller radius have faster velocities. In general, the velocity becomes minimum near the throat wall whereas maximum at the throat center. Therefore, the velocity increases gradually from the surface to the center of the throat.

Fig. 5
figure 5

Global slice diagrams of velocity field in different samples without ABL. The legend represents the range of velocity on global slice. The red arrows indicate flow direction. The fluid flows in the z direction

Full size image

3.3 Influence of ABL viscosity and thickness

We define dimensionless thickness and viscosity as:

$$h_{\text{r}} = \frac{h}{{r_{0} }},\;\mu_{\text{r}} = \frac{{\mu_{1} }}{{\mu_{2} }}$$
(2)

where hr, h, r0 are the relative thickness of ABL, ABL thickness, and throat radius, respectively; μr, μ1, μ2 are the relative viscosity of ABL, ABL viscosity and bulk fluid viscosity, respectively. The viscosity of fluid in a throat can be calculated by Huang (1998):

$$\mu = A_{\text{r}} \mu_{1} + \left( {1 - A_{\text{r}} } \right)\mu_{2} .$$
(3)

The ratio of the cross-sectional area of adsorption boundary region to the total cross-sectional area of throat is (Huang 1998; Huang et al. 2013):

$$A_{\text{r}} = 2\frac{h}{{r_{0} }} - \left( {\frac{h}{{r_{0} }}} \right)^{2} .$$
(4)

Given different values of h and μr, μ can be calculated from Eq. (3). The ABL thickness and viscosity remain constant during our flow simulations. Here, we assume that μ2 is 1 mPa·s; therefore, the value of μr is equal to μ1. In order to quantify the influence of adsorption layer thickness and viscosity, we assume that the range of h is 0–3 μm and μr is 2–10, respectively.

By assuming the thickness of ABL is constant, we investigate the effect of ABL viscosity on the microscopic fluid flow in porous media. The velocity fields of our models with consideration of ABL (h is 0.2 μm, μr is 5) are shown in Fig. 6. Taken the influence of ABL into account, the flow velocities in all samples are generally lower compared with those simulation results without ABL.

Fig. 6
figure 6

Global slice diagrams of velocity field in different samples with ABL. (h is 0.2 μm, μr is 5.) The legend represents the range of velocity on global slice. The red arrows indicate flow direction. The fluid flows in the z direction

Full size image

In this paper, when ABL is considered for a given thickness and viscosity, the outlet flow rate is denoted by Qsim, whereas the flow rate without ABL is expressed as Q0. And \({{Q_{\text{sim}} } \mathord{\left/ {\vphantom {{Q_{\text{sim}} } {Q_{0} }}} \right. \kern-0pt} {Q_{0} }}\) is used to quantitatively analyze the effect of ABL thickness and viscosity. For instance, \({{Q_{\text{sim}} } \mathord{\left/ {\vphantom {{Q_{\text{sim}} } {Q_{0} }}} \right. \kern-0pt} {Q_{0} }} \, = \, 1\) means that no ABL is considered.

From Figs. 7, 8, 9 and 10, with an increase in ABL viscosity, the simulated flow rate decreases and the absolute value of the slope increases. When the viscosity of ABL is low, the initial relationship between \({{Q_{\text{sim}} } \mathord{\left/ {\vphantom {{Q_{\text{sim}} } {Q_{0} }}} \right. \kern-0pt} {Q_{0} }}\) and ABL thickness is approximately linear. However, as the viscosity of ABL increases, the data points gradually deviate from the straight line, presenting a nonlinear relationship. As the viscosity of ABL increases, the drop of \({{Q_{\text{sim}} } \mathord{\left/ {\vphantom {{Q_{\text{sim}} } {Q_{0} }}} \right. \kern-0pt} {Q_{0} }}\) slows down at the same ABL thickness, which indicates the influence of ABL viscosity on fluid flow is limited to a certain extent. The relationships between \({{Q_{\text{sim}} } \mathord{\left/ {\vphantom {{Q_{\text{sim}} } {Q_{0} }}} \right. \kern-0pt} {Q_{0} }}\) and ABL thickness of samples 1–3 show more nonlinear behavior than the Berea sandstone sample, which indicates the thickness and viscosity of ABL have a more significant effect on low-permeability porous media, and this is mainly caused by the larger pore radius and better connectivity in the Berea sandstone sample compared with samples 1–3.

Fig. 7
figure 7

Influence of ABL thickness (0–0.2 μm) on fluid flow in the Berea sandstone sample under different ABL viscosities

Full size image
Fig. 8
figure 8

Influence of ABL thickness (0–0.2 μm) on fluid flow in sandstone sample 1 under different ABL viscosities

Full size image
Fig. 9
figure 9

Influence of ABL thickness (0–0.2 μm) on fluid flow in sandstone sample 2 under different ABL viscosities

Full size image
Fig. 10
figure 10

Influence of ABL thickness (0–0.2 μm) on fluid flow in sandstone sample 3 under different ABL viscosities

Full size image

Figure 11 shows the slope values of fitting curves in Figs. 7, 8, 9 and 10 that are plotted as a function of the relative viscosity of ABL. Note that for a given viscosity, the low-permeable samples have higher absolute slope values compared to the Berea sandstone sample, which indicates that the thickness and viscosity of ABL have greater influence on low-permeability rocks. This observation is consistent with the characteristics of the throat radius shown in Fig. 3b in which the throat radius of samples 1–3 is generally lower compared with those of the Berea sandstone sample. Due to the smaller size of throat radius in low-permeable rocks, when the thickness and viscosity of ABL increase slightly, the flow resistance increases significantly, which results in a lower flow rate at the same pressure gradient.

Fig. 11
figure 11

Slope values of fitting curves in Figs. 7, 8, 9 and 10

Full size image

Figures 12, 13, 14 and 15 show the effects of wide-ranged ABL thickness for given ABL viscosities on microscale flow of Berea sandstone and samples 1–3, respectively. We note that for a given ABL viscosity, there is a logarithmic relationship between ABL thickness and \({{Q_{\text{sim}} } \mathord{\left/ {\vphantom {{Q_{\text{sim}} } {Q_{0} }}} \right. \kern-0pt} {Q_{0} }}\). As the thickness increases, the decline trend of \({{Q_{\text{sim}} } \mathord{\left/ {\vphantom {{Q_{\text{sim}} } {Q_{0} }}} \right. \kern-0pt} {Q_{0} }}\) becomes slower. For a given ABL thickness, as the ABL viscosity increases, the downtrend of \({{Q_{\text{sim}} } \mathord{\left/ {\vphantom {{Q_{\text{sim}} } {Q_{0} }}} \right. \kern-0pt} {Q_{0} }}\) slows down gradually. Our simulation results demonstrate that the increase in ABL viscosity and thickness lead to stronger flow resistance and significantly hinder liquid flow in pore throats.

Fig. 12
figure 12

Influence of ABL thickness (0–3 μm) on fluid flow in the Berea sandstone sample under different ABL viscosities

Full size image
Fig. 13
figure 13

Influence of ABL thickness (0–3 μm) on fluid flow in sandstone sample 1 under different ABL viscosities

Full size image
Fig. 14
figure 14

Influence of ABL thickness (0–3 μm) on fluid flow in sandstone sample 2 under different ABL viscosities

Full size image
Fig. 15
figure 15

Influence of ABL thickness (0–3 μm) on fluid flow in sandstone sample 3 under different ABL viscosities

Full size image

3.4 Effects of ABL on fluid flow in different rock samples

In order to investigate the influence of pore space structure on flow behavior, the simulated flow rates under the effect of ABL for Berea sandstone and low-permeable rock samples are presented in Figs. 16 and 17. Firstly, we assume that the absolute thickness of ABL is a constant value of 0.2 μm and calculate flow rates for different rock samples (Fig. 16). Secondly, the relative viscosity of ABL is also assumed to be a constant (\(\mu_{\text{r}} = 2\)); then, the relationships between flow rate and absolute ABL thickness are obtained (Fig. 17).

Fig. 16
figure 16

Influence of relative viscosity of ABL on simulated flow rate for different rock samples

Full size image
Fig. 17
figure 17

Influence of ABL thickness on simulated flow rate for different rock samples

Full size image

As shown in Figs. 16 and 17, when the thickness and viscosity of ABL increase, the flow rates of all samples get smaller. The flow rate of the Berea sandstone sample is about 100 times larger than that of the other three samples and less influenced by ABL thickness and viscosity. In comparison, for samples 1–3 with ultra-low permeability, the flow rate remains within a very small scale and becomes more sensitive to ABL. Therefore, in low- and ultra-low-permeability rocks, the existence of the ABL will cause great damage to the effective development of reservoirs.

4 Conclusions

The microscale model of porous media based on CT and digital rock analysis can reflect the real pore structure of rocks. Therefore, it can be used to simulate the influence of adsorption boundary layer. Using the digital rock constructed by CT scanning, we meshed finite element grids and performed numerical simulation to study the effect of ABL viscosity and thickness on fluid flow for different rock samples.

With respect to our results, when the viscosity and thickness of ABL are small, the relationship between decrease in flow rate and increase in ABL thickness is approximately linear. On the contrary, it tends to be a nonlinear relationship. The effect of ABL thickness and viscosity for low-permeability porous media is more significant than Berea sandstone. It is due to larger pore radius and better connectivity in Berea sandstone. Compared with conventional reservoirs, the existence of ABL in low-permeability reservoirs causes more damage to reservoirs and makes oil development more difficult.