Effects of Slip Length and Inertia on the Permeability of Fracture with Slippery Boundary Condition

Abstract

Although the slippery boundary condition (BC) has been validated to enhance fracture permeability (k), the coupling effects of heterogeneous slippery BC and inertia on k remain less understood. We used computational fluid dynamics to investigate the competing roles of slippery BC and inertial forces in controlling k evolution with increasing pressure gradient by designing six cases with different slip length scenarios for a two-dimensional natural fracture. Our results suggest that pronounced inertial effects were directly related to and demonstrated by the growth of recirculation zone (RZ); this caused flow regimes transitioning from Darcy to non-Darcy and significantly reduced k, with an identical tailing slope for six cases, regardless of the variability in slip lengths. Moreover, the slippery BC dominantly determine the magnitude of k with orders depending on the slip length. Lastly, our study reveals that the specific k evolution path for the case with a varying slip length was significantly different from other cases with a homogeneous one, thus encouraging more efforts in determining the slip length for natural fractures via experiments.

1. Introduction

Fracture flow is critical for the success of many environmental and engineering applications, including the exploration of oil and gas and geothermal energy [1,2], groundwater contamination remediation [3,4,5,6,7], and geological carbon sequestration [8,9,10,11]. Essentially, the fluid-borne contaminants travelling along with subsurface fluid flow are normally hard to predict, where non-Fickian transport phenomenon is pervasive [3,5]. In such a circumstance, a reliable characterization of subsurface flow fields is the key to accurately capture the transport process [12]. For fractures, although the traditional fracture flow theory, namely Cubic Law [13], or the recently modified Local Cubic Law considering fracture heterogeneity [14,15], have been widely used to predict fracture flow, how inertial effects come to play, and change fracture permeability remain less constraint.

Inertial effects cannot be ignored with increasing pressure gradient, driving flow regime transitioning from Darcy to non-Darcy [16,17,18]. Consequently, permeability declines with increasing inertial forces, which makes it difficult to accurately predict (non-Darcy) fracture flow. To overcome the challenge in describing non-Darcy flow in fractured and porous media, a recent study proposed a universal relationship between inertial and viscous permeability based on thousands of laboratory experimental results [19]. With the established universal relationship, one can easily predict fracture nonlinear flow, but with several orders of magnitude of uncertainties that need to be refined.

We proposed that one of possible reasons resulting in large uncertainty is the slippery boundary condition (BC). Slippery BC could be caused either by the chemical alteration or by the coverage of wetting film on the solid surfaces. For instance, many studies reported that the rock surface can be switched from hydrophilic to hydrophobic after its exposure to a non-wetting fluid [20]. The slippery BC, in fact, has been tested to enhance fracture flow (i.e., permeability) via both laboratory and numerical experiments [20,21]. Note that the enhanced fracture flow might not increase the longitudinal dispersion coefficient within the simplified parallel plate model [21,22] for the first time investigated the competition between the inertial force and slippery BC effects in changing fracture permeability, where inertial force on the one hand reduces permeability, while slippery BC, on the other hand, enhances permeability. By employing a homogeneous slip length along fracture boundary, they found out that slippery effects dominated over inertial effects when the pressure gradient was small, and the inertial effects gradually caught up with that caused by the slippery BC with a further increase in pressure gradient [21]. Nonetheless, how heterogeneous slippery BC (i.e., slip length varies along fracture length) alters permeability with increasing pressure gradient remain less understood.

In this study, we aim to address following questions: (1) how does the heterogeneous slippery BC alter the flow regimes transitioning from Darcy to non-Darcy; (2) is there an equivalent homogeneous slip BC that can be used in lien of heterogeneous slippery BC in terms of permeability alteration; (3) how are the recirculation zones (RZs) that are caused by inertial effects developing and growing that could be further complicated by heterogeneous slippery BC. To this end, we resorted to computational fluid dynamics to quantify RZ growth and permeability evolution with increasing pressure gradient.

2. Methodology

2.1. Fracture Data and Computational Fluid Dynamics

For the purpose of gaining fundamental knowledge without losing generality, we used a two-dimensional (2D) fracture (Figure 1a), extracted from a three-dimensional fracture that was obtained from the Digital Rocks Portal (https://www.digitalrocksportal.org/), where the welded Santana tuff sample has been used to investigate the fluid flow and solute transport processes [5,14,21]. The fracture was scanned by X-ray tomography with a resolution of about 0.25 mm done at The University of Texas at Austin, with a length of about 15 cm and mean aperture of about 6.2 × 10⁻⁴ m. More fracture information can be found at reference [14].

Fluid flow through fractures is essentially governed by the Navier–Stokes equations (NSE), along with a continuity equation:

$$\rho \left(u·\nabla \right)u=-\nabla p+\mu {\nabla }^{2}u$$

(1)

$$\nabla ·u=0$$

(2)

where ρ (1000 kg/m³) is the fluid density, u is the velocity vector, μ (0.001 Pa⋅s) is fluid dynamic viscosity, and p is total pressure. The NSE were directly solved, including an inertial term on the left-hand side to precisely capture the growth of recirculation zones (RZs) with increasing pressure gradient.

In this study, the NSE were solved by computational fluid dynamics via the commercial software COMSOL Multiphysics that is built upon the finite element method. To solve the NSE, a fracture domain was discretized with a minimum mesh size of 30 μm around the fracture walls and mesh size of 60 μm towards to the middle of domain (Figure 1b), resulting more accurate numerical solutions (e.g., velocity field in Figure 1a) around the boundaries. Moreover, we applied a specified pressure boundary condition at both upstream and downstream (i.e., Dirichlet boundary condition), driving fluid flowing from left to right (Figure 1a). We sequentially increased the pressure gradient ($\nabla p$ was up to 6000 Pa/m with an uneven interval of about 50 Pa/m) until RZs were fully developed to be able to capture considerable inertial effects (Figure 1c).

In terms of slippery boundary condition (BC), we followed previous study [21] by specifying a slip length (L_s) along the fracture walls. Different from previous study, we considered six cases with different L_s: (1) Case 1: No slip case with L_s equal to 0 m; (2) Case 2: L_s was homogeneous and set to be the minimum value in the case 3; (3) Case 3: L_s varies along the fracture length that depends on the fracture aperture field following the equation proposed by [21]; (4) Case 4: L_s was homogeneous and set to be the arithmetic mean value in the case 3; (5) Case 5: L_s was homogeneous and set to be the median value in the case 3; and, (6) Case 6: L_s was homogeneous and set to be the maximum value in the case 3. The six cases enable us to further distinguish the competing effects that are caused by the inertial and slippery BC on permeability.

2.2. Automatic Recirculation Zone Quantification

Inertial effects can lead to permeability reduction, which is directly related to development and growth of RZs inside the pores and fractures [16,21]. We resorted to the open source of MATLAB code to be able to capture and quantify the growth of RZs (https://figshare.com/articles/Automatic_detection_of_2D_recirculation_zone_volume_/5811399). More details about the delineation of RZs can be found in previous study [23].

2.3. Fracture Permeability Analysis

Upon solving the NSE, we computed fluxes (Q) for above-mentioned six cases by integrating numerically derived velocity at the downstream of fracture. With knowing Q, the permeability (k) can be estimated by:

$$k=\frac{{\left(\frac{12\mu Q}{\nabla p}\right)}^{2/3}}{12}$$

(3)

where $\nabla p$ is pressure gradient; and, the Reynolds number in this study is defined as:

$$Re=\frac{\rho ⟨v⟩b}{\mu }=\frac{\rho Q}{\mu }$$

(4)

where $⟨v⟩$ is the mean velocity within the fracture domain and b is mean fracture aperture.

3. Results and Discussion

3.1. Fracture Nonlinear Flow Behavior

The flow transitioning from Darcy (linear) to non-Darcy (nonlinear) regimes occurs when the inertial effects are non-negligible [17,19]. Consistent with previous studies, our computational results under increasing ∇p clearly show that $\nabla p$ − Q curves deviate from the linear fitting line (Figure 2a), where Q is less than the linear proportionality as Darcy’ s law prescribes due to additional inertial energy dissipation [19]. Noticeably, the largest degree of flow nonlinearity deviating from linear Darcy regime is the no-slip case, as expected, followed by the case 2 with a minimum L_s along fracture length due to the enhancement of flow field by slippery boundary condition (BC), and further followed by the cases 3–6 with spatially varying, maximum, mean, and median L_s, respectively (Figure 2a).

An alternative way to diagnose the flow regime transition is to check the k evolution path with an increase in $\nabla p$. k remains constant if Darcy’s law holds, and decreases due to additional inertial energy dissipation with increasing $\nabla p$, indicating flow behavior transitioning from Darcy to non-Darcy regimes [16]. We estimated k based on the numerically derived Q, which directly suggests the emergence of nonlinear flow regime (Figure 2b). Congruent with results inducing from$\nabla p$ − Q curves, the case 1 with no-slip BC has the least k, followed by the minimum Ls case 2, and by the cases 3−5, with the maximum k for the case 6. Again, the reduction of k with increasing $\nabla p$ for different cases can be explained by the extra inertial forces with the magnitude of k reduction, depending on L_s. Interestingly, the slopes of k − $\nabla p$ were almost identical with a further increase of $\nabla p$ ($\nabla p$ > 1000 Pa/m), regardless of variability in L_s specification, suggesting that the inertial effects were significant and might dominate over slippery BC effects in determining the k evolution paths.

3.2. Inertial Effects on Reducing Permeability

The inertial effects were more pronounced with increasing $\nabla p$ and, thus, led to the overall k reduction behavior (Figure 2b). We resorted to the estimation of recirculation zone (RZ) volume to be able to quantify inertial effects (Figure 1c), since RZ is directly related to inertial effects [16,21]. Expectedly, RZ volume increased with $\nabla p$ where inertial effects gradually increased (Figure 3), leading to k reduction with $\nabla p$ (Figure 2b). Note that the case 1 with no-slip BC eventually had the least RZ volume when $\nabla p$approached to 6500 Pa/m; the RZ volume increased from the case 2 to case 3, with the largest and almost identical RZ volume for the cases 4−6 (Figure 3).

When considering the intricate competition between inertial and slippery BC effects [21], we replotted k and RZ volume against Reynolds number (Re) (Figure 4a,b). Alike k ~ $\nabla p$ relationship (Figure 2b), k ~ Re curves followed the linear tailing behavior with differences depending on L_s specification. However, unlike the $\nabla p$ ~ RZ volume relationship (Figure 3), there was little difference in the Re ~ RZ volume curves between the cases with varying L_s (Figure 4b), suggesting the RZ volume was fundamentally controlled by the fluid dynamics, as quantified by the dimensionless number Re, where both inertial and slippery effects come to play a role in determining RZ volume. This is congruent with previous study [21], where the roles of inertia and slippery BC are competing with each other in determining k. Moreover, the competing roles might also lead to the variability in the threshold of Re for identifying the transition from Darcy to non-Darcy flow regimes [24], which warrants future study.

3.3. Slippery Effects on Enhancing Permeability

The slippery BC can undoubtedly accelerate the flow velocity by lessening flow resistance along the fracture walls [20] and, thus, enhanced fracture permeability with a given $\nabla p$ (Figure 2 and Figure 3). Consequently, the degree of flow nonlinearity with increasing $\nabla p$ was dampened by an increase of L_s (Figure 2a), resulting in the highest k for the case 6, followed by the case 5 till to the case 1 (Figure 2b) with a given $\nabla p$. This is somehow expected, because case 6 with the largest L_s was supposed to have the largest k with a given $\nabla p$, and vice versa for the case 1 with no-slip BC (Figure 2b and Figure 4a).

While expected for the k enhancement with slippery effects, the RZ growth behavior with increasing $\nabla p$ cannot not fully explained by the order of L_s values (Figure 3); that is, the case 2 with the minimum L_s was expected to have less RZ volume than the case 3 with heterogeneous (varying) L_s, which was not the case here. The reversed order of RZ volume in the cases 2−3 can be reasonably explained by the complex interplay between inertial and slippery effects, because: (1) the RZ favorably grows in the small aperture areas right behind the large area zone [25], where L_s was small corresponding to the small aperture, thus possibly leading to an equivalently-effective L_s in enhancing RZ growth for the cases 2−3; (2) the slippery effects can be dampened or overshadowed by the inertial effects [21], thus a large L_s can possibly have a small RZ volume. Although RZ volume demonstrated a complex relationship with $\nabla p$, the rest of fractures have the expected order of RZ volumes following the magnitude of L_s, because inertial effects were still weak when compared to the slippery effects, suggesting that the overall k magnitude was controlled by the L_s (Figure 2, Figure 3, Figure 4).

We plotted k against RZ volume to further demonstrate this argument and clarify this complexity in the k evolution (Figure 5). Clearly, k decreased with RZ volume (i.e., inertial effects) with different trajectories depending on the L_s; that is, the magnitude of k with RZ volume followed the order of L_s, with the largest k for case 6, followed by the cases from 5 to 1. This sequence was not interfered by the reversed order (cases 2 and 3) of RZ volume with increasing L_s (Figure 3), suggesting the dominant role of L_s in controlling k.

Finally, our results reveal that the realistic case (case 3) with the varying L_s along fracture walls cannot be replaced and investigated using a homogeneous with either minimum, mean, median, or maximum values of L_s, thus encouraging more experimental insights into the determination of L_s for natural fractures.

4. Conclusions

We investigate the roles played by the slippery effects and inertial effects with increasing pressure gradient in controlling fracture permeability (k) by designing six numerical experiments with different slip lengths for a two-dimensional fracture. We found that the inertial effects can significantly reduce k with emergent growing recirculation zones (RZs) inside the fracture. Unlike inertial effects, slippery effects can enhance fluid flow and k, thus competition between inertial and slippery effects was demonstrated by the growth of RZs. Moreover, slippery effects overall controls the magnitude of k over the course of k reduction with increasing pressure gradient, where the more realistic varying slip length case cannot be represented by using any homogeneous one. Therefore, urgent effort is required to measure the spatial distribution of slip length via laboratory experiments to accurately predict fluid flow through fractures considering slippery effects.