Shale oil plays a growing role in the global energy supply chain as conventional oil reserves decline. The existence of numerous pore–throat structures in shale reservoirs causes a continuous alteration in the flow path’s width during the flow of shale oil (Wang et al. 2023; Nelson 2009; Su et al. 2022; Chauhan et al. 2021). Consequently, it is essential to investigate pore–throat structures and their impact on shale oil flow (Tan et al. 2022).

The primary reason for the existence of pore–throat structures in shale is its high heterogeneity (Cai et al. 2018). Shale consists of various inorganic mineral components, such as quartz, calcite, dolomite, and organic matter (Gupta et al. 2017, 2018; Xu et al. 2022; Hu et al. 2023; Sun et al. 2023a, b). The different shapes of each component particle result in interparticle pores having pore–throat structures. For example, Naraghi et al. (2018) found by SEM that feldspar, calcite, and clay in shale are elliptical, while quartz is pentagonal. Moreover, pore size is multiscale in shale reservoirs (Taghavinejad et al. 2021; Chen et al. 2021; Sun et al. 2022). For organic matrix, the pore size ranges from 10 to 500nm, while for inorganic matrix, it ranges from 10 nm to 100 mm (Naragh and Javadpour 2015; Clarkson et al. 2013). Furthermore, micro-fractures and fractures produced by hydraulic fracturing are much larger (Li et al. 2022; Xu et al. 2019). Consequently, the flow path size often changes abruptly when fluids flow between OM and iOM or from pores to fractures (Lee et al. 2016).

In order to facilitate theoretical research, many researchers used simplified regular shapes to represent pore–throats, which can be divided into two categories: (1) Pore throats that abruptly change in size or flow direction. (2) Pore throats that gradually change in size or flow direction. The size of type (1) can be characterized by piecewise function (Gravelle et al. 2013; Li et al. 2019; Chao et al. 2019), as shown in structure IV in Fig. 1. Type (2) can be represented by a variety of mathematical functions and is widely used for mathematical modeling, as shown in structures I, II and III in Fig. 1. Cai et al. (2019) and Müller-Huber et al. (2016) used exponential function (I in Fig. 1) to characterize pore–throat structure and calculated permeability of porous media. Sinha et al. (2013) and Fyhn et al. (2021) used trigonometric functions (II in Fig. 1) to study the immiscible two-phase flow in single capillary and porous media, respectively. Besides, in microfluidic experiments (Xu et al. 2017) and pore network models (Suh and Yun 2018), channels composed of circular particles (III in Fig. 1) are often used to represent pore–throat structures.

Pore–throat structures of shale oil flow and its mathematical approximation

Experiments, numerical simulations and mathematical model calculations show that the pore–throat structures have many complex characteristics compared with the constant-radius pore. Li et al. (2019) found that flow pattern transformations in viscoelasticity fluids are greatly influenced by the pore–throat structures. The molecular simulation results of Lee et al. (2016), Zhang et al.(2020) and Sun et al. (2022) found that hydrocarbon fluids will encounter an energy barrier in pore–throat structures. In addition, due to the bending of the streamline, an additional hydrodynamic resistance called the entrance effect will generate when flowing through the pore–throat structures (Gravelle, et al. 2013). Sampson (1891) first derived the pressure drop caused by the entrance effect of liquid transport through a circular hole in an infinite film. Then Dagan et al. (1982) extended the result to convergent flow into a cylindrical pore. Based on FE simulations, Gravelle et al. (2013) modified Sampson’s model to allow it to calculate the additional hydrodynamic resistance of the right part of structure IV in Fig. 1. Wang et al. (2019a, b) then applied it to structure II in Fig. 1 to calculate the entrance effect. While studies have investigated the influence of hydrodynamic resistance on shale oil flow in a single pore–throat structure, the differences in hydrodynamic resistance between structures I, II, III and IV (Fig. 1) have not been extensively studied yet, and the effect of additional resistance generated by different pore–throat structures on permeability has not been fully understood.

In this work, FE models of four pore–throat structures were constructed to calculate the total hydrodynamic resistance in Sect. 2.1. Then in Sect. 2.2, based on Gravelle’s model, the resistance caused by inner effect and entrance effect were derived and fitted by FE results. In Sect. 2.3, we constructed the permeability model coupling pore–throat structures on a REV scale. Finally, the effects of pore–throat structures on permeability with different pore sizes, OM content and distribution were discussed.

2.1 Finite element model of pore–throat structures

Based on previous research, four pore–throat structures were used in this work, as shown in Fig. 2. Structure I can be expressed as (Cai et al. 2019):

where, \(r_{{\text{t}}}\) is the pore–throat radius, nm; \(\alpha = \frac{{\ln \left( {r_{\max } /r_{\min } } \right)}}{{L_{{\text{t}}} }}\); \(L_{{\text{t}}}\) is the pore–throat length, nm.

Four pore–throat structures constructed in this work. The structures have same r_min, r_max and L_t

Structure II is (Wang et al. 2019a, b):

where \(a = \frac{{r_{\max } - r_{\min } }}{2}\), \(b = \frac{1}{{L_{{\text{t}}} }}\), \(c = \frac{{r_{\max } + r_{\min } }}{2}\).

Structure III can be written as:

where \(R = \frac{{r_{\max } - r_{\min } }}{2} + \frac{{L_{{\text{t}}}^{2} }}{{2\left( {r_{\max } - r_{\min } } \right)}}\) is the circle radius; \(\beta = \frac{{r_{\max } + r_{\min } }}{2} + \frac{{L_{{\text{t}}}^{2} }}{{2\left( {r_{\max } - r_{\min } } \right)}}\). To avoid \(r_{\max } > r_{\min } + R\), structure III should meet:

The part with varying radius of structure IV can be expressed as:

It should be noted that the entrance effect only occurs at the corner (Gravelle et al. 2013). It is therefore necessary to have a part with a constant radius \(r_{\min }\).

To calculate the total hydrodynamic resistance of the four structures, 2D FE models (COMSOL 3.5) were established, as shown in Fig. 3. Table 1 lists fixed parameters used in FE models. In addition, the Navier slip condition was used to characterize widespread velocity slip in shale nanopores (Wang et al. 2019a, b).

The velocity distribution of four pore–throat structures obtained by FE simulations (COMSOL 3.5). Blue represents the minimum velocity, while red represents the maximum velocity

By analogy with Ohm’s law, hydrodynamic resistance is defined as:

where, R_t is the total hydrodynamic resistance, mPa s/nm²; Q is the flow flux, nm²/s; ΔP is the pressure difference, Pa. between the inlet and outlet for structure I, II, and III, while for structure IV is the pressure difference between x = L_p and outlet.

2.2 Mathematical model of pore–throat structures

Based on the N–S equation and considering the Navier slip condition, the flow flux of a 2D pore with a constant radius is expressed as:

where, Q_slip-p is the flow flux considering velocity slip, nm²/s; η is the oil viscosity, mPa s; L_p is the pore length, nm; l_s is the slip length, nm; r_p is the pore radius, nm. Then, the hydrodynamic resistance can be written as:

where the R_slip-p is defined as inner hydrodynamic resistance. For the four pore–throat structures with different radius, the resistance can be expressed as

For the hydrodynamic resistance caused by the entrance effect, Gravelle et al. (2013) gave an empirical formula in 3D structure IV:

where C_r is a numerical prefactor; θ is the angle between the pore wall’s tangent and horizontal direction. Here, we tried to use Eq. (10) to describe the R_ent of structure IV:

Because the θ of structure I, II and III change with x, the R_ent should be expressed as:

The Maclaurin expansion of sin(dθ) is:

The higher order infinitesimal term of dθ is ignored, and Eq. (13) can be simplified as follows:

where \(\theta\) can be calculated by:

Combing with Eqs. (1–5) and (14–15), the R_ent of structure I–IV can be obtained.

Structure I:

Structure II:

Structure III:

Structure IV:

2.3 Apparent permeability model of shale oil

The main characteristics of shale different from traditional reservoirs are: (1) Containing OM, (2) Rich in nanopores, and (3) Different pore size distributions (PSDs) in OM and iOM (Wood 2021). Therefore, to study how the content and distribution of OM and the PSD affect the influence of pore–throat structures on permeability, we use the model established in previous work (Xu et al. 2021). Table 2 lists the basic parameters used, and below are the steps to follow:

1. (1)

   Using Quartet Structure Generation Set (QSGS) algorithm (Wang et al. 2007) to construct the 2D spatial distribution model of iOM and OM. As shown in Fig. 4 a, the OM (red part) are randomly distributed in iOM (blue part).

   

   Spatial distribution model of OM and iOM. a Red represents OM and blue represents iOM b Constructing pore–throat structures between adjacent grids, red represents OM grids and blue represents iOM grids

2. (2)

   Assigning pore size to each grid. In this work, PSDs of OM and iOM were characterized using logarithmic normal distributions, i.e. (Zhang et al. 2019):

   $$f\left( {r_{{\text{p}}} } \right) = \frac{\Delta p}{{r_{{\text{p}}} \sqrt {2\uppi } \sigma }}\exp \left[ { - \frac{1}{2}\left( {\frac{{\ln r_{{\text{p}}} - \mu }}{\sigma }} \right)^{2} } \right]$$

   (20)

   where σ is the logarithmic standard deviation; μ is the logarithmic mean. Based on Eq. (20), the pore size of each grid was set by Monte Carlo sampling method.

3. (3)

   Connecting adjacent grids with pore–throat structures (Fig. 5b).

   

   Total hydrodynamic resistance of pore–throat structures with different slip length. The L_t and r_min are a 20 nm and b 2 nm. The r_min and r_max are c 2 nm and d 25 nm

4. (4)

   Calculating the permeability of each grid. Combining Eqs. (8), (9) and (14), a grid's total hydrodynamic resistance is given by:

   $$R_{{\text{t}}} = N\left( {R_{{{\text{slip}} - {\text{p}}}} + R_{{{\text{slip}} - {\text{t}}}} + R_{{{\text{ent}}}} } \right)$$

   (21)

   where, R_t is the total hydrodynamic resistance, mPa s/nm²; N is the pore number in a grid. The \(R_{{\text{t}}}\) can be written as follows according to Darcy's law:

   $$R_{{\text{t}}} = \frac{\Delta p}{Q} = \frac{\eta L}{{k_{{{\text{AP}}}} H}}$$

   (22)

   where k_AP is the apparent permeability, nD; H is the cross-sectional length. Combing Eqs. (21) and (22), the k_AP can be obtained:

   $$k_{{{\text{AP}}}} = \frac{\phi \eta L}{{r_{{\text{p}}} \left( {R_{{{\text{slip}} - {\text{p}}}} + R_{{{\text{slip}} - {\text{t}}}} + R_{{{\text{ent}}}} } \right)}}$$

   (23)

   where \(\phi\) is the porosity; r_p is the pore radius, nm.

5. (5)

   Upscaling permeability to the REV scale. In this work, the upscaling module of MATLAB Reservoir Simulation Toolbox (MRST) was used (Lie 2019).

3.1 Hydrodynamic resistance of pore–throat structures

In the work of Gravelle et al. (2013), C_r is related to the slip length, r_max/r_min and r_max/L_t. Therefore, in Fig. 5 we compared the pore–throat structures' total hydrodynamic resistance calculated by FE and mathematical model. In general, the mathematical model fits well with the FE simulation results. Accordingly, we used the mathematical model to analyze the effect of pore–throat structures on shale oil permeability in Sect. 3.2. In order to analyze the hydrodynamic resistance difference among four pore–throat structures caused by inner effect and entrance effect, the results of two slip lengths (0 nm and 1000 nm) are shown in Fig. 5. According to Eq. (9), inner hydrodynamic resistance is negatively related to the slip length. Therefore, the inner effect can be neglected when l_s = 1000 nm, while its contribution to the total hydrodynamic is dominant for l_s = 0 nm. Figure 6 illustrates the proportion of R_ent to R_t (taking structure II as an example). It is found that the contribution of R_ent can exceed 40% even if l_s = 0 nm. Therefore, when calculating the flow flux of pore–throat structures, it is impossible to ignore the entrance effect.

The ratio of hydrodynamic resistance caused by entrance to the total hydrodynamic resistance (Structure II)

It can be seen in Fig. 5 that the total resistance between different pore–throat structures varies significantly. As shown in Fig. 5a, the R_t of structure III is 133% larger than that of structure IV on average, and in Fig. 5c, structure III is 283% larger on average. It follows that accurate simulation of shale oil flow requires accurate characterization of pore–throat structure. In addition, it should be noted that the magnitude relationship between the two kinds of resistance is opposite. For inner hydrodynamic resistance, structure III > I > II > IV, while for R_ent, structure IV > II > I > III. Figure 7a visually compares the average size of four structures (\(\overline{r}_{{\text{t}}}\)). Because the size of structure III < I < IV and inner hydrodynamic resistance is negatively correlated with the slip length, the resistance shows the characteristics of IV > I > III. The size of structure II is smaller than structure IV for [0, L_t/2) while larger for (L_t/2, L_t], which means that the first half contributes more to inner resistance. However, the average size of structures II and IV is the same as shown in Fig. 7b. Therefore, we can safely conclude that even for inner hydrodynamic resistance, replacing the pore–throat structure with a constant-radius pore with equal average size will overestimate the flux flow.

Size comparison of different pore–throat structures. a intuitive comparison of sizes b average size comparison of structure II and IV

According to Eq. (14), R_ent is positively correlated with \(C_{\text{r}} \cdot \int_{0}^{{L_{\text{t}} }} {\text{d}\theta }\), while negatively correlated with r_t. It can be found in Fig. 8a that the \(C_{\text{r}} \cdot \int_{0}^{{L_{\text{t}} }} {\text{d}\theta }\) of structure I > II > III > IV. However, the average pore size when calculating R_ent is structure IV (r_min) < III < I < II. Therefore, results of the R_ent (IV > II > I > III) are more complicated and should be treated carefully when calculating flux flow. In addition, in Fig. 8b we analyze the prefactor C_r of four structures. It can be found that with the increase in slip length, C_r increases up to a plateau value.

Comparison of parameters affecting the entrance effect

3.2 Effect of pore–throat structures on permeability

In this part, pore–throat structures III and IV were used to analyze the effect of pore–throat structures on permeability. In addition, the slip length in shale nanopores is about 1 nm (Xu et al. 2021). Therefore, the C_r of l_s = 1 nm was used. For structure III C_r = 0.65, while for structure IV C_r = 1.11.

3.2.1 Different pore size distribution

According to Fig. 5, the hydrodynamic resistance of pore–throat structures is closely related to r_max/r_min. Therefore, the effect of PSD should be considered. In Fig. 9a, the average pore size of iOM was increased while that of OM was maintained. It can be found that the effect of pore–throat structures on AP steadily increases with the μ_iOM. There are two reasons: (1) increasing μ_iOM increases the r_max/r_min for iOM grids adjacent to OM grids. (2) the inner hydrodynamic resistance of pores with a constant radius decreases. Similarly, Fig. 9b shows that the increasing of σ enhances the effect of pore–throat structures. This is because increasing σ increases the heterogeneity of pore size, which increases the probability of large r_max/r_min. In addition, when exp(μ_iOM) = 300 nm, the AP without pore–throat structures is 59.6% and 40.8% larger than the AP with structures III and IV, and the AP without entrance effect is 17.6% and 18.8% larger than that with pore–throat structures III and IV. It is indicated that the effect of pore–throat structures on AP cannot be neglected on the REV scale.

Effect of pore–throat structures on apparent permeability with different pore size distributions. a Different μ of iOM b Different σ of OM and iOM. AP_noThroat represents permeability without pore throats, AP_noEntrance represents permeability with pore throats but no entrance effect, AP_t represents permeability considering all factors

3.2.2 Different organic matters content

The TOC of shale can vary from less than 1% (Tran and Sinurat 2011) to 50% (Yao et al. 2018). Therefore, the influence of pore–throat structures on AP with different TOC was evaluated. As shown in Fig. 10a, when TOC < 0.4, permeability is increasingly influenced by pore–throat structures with increasing TOC. This is because there are more iOM grids to establish pore–throat structures with OM grids, which increase the r_max/r_min. However, the effect of pore–throat structures decreases when TOC > 0.4. This may be because the OM grides adjacent to the iOM grids is adjacent to the OM grids after the TOC increases. To prove this point, Fig. 10b shows the effect of pore–throat structure with different OM numbers (TOC = 0.5). It can be seen that the pore–throat structures contribute more with increasing OM number. As shown in Fig. 11, more particles of OM mean a larger area where OM can make contact with iOM. The r_max/r_min of pore–throat structures formed by iOM and OM grids is larger than that of OM and OM grids. The effect of pore–throat structures on AP can therefore be strengthened by improving the dispersion of OM.

Effect of pore–throat structures on apparent permeability with different OM properties. A Different TOC; B Different OM number (TOC = 0.5)

OM distribution models with different particle numbers. Note: The TOC = 0.5

The model can be used to quantify the impact of pore throat structure on shale oil transportation capacity. This study, however, has some limitations due to its empirical basis. numerical prefactor C_r needs to be determined by numerical and mathematical model fitting, while calculation requires integration, which can lead to relatively complex applications. This study was conducted as a preliminary attempt to investigate the effects of different pore throat structures on oil transport. Future research will be conducted.

In this work, the total hydrodynamic resistance of four pore–throat structures was calculated by FE model. Then we established a mathematical model to fit the FE results and analyze the effect of pore–throat structures on AP. Results indicate that both the hydrodynamic resistance due to the entrance effect and the total resistance vary greatly among different pore–throat structures. In addition, the effect of pore–throat structures on AP cannot be neglected, otherwise the error can exceed 60%. As shale contains numerous pore–throats, accurate characterization of pore–throat structure and seepage resistance is needed to calculate permeability. As a result of the findings of this study, we can better understand the effect of pore throat structures on shale oil flow and recovery.