Apparent Liquid Permeability in Mixed-Wet Shale Permeable Media

Apparent liquid permeability (ALP) in ultra-confined permeable media is primarily governed by the pore confinement and rock-fluid interactions. A new ALP model is required to predict the interactive effect of the above two on the flow in mixed-wet, heterogeneous nanoporous media. This study derives an ALP model and integrates the compiled results from molecular dynamics (MD) simulations, scanning electron microscopy, atomic force microscopy, and mercury injection capillary pressure. The ALP model assumes viscous forces, capillary forces, and liquid slippage in tortuous, rough pore throats. Predictions of the slippage of water and octane are validated against MD data in the literature. In up-scaling the proposed liquid transport model to the representative-elementary-volume scale, we integrate the geological fractals of the shale rock samples including their pore size distribution, pore-throat tortuosity, and pore-surface roughness. Sensitivity results for the ALP indicate that when the pore size is below 100 nm strong pore confinement allows oil to slip in both hydrophobic and hydrophilic pores, yet it restricts the ALP due to the low intrinsic permeability. The ALP reduces to the well-established Carman-Kozeny equation for no-slip viscous flow in a bundle of capillaries, which reveals a distinguishable liquid flow behavior in shales versus conventional rocks. Compared to the Klinkenberg equation, the proposed ALP model reveals an important insight to the similarities and differences between liquid versus gas flow in shales.


Introduction
Flow enhancement of liquids in confined hydrophilic and hydrophobic nanotubes is often observed in experiments where the liquid flow rate is reported to be several orders of magnitude more than that predicted by the classic Hagen-Poiseuille equation [20,74,44,55,57]. Molecular dynamics (MD) simulations are often used to understand the fluid structure and the fast transport mechanisms under confinement [44,39,52,65]. Physical properties (e.g. viscosity and density) of liquid near the tube wall can be different from the bulk liquid due to liquid-solid interactions [56,24]. Fast transport of non-wetting liquid is attributed to the hydrogen bonding of the liquid, which results in the recession of liquid from solid surface [39], the formation of "a nearly frictionless vapor interface" between the surface and the bulk phase [52], or fast ballistic diffusion of liquid [65]. Fast transport of wetting liquid is attributed to the presence of excessive dissolved gas at the liquid-solid interface [20] or the capability of water migrating from one adjacent adsorption site to another [37].
The computational effort of MD simulations can be intensive and time-consuming. Quantitative analytical models have so far been able to predict the flow enhancement of the confined liquid. The flow enhancement factor ( f ), defined as the ratio of the measured (apparent) volumetric flow rate (Q app ) to the intrinsic volumetric flow rate predicted by the Hagen-Poiseuille equation (Q), is usually applied to evaluate flow enhancement through nanotubes. Table 1 summarizes the evolution of some analytical models for flow enhancement. Main differences between these models are how viscosity is modeled and who is the contributor to the flow enhancement. The Tolstoi model [7], one of the earliest quantitative attempt to model liquid slippage along a capillary of radius r, assumes that the average liquid viscosity remains constant along the radial direction of the flow and is not affected by the wall. In reality, the viscosity of the confined fluid is affected by the presence of a gaseous film between the bulk liquid and the wall [20,39,66]. Figure 1 presents a conceptual model of two-viscous flow regions and slippage in a ultra-confined channel or pore. A similar spatially variable viscosity was included by Thomas & McGaughey [66] and Myers [55]. Mattia & Calabrò [54] separate the flow enhancement due to solid-liquid interactions from the depletion region of reduced viscosity. They defined liquid slippage as a function of surface diffusion and work of adhesion on the wall. This definition characterizes the migration of liquid molecules on the surface in addition to the viscous flow. Shale rocks are ultra-confined permeable media with a typical intrinsic permeability of less than 0.1 mD [59,64]. The shale constituents are primarily divided into organic matter and inorganic minerals each having evidently different wettability. The presence of organic and inorganic pores causes the mixed-wettability characteristic of shales. Inorganic clay minerals, e.g., kaolinite, illite, and smectite, are usually hydrophilic, while the wettability of organic matter, e.g., kerogen and bitumes, varies from highly hydrophobic to mixed wet based on rock thermal maturity [59]. Experimental study of oil and brine transport in mixed-wet limestones shows that wetting phase can slip in mixed-wet rock, and that the slip length increases with a decreasing pore size [16]. Recent studies attempted to model the apparent liquid permeability in shale rocks by Table 1 Quantitative models of flow enhancement.

Authors
Flow enhancement models Tolstoi [7] f = 1 + 4l slip r (1a) l slip = δ 0 e αS (Wl−Wls)/kBT − 1 (1b) Myers [55] Nomenclature: l slip is the slip length. δ 0 is the distance between the centers of the neighboring liquid molecules. δ w is the near-wall region thickness. W l and W ls denote the work of adhesion of the liquid and the liquid-solid, respectively. S and α s are the surface area and the fraction of the available sites for liquid migration, respectively. µ(r) is the weighted-average viscosity over the cross-sectional area fraction of the near-wall region (denoted as A w ) and the total flow region (denoted as A t ), where their viscosity is denoted by µ b and µ w , respectively. l slip,∞ is the slip length of a liquid on a flat surface (no tube confinement). C is a fitting parameter. L s is the length of the nanotube (straight length).
applying the aforementioned flow enhancement models [18,86,69,82,68,31,29,17]. Cui et al [18] studied liquid slippage and adsorption in hydrophobic organic pores of shales and highlighted the importance of adsorption layer for oil flow in organic pores of size 500 nm. Zhang et al. [86] modeled liquid slippage in inorganic pores and liquid adsorption in organic pores and predicted that when pore size is less than 10 nm, liquid slippage in inorganic pores is considerable. They implied that wettability differences of inorganic and organic pores can increase apparent permeability of inorganic pores as much as four orders of magnitude more than that of organic pores.
Differences in pore structures of inorganic versus organic matters in shales can also impact transport behaviors [26]. The average diameter of organic pores is usually at least one-order of magnitude less than that of inorganic pores. Organic pores are more uniformly distributed by size than inorganic pores, e.g., 18∼438 nm for organic pores [51] versus 3 nm∼100 µm for inorganic pores [11]. The surface roughness of pores is pore-size scale-dependent: the relative roughness, defined as the ratio of the roughness height divided by pore diameter, is often observed smaller in organic pores than inorganic ones [41]. In the literature, the impact of roughness on liquid slippage is not unique: slippage may be reduced due to stronger hydrogen bonding on rougher surfaces [44,55] or enhanced due to the nano-scale 'lotus effect' [8].
To date, studies on liquid transport behavior especially oil in mixed-wet porous media are limited. Understanding is still insufficient on the overall impact of pore structure and liquid-solid interactions of the rock permeability at the representatively-elementary-volume (REV) scale, i.e., the smallest volume of which the measured permeability and porosity can represent the whole porous medium. This study develops a new apparent liquid permeability (ALP) model based on the experimental and MD simulation results presented in the literature and offers an insight to the microscopic and laboratory-scale data integration for estimating the ALP in a chemical and spatially heterogeneous permeable media. The data include core measurements, scanning electron microscope (SEM) images, mercury injection capillary pressure (MICP), lattice Boltzmann (LB) simulations, MD simulations, and atomic force microscopy (AFM) results. The proposed ALP model presents the following contributions: 1. The ALP model quantifies liquid slippage contribution to flow wetting surface compared to non-wetting surfaces. 2. The ALP model accounts for REV-scale heterogeneity in pore size and pore throat tortuosity, and porescale roughness on liquid slippage. 3. The ALP model clarifies the analogies and differences between the shale permeability model and classic permeability models. A critical comparative analysis between the proposed ALP model and the Carman, the Carman-Kozeny, and the Klinkenberg equations is presented to highlight the advantages and features of the ALP model. 4. The ALP model compiles MD data readily via detailed workflows proposed in this work.

Method
This section presents the derived liquid slippage and ALP model for heterogeneous, tortuous, rough, and mixed-wet porous media at the REV scale. The ALP combines the effect of pore structure, near-wall flow regions, and liquid-pore surface interactions.

Flow enhancement model
In shale rocks, oil exists in a free state and an adsorption state [2]. Adsorption oil is abundant in hydrophobic organic pores, while free oil is stored in both inorganic and organic pores. The flow in the pores consists of two viscosity flow regions: a cylindrical bulk flow region of viscosity µ b and an annular near-wall region of thickness δ w and viscosity µ w . The properties of near-wall regions are different in inorganic and organic pores due to wettability.
In this work, the near-wall region starts from the location where fluid density starts to deviate from the bulk density values and ends at the wall surface. In the literature, we find that the "depletion region" is not consistently defined and sometimes is confused with the "near-wall region" in some bi-viscosity models. A widely-accepted definition of the "depletion region" is the region where local density is less than 2%∼5% of the bulk density [44,30]. The formation of a water depletion region on hydrophobic surfaces is due to the repulsive electrostatic interactions between water molecules and nonpolar surfaces. The typical water depletion thickness was found around one water molecule layer, i.e., ∼2.75 Å [43,53,21,40,61]. Figure 2 shows the velocity profiles of octane in a 5.24-nm silica slit and water in a 2.17-nm carbon nanotube (CNT). In the CNT-water system (Figure 2(b)), at the near-wall region water concentration decreases intensively. This low concentration region from the surface is identified as the depletion region for water. This region corresponds to "velocity peak" and "velocity jumps" in radial and axial velocity [44,62]. In previous flow  [72] and (b) water transport in a 2.17-nm CNT [44]. "DR" denotes the depletion region. . enhancement models for confined water [55,54], δ w is often assumed to be 7 Å, which is more than 2.5 times of the typical water depletion thickness. Through our investigation, the value of 7 Å turns out to reflect a wider region -the near-wall region of fluctuating density [44]. Therefore, the terminology of the "depletion region" may be misleading in the flow enhancement models. Compared to water, velocity jump is not observed for octane in Figure 2(a) and octane viscosity near the solid wall oscillates from positive to negative, i.e., from resisted flow to propelled flow, which makes the viscosity effect on flow contribution less obvious to interpret than water. Liquid slippage is modeled by the Ruckenstein's slip length in Equation (4b), where the key slippage mechanisms are surface diffusion and adhesion. By incorporating Equation (4b), the apparent volumetric flow rate is solved as in Equation (B.1). Flow enhancement is estimated by the ratio of the apparent volumetric flow rate to the intrinsic volumetric flow rate (Equation (A.1)): where is the pore-structure factor and is the slippage factor [66,54]. r p is the pore radius. L s is the straight pore length. W A is the work of adhesion which quantifies the energy of liquid adhesion per solid surface area. We propose the apparent slippage factor λ s,app in Equation (8) to characterize the slippage effect in a tortuous and rough pore. The major algorithm follows: First, λ s is multiplied with a roughness term (1 − ε) 4 to honor the effect of roughness on slippage where ε is the relative roughness presented in Equation (A.7). Then, L s is substituted with the tortuous pore length (L p ) with the tortuosity fractal factor (D T ) in Equation (A.3). Last, the weighted-average pore radius (r p ) over the REV is applied. where and N t is the total number of pores in an REV. Assuming a fractal distribution of pore diameters in the REV, N t can be estimated by γ −D p , where γ = (d p,max /d p,min ) D p is the pore-size heterogeneity coefficient and D p is the pore size fractal dimension. Derivation of N t is referred to Appendix A. Equation (8) is an important modification to Equation (7) because it allows ones to quantify liquid slippage through a non-straight, rough pore throat. When D T = 1 and ε = 0, Equation (8) reduces to the slippage factor for a straight and smooth pore. The tortuosity's impact on slippage factor is further discussed in Section 3.4.

Apparent liquid permeability (ALP)
The intrinsic permeability is characterized by fractal features in pore size distribution, tortuosity, and surface roughness, quantified as where ξ(D T , D p , ε, γ) is the fractal function that combines surface roughness and pore size distribution information. d p,max is the maximum pore diameter in the REV. φ = γ 3−D p is the fractal porosity [84,73]. τ(d p,max ) = (d p,max /L s ) −2D T +2 is the fractal tortuosity of the maximum pore diameter. Pore-scale fractal models and the derivation of Equation (10) are presented in Appendix A. Combining Equation (10) and (5) gives apparent permeability: A modification of λ s to be λ s,app /(1 − ε) 4 yields the final ALP model: The ALP model above is applied to derive apparent permeability in inorganic pores (k i,app ) and organic pores (k o,app ).

ALP estimation workflow
We provide a workflow to estimate APL parameters by using lab experimental and MD results. Table 2 summarizes the data for the key fractal and transport parameters reported in the literature. Figure 3 illustrates this workflow, summarized in three major steps: Step 1. Quantify pore structure to calculate k in Equation (10).
Step 2. Quantify liquid transport, where we model the bulk flow region, the near-wall region, and strength of liquid-solid interactions to calculate f in Equations (5) through (8).
Step 3. Couple Steps 1 and 2 to derive k app as in Equation (13).

Pore size distribution (PSD)
Mercury injection capillary pressure (MICP) test is classically used to estimate PSDs [45]. Figure 3(a) shows the MICP results of a PSD for shale samples. The PSD is divided into two PSDs: a wide-spread inorganic PSD (blue) and a narrow-spread organic PSD (orange). The PSDs are approximated by the pore size fractal distributions.

Pore-surface roughness
Equation (A.7) is used to estimate the relative roughness on pore surfaces. Figure 3(c) is an example SEM image of inorganic matters [83] in shale samples. Figure 3(c) also shows the schematic of the conical nanostructures when the pore surface spreads out as a plane. Key parameters such as the areal ratio (α) and conical height ((h c ) d p ) in Equation (A.7) are estimated via length and areal calculations of the structures observed in the SEM images. The fractal dimension of conical base size distribution (D c ) is estimated via interpreting conical base size and number of cones.  [72]. The density fluctuation signifies the near-wall region (blue) from the bulk flow region (orange). Based on MD results, the thickness fraction of depletion region (δ w /r pi ) and the factor λ bi are estimated by Equation (6). Similar procedure is conducted for estimating parameters for octane transport through an organic pore.  2.3.5 Surface diffusion, work of adhesion, , slippage factor, flow enhancement factor Surface diffusion coefficient (D s ) is derived from MD simulations by evaluating the self-diffusion coefficient parallel with the wall in which the coefficient in the first molecular layer is regarded as the surface diffusion coefficient. Work of adhesion can be obtained via atomic force microscopy (AFM) mapping. Figure 3(e) presents AFM map of force versus distance for tip approach (orange) and withdrawal (blue) [3,34,48]. The encompassed gray area estimates the work of adhesion W A . The apparent slippage factor (λ s,app ) is estimated by D s and W A . The flow enhancement factor ( f ) is calculated based on λ s,app and λ b .

Literature data for confined oil transport
In Table 3, we summarize key transport properties in literature of hydrocarbon liquid transport on hydrophilic and hydrophobic surfaces. Slip length of octane varies widely from 0 to >130 nm in different MD models, which implies a dependence of liquid slippage on the substrate type, driving force, and substrate surface Fig. 3 Flowchart of the ALP model for nanoporous shale [28,25]. (a-c) extraction of pore structure information, i.e., pore size, tortuosity, and surface roughness via of MICP experiments, LB simulations [14], and SEM images [83], respectively. (d-f) quantification of oil transport properties, i.e., near-wall region thickness, viscosity, near-wall region, surface diffusion, and work of adhesion via MD simulations and AFM force mapping, respectively. Intrinsic permeability (k), flow enhancement factor ( f ), and ALP (k app ) are estimated. .
roughness. The thickness fraction of the near-wall region to the slit aperture is also presented in Table 3. The near-wall thickness of octane depends on the pore confinement: In narrow hydrophilic slits, the fluctuation of near-wall viscosity may not stabilize at the slit center, which diminishes the bulk region and cause 2δ w /H → 1. In narrow hydrophobic slits, the bulk-density may not present, which is due to the superimposition of the interaction potentials as well as the adsorption layers from substrate surfaces. Compared to the near-wall thickness in hydrophilic slits, total adsorption layer thickness fraction in hydrophobic slits is more consistent for 2-nm ∼ 5-nm slits, i.e., around 40%∼50% of the entire flow region. A recent MD study [81] indicates that adsorption layer thickness increases with the increase of pore aperture in kerogen, indicating a constant adsorption layer fraction. Still, W A and D s data are generally lacking from the current body of the literature, which may highlight the importance of proper estimation of W A and D s if not accessible. Table 3 Literature data (MD and AFM) for hydrocarbon liquid physical and transport properties on hydrophilic and hydrophobic surfaces.  3 Results

Validations against MD data
Recent ALP models on liquid slippage in shale matrices have shown their ability to predict the enhancement of the confined water transport in straight nanotubes via MD data, yet their capability of predicting realistic liquid (including oil and water) transport in tortuous, rough nanopores is unknown [18,86,29,67]. Here, we validate our model against a series of MD data in the literature for 1. confined octane transport in straight slit pores, estimated by Equation (14); 2. confined liquid transport in tortuous cylindrical pores, estimated by Equation (8); 3. confined liquid transport in rough cylindrical pores, estimated by Equation (8).

Confined octane transport in straight slit pores
As discussed earlier, Ruckenstein's model in Equation (4b) has been proposed and applied to confined water flow through nanotubes [54]. Recently proposed ALP models [18,86,29,67] To validate Equation (14) for shale oil transport, we compiled the MD data in Table 4 for octane flow through a straight, silica slit [72]. In Equation (14), we assume that: 1. L s is the length of the slit in the axial direction. 2. Slit confinement has little impact on the liquid adsorption, i.e., W A is independent of H. 3. D s and η w depend on H according to MD simulation results [22,71,87].
The following algorithm is implemented to estimate l slip for different slit apertures: Step 1. Estimate l slip from the MD velocity profile for the 5.24-nm slit. The value of l slip is estimated by extrapolating the MD velocity beyond the liquid-solid interface until the liquid velocity becomes zero, where l slip = −v slip /(dv/dz) wall , v slip is the slip velocity at the wall, z is the direction perpendicular to the wall [49].
Step 2. Estimate δ w based on the MD density profile for the 5.24-nm slit.
Step 3. Estimate η w for the 5.24-nm slit from MD viscosity profile via either of the following methods. One can estimate η w via averaging the liquid viscosity in the identified the near-wall region from Step 2. An alternative method is to calculate η w from the effective viscosity data (η e f f ) if available. The effective viscosity is the weighted-average based on the fraction of the cross-sectional areas of the bulk flow and the near-wall region, similar to Equation (2b): where A w = 2δ w L and A t = HL are the cross-sectional areas of the near-wall region of thickness δ w and the entire flowing region in an H-aperture slit. η b is assumed as the bulk viscosity. In this way, Step 4. Estimate W A by Equation (14).
Step 6. Estimate D s for different slit apertures. In literature, the self-diffusion coefficient (D sel f ) of water, n-octane, octanol, dimethyl sulfoxide as well as supercritical methane were found to increase with confinement; this relationship can be described in two piece-wise linear relations where the turning point is around 10 nm [22,71,87]. For the cases studied here in Step 5, slit aperture is less than 12 nm, we assume that the linear relation holds for the lateral self-diffusion of octane at the silica surface, i.e., surface diffusion. Therefore, two D s data including the obtained D s (H=5.24 nm) are required to estimate the linear relation of D s versus H. For example, 7.61 nm slit is chosen to estimate the second D s data. Given W A , δ w , l slip in Table 4 we estimate D s (H=7.61 nm) = 4.24e-9 m2/s. Now with D s (H=5.24 nm), the relation of D s = 0.57H − 0.1 is obtained, where H is in nm and D s is in 1e-9 m 2 /s. A series of D s for H ≤12 nm is estimated accordingly.
Step 7. Repeat Step 2 and 3 to estimate δ w and η w for different slit apertures based on their density profiles.
Density data can be found in Ref. [72].
Step 8. Calculate l slip for different slit apertures by Equation 14. In Figure 4(a), the estimated l slip from Step 8 is plotted against MD predictions from Step 1. Equation (14) can slightly overestimate the slip length which may be due to the simplified approximation of D s values; direct calculations from MD for different apertures, if available, can improve the reliability of the D s input for Equation (14). Yet given that the differences between MD data and predictions are small (≤ 9%), Equation (14) can reasonably capture the octane slippage in silica slits of apertures ≤ 12 nm.

Confined liquid transport in tortuous cylindrical pores
We propose Equation (8) to estimate apparent liquid slippage on tortuous, rough cylindrical nanopores. Assuming that the impact of tube roughness ε on liquid slippage (λ s,app ) is much smaller than tube tortuosity, Equation (8) reduces to a function of the tortuosity dimension (D T ). To ensure that D T is physically meaningful, i.e., 1 ≤ D T ≤ 3, Equation (8) is valid if the sample tortuosity Therefore, the validity of Equation (8)  To testify its ability to predict liquid slippage in tortuous pores, we adopt MD data for confined water transport in bent and tilted CNTs. The following algorithm is performed: Step 1. Estimate τ via bending angles (α * ) or tilting angles (β * ): τ = (sin(α * /2)) −2 or (sin(β * )) −2 .
Step 3. Estimate D T via τ: D T = 1 − ln(τ) 2ln(d p /L s ) . Fig. 4 (a) Comparison of the MD data [72] for octane transport in silica slit versus predictions from Equation (14). Relative differences are shown for prediction deviations from the MD data. (b) Comparison of the MD data versus λ s,app predictions from Equation (8). MD dataset 1 [58] corresponds to CNT 1 with a bending angle α * . Tortuous length and tube size are fixed as L p = 3.8 nm and d p = 0.777 nm. Different tube tortuosity is achieved via alternating α * . MD dataset 2 [50] corresponds to CNT 2 with a tilting angle β * . Straight CNT configuration is shown as CNT 3. Tortuous length and tube size are fixed as L p = 3.824 nm and d p = 0.782 nm. Different tube tortuosity is achieved via alternating β * . Both models were conducted under 300 K. Bulk viscosity is η b ≈0.85 mPa·s [46]. Work of adhesion W A =97 mJ/m 2 and surface diffusion coefficient D s = 4e-9 m 2 /s are used [54]. (c) Comparison of the MD data [60], a slip model for smooth CNTs [86], versus λ s,app predictions from Equation (8).
Step 5. Estimate λ s,app for tortuous CNT type 1&2 by Equation (8) with inputs of D T from Step 3. (8) (with ε → 0) against MD data. Detailed descriptions of the dataset are captioned in Figure 4(b). The result shows that Equation (8) is a good predictor of the slippage factor in tortuous nanopores. For the CNT type 1&2, d p and L p are known as constant and τ = (sin(α * /2)) −2 or (sin(β * )) −2 . To meet the requirement of 1 ≤ D T ≤ 3, the bending angle α * ≥ 360 π arcsin[(  Figure 4(b)) follows the trend of the MD data, Equation (8) should not be used to predict the slippage factor.

Confined liquid transport in rough cylindrical pores
Another feature of Equation (8) is the consideration of nanopore roughness. Relative roughness is a scaledependent parameter, the importance of which can increase with the decrease of pore diameter. In Equation (8), roughness in modeled as resistance to liquid slippage. To validate this resistance effect, we adopt the MD data for confined water in straight, rough CNTs. Figure 4(c) compares the MD data, the apparent slippage factor estimated by Equation (8) (with D T = 1 for the straight CNT), and a slippage model for smooth CNTs [86]. The results show that the apparent slippage factor agrees with the MD results better than the model for smooth CNTs, which signifies the pore-surface roughness impact. Relative roughness is estimated to be ε = 0.07 through data matching.

Governing factors of confined oil transport
Input data are adopted from experimental and MD results [23,75,34,45,54,79,14,78] and presented in Table  2. Some parameters of ALP are estimated by the equations listed in Table 2.  Figure  5 shows that a thicker near-wall region in inorganic pores improves the flow capability while in organic pores it reduces the flow capability, although their influences are generally small. This is expected as interactions between oil and inorganic pore surfaces are weaker than oil and organic pore surfaces. Figures 6(a) and 6(c) respectively present the impact of D s on the ALP in inorganic and organic pores. The ALP increases with the rise of D s . Oil slippage is more evident in inorganic pores (with a higher D s and a lower W A ) than in organic pores of the same diameter, which qualitatively agrees with MD simulations of octane through muscovite and kerogen pores [36]. Figures 6(b) and 6(d) show the impact of W A on the ALP for inorganic and organic pores. The ALP decreases with the increase of W A due to strong adhesion between liquid and pore surface and the resultant weaker slippage.

Pore confinement & surface roughness
The impact of pore confinement on slippage is demonstrated in Figure 7 66/1 = 66 can enhance the slippage factor by 10 folds for a pore radius of 1 nm-100 nm. Slippage is also influenced by pore size. Slippage becomes more noticeable when the pore radius decreases. Figure 7(b) shows that the slippage factor decreases exponentially as the pore radius increases. When the pore radius reaches 100 nm, slippage decreases until no flow enhancement is observed (λ s,app → 1). The sensitivity of ALP to surface roughness depends on pore type since organic pore surfaces are generally smoother and is more uniform than inorganic ones. The "resistance" effect of surface roughness is therefore not as evident in smoother organic pores as in rougher inorganic pores. Figure 7(a) presents the impact of relative roughness on the apparent slippage factor (λ s,app ). A higher relative roughness entails a lower λ s,app . The lower slippage of fluid along rough surfaces is because the molecules in the fluid might prefer rougher surfaces or that the fluid is "pinned to" the irregular wall surface [33]. Figure 8 shows the overall effect of pore confinement on the apparent and intrinsic liquid permeability, where important observations follow.

Apparent versus intrinsic liquid permeability
1. With decrease in pore tortuosity and increase in pore size, the intrinsic permeability increases. 2. When the pore size increases, the gap between apparent permeability (k app in solid lines) and intrinsic permeability (k in dashed lines) becomes narrower and lines eventually overlap. This implies that the effect of pore size on flow enhancement decreases until it is non-existent as the pore size increases. 3. When the pore tortuosity is increased, the gap between apparent permeability and intrinsic permeability decreases for a certain pore size. This suggests that liquid slippage is weakened. For the most tortuous condition, the apparent permeability is enhanced to its maximum value, even for larger pore radii.
In Figure 8, the point at which lines start to overlap marks the onset of reduced slippage contribution. Flow enhancement due to liquid slippage is effective when r p < 100 nm. Pore radius starts to exert a positive effect on apparent permeability when r p reaches 100 nm. For a straight pore (τ = 1) and intermediatelytortuous case (τ = 66), this positive effect is honored. For the highly-tortuous case (τ = 4518), the effect of r p is reversed when r p ≥ 100 nm: A slight decrease of the solid line is observed, rather than an increase Presented are two representative pore-confinement conditions with average pore sizes and average pore tortuosities: a weakly confined pore (top) and a strongly confined pore (bottom). We assume that the representative pore confinement can reflect the intrinsic permeability of porous media.
in the less confined cases. Also, pores with lower intrinsic permeability always has lower apparent permeability. This is because although considerable slippage occurs in highly-confined pores, the strong effect of confinement on intrinsic permeability limits the effect of slippage on its apparent permeability.

Liquid transport mechanisms in shales
Multiple structural and fluid factors affect apparent liquid permeability and liquid slippage. Equation (8) indicates that liquid slippage is affected by macroscopic liquid-solid interaction parameters, D s and W A , local pore confinement properties, such as average pore size (r p ) and pore tortuosity (τ). Pore confinement has opposite effects on intrinsic permeability and liquid slippage. Confinement decreases intrinsic permeability, but increases slippage. However, strong liquid slippage may not boost high apparent permeability when the shale intrinsic permeability is in tens of nanodarcy or less range.
Dissimilarities in wettability, average pore size, and pore tortuosity for pores of different types in mixedwet porous media lead to a complicated slippage mechanism, intrinsic permeability, and apparent permeability. Table 5 summarizes the sensitivity conclusions drawn from this work. The governing factors of ALP are the pore size distribution and pore throat tortuosity.
Liquid viscosity near the pore wall can be different from the center of the pore due to the solid-liquid interactions. If the work of adhesion is strong enough where liquid tends to stick to the pore surface, the nearwall viscosity is higher than the viscosity in the pore center. The sensitivity results indicate that viscosity variation near the pore wall may not have a significant impact on the flow enhancement unless the pore diameter is ultra-small, i.e. within an order of liquid molecule size. Given that in the shales, the largest connected pores have the most share of the contribution to the overall flow, one may conclude that the fluid viscosity change near the pore wall has a negligible effect compared to surface diffusion and wettability. A quantitative comparison between the estimated range for the intrinsic permeabilities and apparent permeabilities (by using Equations (10) and (13), and Table 2) suggests that flow enhancement, mostly due to liquid slippage, can reach nearly 300 in both wetting and non-wetting pores. Here, the dual effect of liquid-solid interaction (wettability, adhesion, and surface diffusion) and pore confinement (pore size and pore tortuosity) renders a comparable flow enhancement in inorganic and organic pores.

The ALP model comparison with the Carman & the Carman-Kozeny equation
It is instructive to understand the relation between the ALP and the fluid equations that predict the pressure drop of fluids through permeable media. We investigate two classic equations, namely Carman [10] and Carman-Kozeny [9,1] and show that under reasonable assumptions, the ALP will reduce to the spirit of these two classic equations.
To derive the Carman equation, we begin with a simple version of the ALP, as in Equation (A.6). Apply the assumptions of the Carman equation to a constant viscosity distribution (µ b = µ w ), no surface diffusion (D s = 0), no pore roughness (ε = 0), and uniform pore diameter (d p ≡ d p,max ) and set the fractal parameters to D T = 1 and D p = 2, the ALP reduces to the Carman equation (Equation (A.2)) in the limit.
The Carman-Kozeny equation is used for predicting a fluid flowing through permeable media packed with spherical, smooth, and solid grains. A generalized version of the CarmanâȂŞKozeny equation is [1] where k CK is the Carman-Kozeny permeability, F s is the pore shape factor, S gv is the ratio of grain surface area to the grain volume (S gv = 4φ/d p (1 − φ) for spherical grains). Equation (16) accounts for the porosity and geometric properties of grain and pore. The product of F s τ is referred to as the Kozeny constant and is a strong function of grain size distribution. The Kozeny constant is often fitted to the experimental data to obtain the best predictor of permeability based on the porosity for different hydraulic units. Equation (10) has the spirit of Carman-Kozeny equation in Equation (16) where the fractal function (ξ) is the inverse of half the pore shape factor, i.e., ξ = 2/F s . For a bundle of identical straight cylinders, i.e., F s = 2 and correspondingly ξ = 1, Equation (10) becomes identical to Equation (16).

The ALP versus Klinkenberg equation
A fundamental insight into the fluid flow in nano-scale confined medium such as shale rocks, is the presence of fluid slippage, regardless of the phase type. Gas slippage, also known as the Klinkenberg effect occurs due to the rarefaction. Gas rarefaction is caused by decrease of gas pressure, reduction of characteristic length or pore size. Either of these factors increases the dimensionless Knudsen number (Kn). With the increase of the Knudsen, gas flow becomes more rarefied and transits from slip flow to transitional flow, and eventually to the free molecular flow [42,19,27]. As for liquid flow, the mean molecular free path is far less than the pore size, and the Knudsen number is small. Therefore, the Knudsen number does not characterize the transition of the liquid flow regimes. The work of adhesion for liquid has a similar role as the Knudsen number for gas. Work of adhesion characterizes the liquid-solid interactions, a phenomenon that governs liquid slippage in confined pores. Work of adhesion quantifies the energy required to overcome free energies per area of three-phase interfaces of liquid-solid, solid-vapor, and liquid-vapor and subsequently separate liquid phase from the solid phase, W A ≈ γ LV (1 + cosθ) [89], where γ LV is surface tension between liquid and saturated vapor in the unit of N/m, and θ is the contact angle between liquid-vapor and liquid-solid interfaces. A small value of θ indicates an increased work of adhesion, which is a characteristic of a liquid-wetting solid surface. A liquid-wetting surface requires more energy per surface area to separate the liquid from the local site, the liquid is less willing to flow near the wetting surface, and subsequently slippage is small. The impact of pore surface wettability is strong in ultra-confined pores because the fluid flow inevitably occurs near the pore wall. The interfacial interaction between the liquid and vapor phase, quantified by γ LV , is related to the near-surface fluid viscosity, contact angle (wettability), and surface diffusion. The cause-and-effect relationship between these parameters is not identified to this date; nevertheless, the liquid slippage phenomenon is a synergic effect of all the above parameters.
Based on the sensitivity results, we find that the effect of the viscosity near the wall is much smaller than that of the surface diffusion and the work of adhesion. We, therefore, drop the flow contribution from the viscosity. Similarly, the roughness term can also be dropped as pore surface roughness is less dominant than the pore size and tortuosity effect. Given that the fractal intrinsic permeability can essentially represent the Carman-Kozeny equation, the derived ALP in Equation (13) can be arranged in terms of k CK and τ as The ALP model (Equation (17)) for liquid presents some interesting analogies as the Klinkenberg equation for gas [47]. First, the liquid slippage is inversely proportional to W A whereas the gas slippage is inversely proportional to the gas pressure. Second, the term b, defined here as the liquid slippage constant, determines the flow enhancement contribution upon the intrinsic k CK . It is a function of pore confinement, i.e., the average pore diameter and the average tortuosity, and surface diffusion of liquid. Interestingly, the term b is similar to the gas slippage constant (b ) in the Klinkenberg equation as b is found to be a strong function of tortuosity and the tangential momentum accommodation coefficient (TMAC) [77], where the TMAC characterizes the how gas molecules are reflected in terms of diffuse reflection and specular reflection on the wall after the gas-wall collision [4].
In a more general manner, the gas slippage is proportional to the Knudsen number (Kn) as a ratio of the mean free path (λ) of the gas to the average pore diameter (d p ). For liquid, λ is much smaller than d p , Kn cannot be a proper indicator of liquid slippage. By assuming that liquid would follow the hydraulic pathways of the pores, we define a dimensionless parameter, the liquid confinement number (Cn), as the ratio of the tortuous path length to the characteristic length, to characterize the liquid slippage. In practice, the average pore diameter is applied as the characteristic length, therefore Cn = L s √ τ/d p . Equation (17) then reads k CK,app = k CK (1 + α · Cn) where the dimensionless parameter quantifies the surface diffusion of liquid of viscosity µ b in the a straight pore of the diameter d p by overcoming the work of adhesion W A . Equation (19) shares the similar structure as the gas slippage model due to gas rarefaction [6]. The derived ALP model in Equation (13) along with its transformation in Equations (17) and (19) delivers a more comprehensive description of the liquid flow in tortuous, heterogeneous porous media, and under proper restricting assumptions, reduces to the spirit of the Carman-Kozeny equation.

Conclusions
Liquid transport in shale rocks is governed by local pore confinement, liquid-solid interaction, and poresurface roughness. We proposed an apparent liquid permeability (ALP) model for heterogeneous and rough nanoporous shale matrices, and a workflow for the ALP estimation. Major conclusions follow: Substituting A in Equation (A.6) with (A.9), one can derive the intrinsic permeability as in Equation (10).

B Derivation of flow enhancement
The apparent volumetric flow rate in a pore is solved as [86] where D s is the surface diffusion coefficient. A high value of D s reflects a fast diffusion of liquid molecules on the surface. Measured D s values for oil on different wettability surfaces are reported in the order of 1e-9 m 2 /s to 1e-8 m 2 /s [23,54]. W A is the work of adhesion. A low value of W A implies low adhesion, and a simultaneously strong slippage at the liquid-solid interface. The value of W A ranges from 3e-3 J/m 2 for a strongly non-wetting surface, to 3.5e-1 J/m 2 for a strongly wetting surface [34,54]. The flow enhancement is derived through deviding Q app by the intrinsic flow rate (Q), where Q is derived by incorporating δ w = 0 and D s = 0 in Equation (B.1):