Stochastic-based liquid apparent permeability model of shale oil reservoir considering geological control

Shale is a complex porous medium composed of organic matter (OM) and inorganic minerals (iOM). Because of its widespread nanopores, using Darcy’s law is challenging. In this work, a two-fluid system model is established to calculate the oil flow rate in a single nanopore. Then, a spatial distribution model of shale components is constructed with a modified quartet structure generation set algorithm. The stochastic apparent permeability (AP) model of shale oil is finally established by combining the two models. The proposed model can consider the effects of various geological controls: the content and grain size distribution of shale components, pore size distribution, pore types and nanoconfined effects (slip length and spatially varying viscosity). The results show that slip length in OM nanopores is far greater than that in iOM. However, when the total organic content is less than 0.3 ~ 0.4, the effect of the OM slip on AP increases first and then decreases with the decrease in mean pore size, resulting in that the flow enhancement in shale is much smaller than that in a single nanopore. The porosity distribution and grain size distribution are also key factors affecting AP. If we ignore the difference of porosity between shale components, the error of permeability estimation is more than 200%. Similarly, the relative error can reach 20% if the effect of grain size distribution is ignored. Our model can help understand oil transport in shale strata and provide parameter characterization for numerical simulation.


Introduction
Advances in hydraulic fracturing and horizontal drilling technologies have sparked a surge of interest in exploring and exploiting shale oil (Honarpour et al. 2012;Hou et al. 2016;Zheng and Sharma 2020). Permeability, the most critical property for evaluating the capability of oil flow in shale strata, is essential for the development of shale oil reservoirs ). However, precisely determining oil permeability is challenging, largely due to the shale's various components, their multi-scale pore structure and the special flow mechanisms defying Darcy's law (Tahmasebi et al. 2015a;Zheng et al. 2021).
Two major components exist in shale samples: organic matter (OM) and inorganic minerals (iOM) (Naraghi et al. 2018). The OM mainly consists of kerogen, which can be further divided into three main types (Chen et al. 2015;Kelemen et al. 2007;Ungerer et al. 2014). The iOM majorly consists of quartz, calcite, dolomite, clay, feldspar and pyrite, etc. Both the types and contents of OM and iOM can vary dramatically between different shale plays and even within a play (Xu et al. 2020). Bakken, the largest shale oil plays in the USA, is mainly type I and II kerogen (Onwumelu et al. 2019). The total organic carbon (TOC) in upper and lower members can reach up to 36%, while TOC in the middle member is no more than 1% (Tran et al. 2011). Minerals in Bakken shales mainly consist of quartz, dolomite and calcite (Fakhry and Hoffman 2018). Using laboratory petrophysical measurements, researchers classified shale rocks in Eagle Ford, Barnett, Woodford and Wolfcamp formations into three types, respectively (Gupta et al. 2017(Gupta et al. , 2018. TOC of three rock types in Eagle Ford and Barnett is 1.5-5.5%, and that of Woodford and Wolfcamp is 0.5-9.0% and 1.0-4.0%, respectively. Calcite, dolomite, quartz and clays are the dominant minerals, which can be up to 80-90%. Another complex characteristic is that each component has its distinct pore topology and morphology, such as porosity, pore size distribution (PSD) and pore connection (Chen et al. 2015;Naraghi et al. 2018). Loucks et al. (2012) systematically classified pores of shales into OM pores, mineral interparticle pores (InterP), and mineral intraparticle pores (IntraP). Experiments, including nitrogen adsorption tests, smallangle and ultra-small-angle neutron scattering (SANS and USANS), high-pressure mercury intrusion and nuclear magnetic resonance (NMR), etc., revealed that the OM pores size (10-500 nm) is usually an order of magnitude smaller than InterP and IntraP size (10 nm-100 μm), causing a bimodal pore size distribution (Clarkson et al. 2013;Naraghi and Javadpour 2015;Rafatian and Capsan 2015). However, the porosity in OM can reach as high as 50% (Curtis et al. 2010), while minerals like calcite and quartz usually have ultra-low porosity (Houben et al. 2014).
Due to the aforementioned complexity manifested in shale, stochastic reconstruction methods have been recognized as one of the most reliable ways to produce variability (Tahmasebi et al. 2016b;Wu et al. 2020a). Naraghi et al. (2018) divided the stochastic methods into three categories: object-based, process-based and statistical methods. Objectbased methods try to model the porous medium as a collection of stochastic objects. The properties, such as proportion, shape, and other morphological data, are determined by prior knowledge (Tahmasebi et al. 2016b). Naraghi and Javadpour (2015) proposed a 2D stochastic apparent permeability (AP) model for shale gas. They modeled the OM with squares randomly dispersed in inorganic matrix for the first time. Followed by Cao et al. (2017), Xu et al. (2019) and Feng et al. (2019) extended the 2D model to 3D, and the OM shape was correspondingly extended to block. However, these models did not distinguish the various minerals in iOM. To overcome the limitation, Naraghi et al. (2018) then constructed a new object-based model. Using mineralogy derived from X-ray diffraction and scanning electron microscopy (SEM), they modeled mineral components with different geometry. By honoring mineral types and shapes, they find that the calculated AP of shale gas is closer to the measured value. Process-based methods try to mimic the formation process of porous media (Tahmasebi et al. 2016a). Using the elementary building block (EBB) model, Chen et al. (2015) reconstructed the structure of shale matrix, in which quartet structure generation set (QSGS) was used to mimic the growth of grain size. In their model, clay, calcite, pyrite and OM are considered as an EBB, respectively, and in each EBB, InterP and IntraP are locally defined. Yao et al. (2018) developed a multi-scale pore network model to estimate the intrinsic permeability (IP) of shale. In their study, geological processes, including sedimentation and diagenesis (compaction, cementation, generation of pores inside organic matter and clay), were simulated. Wu et al. (2020a;2020b) proposed a new dynamic and process-based model to consider the geological evolutions that occurred in shale. In their model, with the increase in thermal maturation of shale, three physical processes, including sedimentation, compaction and cementation, are mimicked by QSGS, seam carving and dilation algorithms, respectively. Also, the model can consider the PSD and components content in shale. Therefore, the process-based model can be used to investigate the effects of various geological controls, including total organic matter content (TOC), minerals content, PSD, etc., on the intrinsic properties (like permeability) of shale. The statistical method uses statistical properties and statistical algorithm to generate stochastic shale models. The statistical properties can be obtained from 2D SEM images and experimental data. Based on cross-correlationbased simulation, Tahmasebi and his co-authors proposed a novel and accurate reconstruction method to develop 3D shale models by directly using 2D images (Tahmasebi et al. 2015a(Tahmasebi et al. , 2015b(Tahmasebi et al. , 2016a(Tahmasebi et al. , 2016b. In terms of OM and PSD, these models are in excellent agreement with the original data. However, the components in iOM were not carefully divided. Later, Li et al. (2018) and Ji et al. (2019aJi et al. ( , 2019b further inserted different minerals into shale models with the energy-dispersive spectroscopy (EDS) images. Although statistical methods result in accurate reconstruction, they are difficult to build a realization at the representative elementary volume (REV) scale due to the lack of images at the REV size (Naraghi et al. 2018). What is more, hours or days that may be taken to generate models is also a limitation (Naraghi et al. 2017).
Although an extensive body of studies have discussed the reconstruction of realistic shale, understanding that the liquid flow, especially oil in their dominant nanopores, is still in its infancy (Falk et al. 2015;Fan et al. 2020;Sun et al. 2019). Classical fluid mechanics, like the Hagen-Poiseuille (H-P) equation, usually assumes that there is no velocity slip between liquid and solid wall. However, this convention is no longer valid when the flow path size decreases to the nanoscale. Experimental results have shown that the liquid flow rate (water, hexane, decane, etc.) through carbon nanotubes (CNTs) can be 4-5 orders of magnitude higher than that predicted by the no-slip H-P equation (Majumder et al. 2005;Mattia and Calabrò 2012). To calculate the slip velocity theoretically, based on molecular-kinetic theory, Blake (1990) quantitatively related the slippage to the contact angle between liquid molecules and solid wall. With the same premises, Mattia and Calabrò (2012) further incorporated the surface diffusion and work of adhesion. What is more, they developed a two-fluid system (TFS) model to calculate the water flow rate in CNTs. The near-wall region was characterized by slip length and reduced viscosity. Slip length is defined as the extrapolated distance from a solid wall where the tangential velocity component disappears (Zhang et al. 2018). Followed by Cui et al. (2017), Zhang et al. (2017) and Sun et al. (2019) extended the model to oil flow in tortuous organic and inorganic nanopores. However, validation of this extension is in need because the structure of oil components, like n-alkanes, is very different from water. Recently, based on the different mobility between the first-layer molecules and those in the bulk, Wu et al. (2019) derived the slip length when n-alkanes flow in nanotubes. However, parameters in their equations, like natural relaxation time, are difficult to determine. In recent years, molecular dynamics simulation (MDS) has become an effective tool to obtain the slip length of hydrocarbon flow in nanopores (Fan et al. 2020;Zhan et al. 2020). Wang et al. (2016a;2016c) used MDS to study the octane flow in OM (simplified as graphene), quartz and calcite nanopores. Due to the difference of wall roughness and oil-wall molecular interactions, they found that slip length in different nanopores varies dramatically. Slip length in a 5.24 nm graphene nanopore can be up to 132.48 nm, while in a calcite nanopore of 3.34 nm, it even becomes a negative one (−0.87 nm). Although the slip length in quartz (~ 1.15 nm for a 5.24 nm nanopore) and calcite nanopores is small, ignoring the slip boundary conditions will lead to large errors. Similarly, Sui et al. (2020) explored the octane flow mechanisms in dolomite nanopores. With the same temperature, pressure and driving force as the quartz nanopores, the slip length for a 5.7 nm dolomite nanopore is only 0.44 nm. Therefore, the large difference of slip length between OM and different iOM nanopores further demonstrates the necessity to consider different shale components when calculating the AP of shale oil.
The rest of the paper is arranged as follows. In "Two-fluid system model of shale nanopores" section, a TFS model to calculate the oil flow rate in different shale nanopores is also developed, but slip length and varying viscosity are obtained based on the MDS data. In "Spatial distribution model of shale components" section, a 2D spatial distribution (SD) model of shale components (OM and iOM minerals) with a modified QSGS algorithm is built, which can incorporate the grain size distribution (GSD) of shale components. Then, in "AP model of REA-scale" section, the TFS model is integrated into SD model using the capillary bundle model, and the AP of shale oil is scaled up to the representative elementary area (REA) scale through this way. In "Results and discussion" section, based on the REA-scale AP model, the effect of nanoconfined effects, PSD, pore type, porosity distribution, and GSD on AP are discussed. Finally, important conclusions are summarized in "Conclusions" section.

Two-fluid system model of shale nanopores
In terms of water flow in CNTs, Mattia and Calabrò (2012) developed a TFS model to account for the viscosity change in an annular region near the wall. Similarly, when oil flows in a single nanopore, MDS has demonstrated that about 2-4 layers of oil molecules near the wall form an ordered molecular shell, exhibiting a significant viscosity change (Wang et al. 2016c). Therefore, considering the similar flow behaviors compared with water, the TFS model is employed to calculate oil flow rate in shale nanopores. This model can incorporate both the viscosity change and the slip length (positive, zero and even negative), as shown in Fig. 1.
Based on the no-slip Hagen-Poiseuille equation in a circular pore, velocity in two regions can be expressed as where v b , v nw and b , nw are the oil velocity and the dynamic viscosity in the bulk region and the near-wall region, respectively; Δp L is the externally applied pressure gradient; and R is the nanopore radius. To guarantee the continuity of velocity and shear stress, as well as the mass conservation, boundary conditions of each region should satisfy where l s is the slip length. By solving Eq. (1)-(2) with the boundary conditions given in Eq. (3)-(6), the velocity of each region can be expressed as Assuming the viscosity in each region is uniform and then integrating Eq. (7)-(8), the oil flow rate in each region can be obtained by The corresponding total volumetric flow rate is given by Assuming an effective viscosity can characterize the effect of viscosity change, i.e., eff = b = w , then Eq. (11) now becomes where eff is the effective viscosity and l s ′ is the slip length without considering viscosity change (Fig. 1). Set l s ′ =0 Eq. (12) will convert to the well-known no slip H-P equation To characterize the real nanopores in shale more accurately, tortuosity is introduced into Eq. (11) where is the tortuosity, which is defined as the ratio of the length of a sinuous track to the length of a sample (Srisutthiyakorn and Mavko 2017).
Fig. 1 TFS model for oil flow of shale nanopore modified from Zhang et al. (2017). R is the radius of nanopore, is the thickness of near-wall region, and ls' is the slip length without considering viscosity change. ls is the slip length considering oil viscosity change in near-wall region. From A to C, the molecular interactions between oil and shale nanopore wall are decreasing, which leads to the increase in slip length from a negative one to a positive one The key to solve Eq. (11) is to determine the values of nw and l s . In this study, MDS data from Wang et al. (2016a;2016c) and Sui et al. (2020) are analyzed to obtain these two parameters in graphene, quartz, calcite and dolomite nanopores, respectively. Octane is used in their simulations to represent oil. In terms of light shale oil, it is a good simplification (Wang et al. 2015). However, in their works, an effective viscosity model (Eq. (12)) was used to obtain the slip length, which means that they cannot be used in our work directly. Therefore, we re-fitted the slip length in different nanopores by piecewise fitting the MDS data. With the increase in pore size, at most about four layers of octane are affected by the strong oil-wall interactions, and the thickness of each layer is about 0.45 nm (Wang et al. 2016c). Thus, when the radius of nanopore is greater than 1.8 nm, the in Eq. (7-11) is set to 1.8 nm; otherwise, = R . The fitting process is provided in Supplementary Information (SI), and the results are listed in Table 1 (Fig. 2).

Spatial distribution model of shale components
In this section, a modified QSGS algorithm is used to construct the SD model of shale components. QSGS (Wang et al. 2007) has been demonstrated a powerful tool to mimic the forming processes of porous media (Chen et al. 2015), in which the volume fraction, average particle size and particle shape of each phase can be taken into account. However, size distributions of components in real shale are not uniform. For example, Naraghi and Javadpour (2015) found that the size of OM in Eagle Ford conforms to the normal distribution. Therefore, we modify the QSGS algorithm with an objected-based method to consider the size distribution of shale components. The flowchart to construct the SD model using modified QSGS in a grid system is shown in Fig. 3, and the steps are described as follows: (1) Obtaining data: Before initiation, one should determine the types, area fraction and size distribution of shale components, as well as the REA size of shale. These can be analyzed by SEM-EDS images and X-ray powder diffraction (XRD) (Houben et al. 2014;Loucks et al. 2012). In this work, quartz, calcite and dolomite are selected to represent the iOM. The basic input properties of these four components are summa-     Table 2 and Fig. 4. It should be noted that for the convenience of sensitivity analysis, we assume that the grain sizes of all shale components conform to the same normal distribution. The REA size used here, 150 μm × 150 μm, is a typical value of fine-grained materials (Houben et al. 2014). Grid size is set as 0.5 μm × 0.5 μm (Naraghi and Javadpour 2015), which is smaller than most grain size shown in Fig. 4. (2) Classifying: Classical QSGS algorithm usually regards the pore as a nongrowing phase, while OM and each mineral as a growing phase. To consider the size distribution of shale components, we regard every grain of OM and minerals as a growing phase. The size of each growing phase is constrained by the GSD obtained in step (1). To consider the PSD and facilitate the use of TFS model, instead of regarding pores as nongrowing phases, we use capillary bundle model to generate them in the grid system in "AP model of REA-scale" section. In other words, the modified QSGS algorithm is only used to generate shale components, but not shale pores (Fig. 5).
(3) Seeding and growth: Randomly select a grain size from the size distribution curve of component I and locate the core of the specific grain in a random grid G i,j . Then grow the core to its neighboring grids in eight directions based on the growth probability D i (Fig. 2), until the grain size in grid system reaches the selected value. (4) Repeating: Repeat step (3) until the area fraction of component I reaches the pre-obtained value in step (1) and continue with step (3) for the remaining component. (5) Identifying the InterP grids: After the fraction of every shale component reaches its target value, we regard the interface grids between different mineral grains as InterP grids with another given probability, which will be used to generate InterP in "AP model of REA-scale" section.

AP model of REA scale
Based on the parameters listed in Table 2, we construct a SD model of shale components as shown in Fig. 6. Then, we will generate IntraP and InterP and calculate AP in each grid using TFS model and capillary bundle model. Finally, the AP of shale oil is scaled up to REA scale. The flowchart is shown in Fig. 3, and steps are described as follows: Obtaining data: PSD and porosity can be obtained by SEM images and various experiments cited in "Introduction" section. Detail parameters can be found in Table 2 and Fig. 6. The PSD of OM and iOM is derived from an Eagle Ford Shale sample, which conforms to a logarithmic normal distribution, which can be expressed as where r is the pore radius and and are the mean and standard deviation of natural logarithm of radius.
It should be noted that the porosity given in Table 2 is the area fraction of pores occupying the entire SD model. Therefore, the porosity in each grid should be converted properly. For example, the number of InterP and IntraP grids of quartz in Fig. 6a is 23260 and 3740, respectively, so the porosity in each IntraP grid is 3%×90000/ (23260+3740×2)=8.78% and the porosity in each InterP grid is 17.56%.
Generating pores: Traverse the grid system; in a grid G i,j , the component type (OM, iOM minerals) and pore type (OM pores, InterP and IntraP) have been determined previously. Therefore, we can select a pore size from the PSD of the corresponding component, and then generate pore in G i,j .
Calculating AP in each grid: Based on capillary bundle model and Eq. (14), the volumetric flow rate in each grid can be expressed as where n i,j is the number of pore in grid Gi,j. According to Darcy's law, the volumetric flow rate can also be given by . 6 a SD model of shale components. There are 90,000 grids in the grid system, including 9000 OM grids, 27,000 quartz grids (3740 InterP grids), 27,000 calcite grids (3889 InterP grids), 27,000 dolomite grids (4126 InterP grids). b AP distribution in the grid system. c Velocity distribution in the grid system. The AP of REA scale is 5.86μD ▸ where k AP i,j is the apparent permeability of grid Gi,j; A g is the cross-sectional area of a grid. Combining Eq. (16) and Eq. (19), AP in a grid can be expressed as Upscaling AP to REA scale: In this work, the upscaling module in MATLAB Reservoir Simulation Toolbox (MRST) is used to upscale the AP (Lie 2019). Because it is a single-phase flow in 2D, the upscaling can be completed in a few seconds. The results are shown in Fig. 6.

Model validation
To validate the accuracy of our model, we compare the calculated AP with Zhang et al. (2019). In their study, a generalized lattice Boltzmann model (GLBM) is proposed to calculate the AP of shale oil. Parameters used in their simulation are listed in Table 3. It should be noted that only two components are divided in their work, and slip length is calculated based on the effective viscosity model. Results are compared in Fig. 7. It can be seen that our model slightly overestimates the AP compared with the results of GLBM, but the average relative error is only 3.1%. Therefore, our model is reliable to evaluate the AP of shale oil.

Results and discussion
In this section, based on the AP model, the effects of different geological controls on AP are analyzed, including nanoconfined effects (slip length and spatially varying viscosity) (Wu et al. 2017), PSD, pore type, porosity and GSD. Unless otherwise specified, the model uses the parameters listed in Table 2.

Effects of pore size distribution and nanoconfined effects
The slip length listed in Table 1 is closely related to the pore size. Hence, it is attractive to analyze the relationship between AP and PSD. Figure 8 shows the sensitivity of AP to PSD of iOM and OM, respectively. Both the increase in mean pore radius (e μ ) in iOM and OM will increase the AP, but the slope change trend is opposite. In iOM, the slope AP also increases as the mean pore radius increases, while in OM, with the increase in mean pore radius, the slope of AP decreases. This is because the low area fraction of OM (10%) limits its contribution to AP (Zhang et al. 2017). The conclusion is consistent with the study of Naraghi and Javadpour (2015). A larger standard deviation (σ) will increase the uncertainty of pore radius and decrease the contribution of pores with larger radius. Therefore, AP decreases as the increase in σ in iOM as shown in Fig. 8a, while the change of σ in OM has minor effects on AP. Figure 9 compares the enhancement factors ( ) when considering nanoconfined effects of different shale components, which is defined as: where is the enhancement factor and q HP is the noslip flow rate. IP is the intrinsic permeability ignoring nanoconfined effects in all shale components. For the cases of considering nanoconfined effects only in iOM, it can be seen that the enhancement factors changes monotonously with the increase in mean pore radius. With the same pore radius, the slip length in quartz nanopores > dolomite nanopores > 0 > calcite nanopores, so ε_quartz > ε_dolomite > 1 > ε_calcite. The difference between them becomes more pronounced with the decrease in mean pore radius in iOM as shown in Fig. 9a, while the change of mean pore radius in OM has minor effects as shown in Fig. 9b. Interestingly, with the decrease in mean pore radius in iOM or OM, the ε and ε_OM show a trend of increasing first and then decreasing. As for Fig. 9a, IP decreases gradually as the decrease in mean pore radius in iOM, but the slip length in OM does not change, so the ε_OM increases gradually at first. At the same time, the slip length in OM is much larger than that in iOM, which leads to the nanoconfind effects in OM plays a leading role. Therefore, the trend of ε is more close to ε_OM. However, the area fraction of OM is only 10%, when the mean pore radius in iOM is small enough, 90% smaller AP grids in iOM play a leading role, which will reduce the contribution of nanoconfined effects in OM pores to AP. Therefore, ε_OM gradually decreases to one with the further decrease in mean pore radius, while ε gradually becomes less than one due to the larger contribution of the nanoconfined effects of calcite pores in iOM components. As for Fig. 9b, with the decrease in mean pore radius in OM, IP gradually decreases, while the slip length in OM gradually increases, so the ε_OM also increases gradually at first. When the mean pore radius in OM is small enough, 90% larger AP grids in iOM reduce the contribution of nanoconfined effects in OM pores to AP. Therefore, ε_OM and ε then start to decrease. Figure 9 (c) shows the relationship between ε and mean pore radius of OM under different TOC. The porosity in each grid is 10%, and all other parameters are listed in Table 2. When ε begins to show the trend as shown in Fig. 9b, the TOC is between 0.3 and 0.4. Figure 9c confirms that the decrease in ε is mainly attributed to the low area fraction of OM. The TOC in different shale plays is no more than 0.4 as cited in "Introduction" section, which may be one of the reasons why the mearsured ε in shale porous media is much smaller than that in single nanotube (Majumder et al. 2005).

Effects of InterP and IntraP
During the shallow burial period of shale, there are large number of InterP and few of IntraP (Wu et al. 2020a).
With the increase in thermal maturation, InterP gradually decreases, while IntraP gradually increases due to the dissolution of acidic fluid. In Fig. 10, we change the ratio of InterP porosity to Intra porosity (R Inter_Intra ) to analyze the effects of InterP and IntraP on AP. A higher ratio means a   lower thermal maturation of shale. The total porosity is kept as 10%. It can be seen that the change of ratio has little effect on the enhancement factor. However, the effect of ratio on AP showed two trends: when the ratio is greater than one, AP increases with the decrease in the ratio, while when the ratio is less than one, AP gradually decreases. This phenomenon indicates that a more uniform distribution of InterP and IntraP porosity will reduce the flow resistance. Accordingly, this may also explain why the appropriate high shale maturity is more conducive to the development of shale oil.
In order to further study the effect of heterogeneity of porosity distribution on AP, we set the porosity in different shale component grids to different values. The generation and expulsion of hydrocarbons will form numerous pores in OM. In a single OM grain, the porosity can even reach 50% (Yao et al. 2018). Therefore, we set the porosity of OM grids as the highest. Loucks et al. (2012) classified iOM minerals based on chemical stability. Quartz is not easy to be dissolved, so it has less IntraP. Calcite and dolomite are chemically unstable minerals, so they have more IntraP. Therefore, we set the porosity in the IntraP grids of quartz as the lowest, followed by the IntraP grids of calcite and dolomite. For a specific iOM shale component, the porosity of InterP grids is still twice that of IntraP grids. The total porosity is still kept as 10%. As shown in Fig. 11, with the increase in porosity ratio, AP gradually decreases. AP with InterP porosity/IntraP porosity=9 Fig. 12 Effects of grain size distribution on AP. a Sensitivity of AP to the mean grain size, the standard deviation is 10μm 2 . b Sensitivity of AP to the standard deviation, the mean grain size is 50μm 2 a porosity ratio of one (6.55μD) is nearly three times as high as AP with a porosity ratio of ten (2.41μD). The results confirm the necessity of dividing iOM into different mineral components in the estimation of AP. If we ignore the difference of porosity between different mineral components, the estimated AP could be erroneous.

Effects of grain size distribution
In this section, we study the effects of GSD on AP. As shown in Fig. 12, when the R Inter_Intra equals to one, GSD has minor effects on AP. These results are consistent with the study of Zhang et al. (2019). In their work, they found the influence of OM distribution on AP is nearly negligible. However, when the R Inter_Intra is greater than one, AP will increase with the increase in the mean and standard deviation of GSD. This is because a GSD with larger mean and standard deviation will have fewer InterP grids, resulting in more uniform porosity distribution. It should also be noted that, different from the discussion in "Effects of pore size distribution and nanoconfined effects" section that the large standard deviation of PSD will reduce the contribution of larger pores to AP, the large standard deviation of GSD would increase the contribution of larger grain to AP. As shown in the green dotted line of Fig. 12a, the AP with the mean grain size of 100μm 2 is nearly 20% larger than that with the mean grain size of 10μm 2 . What is more, the AP with the standard deviation of 100μm 2 is nearly 10% larger than that with the standard deviation of 10μm 2 . Therefore, the effect of GSD should be taken into account to estimate AP correctly, especially when the OM maturation is low.

Conclusions
In this work, based on MDS data, we develop a TFS model to calculate the oil flow rate in shale nanopores. Then, a SD model considering different shale components (OM and iOM) is established with a modified QSGS algorithm. After that, we integrate the TFS model into SD model to construct the REA-scale AP model. Based on the AP model, effects of nanoconfined effects, PSD, pore type, porosity distribution and GSD on AP are studied. The main conclusions of this work are as follows: (1) Both the increase in mean pore radius in iOM and OM will increase AP, but AP is less sensitive to the mean pore radius in OM than in iOM due to the low area fraction of OM. Under the same radius, the slip length in nanopores of different shale components is: OM > > quartz > dolomite > 0 > calcite. Therefore, when TOC is greater than 0.3 ~ 0.4, the main contribu-tion to AP comes from the slip in OM. However, when TOC is less than 0.3 ~ 0.4, the contribution of slip in OM increases first and then decreases with the decrease in mean pore size. This implies the enhancement factor in shale porous media is far less than that in single nanotube.
(2) When the R Inter_Intra is greater than one, AP increases gradually with the decrease in ratio. Therefore, we may safely draw a conclusion that shale oil development is most favorable when the maturity of shale gradually increases to the point where the number of IntraP is comparable to that of InterP. In addition, the heterogeneity of porosity distribution can significantly reduce the AP. If we ignore the difference of porosity between different shale components, the estimated AP could be erroneous.
(3) When the R Inter_Intra is greater than one, both the increase in mean and standard deviation of GSD will increase the AP. The results suggest that the relative error of AP estimation based on average grain size can be up to 20%. Therefore, effects of GSD on AP cannot be ignored, especially when the OM maturation is low.
Funding This study was supported by the National Science Foundation (51804328).

Conflict of interest
On behalf of all the co-authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.