Recirculation zones (RZs) in rock fractures have been widely observed by experiments and numerical simulations. While previous studies focused on the effects of RZs on flow regimes and solute transport, limited attention has been given to their evolution across a wide range of flow velocities and the associated impacts on fracture permeability. In this study, numerical simulations were conducted to investigate the evolution of RZs over a wide range of Reynolds numbers (Re) and their effects on the viscous (kv) and inertial (ki) permeabilities of single fractures. A three-stage evolution of RZ across a wide Re range was detected: Stage I (rapid growth): During the initial formation of RZs, their volume (Sv′) increases rapidly with Re; Stage II (slow growth): As Re increases, Sv′ continues to grow, but dSv′/dRe gradually decreases. Stage III (fully developed): At higher Re, Sv′ becomes insensitive to further increases in Re, with dSv′/dRe ≈ 0. During the transition from Stage I to Stage II, the expanding Sv′ compresses the main flow channel (MFC), reducing its nonlinearity. This leads to a decrease in viscous permeability (kv) and an increase in inertial permeability (ki) as Re increases. In Stage III, RZs become fully developed and independent of Re, resulting in stable kv and ki as RZs and MFCs reach a highly differentiated and stable configuration. A critical Re (Rec,stable) was defined to capture the stable kv and ki, referred to as kvglobal and kiglobal, respectively, encapsulating the overall evolution of hydraulic conductivity in rock fractures. Additionally, quantitative models for kvglobal and kiglobal were derived and validated.

Understanding residual non‐aqueous phase liquid (NAPL) architectures in fractures is crucial for predicting NAPL dissolution and depletion in fractured aquifers. Although NAPL infiltration and dissolution in porous media have been widely studied, residual NAPL architectures and associated mechanisms in fractures remain poorly understood. To address this, we conducted 162 water‐NAPL displacement experiments using microfluidics. Three distinct residual NAPL architectures were identified, Pools Pattern (PP), Ganglia Pattern (GP), and Mixed Pattern (MP), affected by aperture heterogeneities and flow rates. Experiments combined with numerical simulations revealed that aperture anisotropy, roughness and flow rate determine balance between local capillary and viscous forces, driving NAPL entrapment mechanisms that form observed architectures. Capillary‐trapping leads to NAPL pools (PP), viscous‐trapping leads to ganglia (GP), and a combination of both mechanisms leads to MP. This work elucidates the mechanism of residual NAPL architectures in fractures and lays a foundation for modeling residual NAPL dissolution in fractured aquifers.

Understanding the anomalous solute transport in single fractures is important for many hydrogeologic processes and subsurface applications. Recirculation zones (RZs) and corresponding main flow zones (MFZs) have been widely recognized as low-velocity regions and preferential pathways that could explain the simple anomalous solute transport, i.e., heavy tailings and early arrival. However, the direct relation between RZs and more complex anomalous transport phenomena, e.g., multi-modal peaks and fluctuating tailings, has been elusive. This may be due to the limited understanding of the evolution of RZs and the mass transfer process between RZs and MFZs, i.e., the monotonically increasing RZs volume (Sv) and defaulted diffusion-dominated mass transfer. In this study, we systematically generate a series of 2D/3D rough single fractures with different geometric properties to investigate the evolution of RZs and its influence on anomalous transport across a wide Re range of 0–426.88. Three-stage evolution of RZs with increasing Re was identified by using the growth rate of Sv (dSv/dRe), the rapid growth stage (Stage I) where dSv/dRe increase, the slow growth stage (Stage II) where dSv/dRe decrease, and the fully developed stage (Stage III) where dSv/dRe is a constant. The mass transfer mode between recirculation and main flow zones is shifted from diffusion-dominated in Stage I to convection-dominated in Stage II due to the enhanced convection in RZs. This shift of mass transfer mode enhances the mass transfer rate (α) between RZs and MFZs by 5–20 times. In Stage II, the solute was trapped around the interface between RZs and MFZs before entering RZs, i.e., the solute “film”. The coexistence of the solute “film” and the solutes trapped by RZs induces multi-modal peaks and strengthened tailings of BTCs. In Stage III, the solute “film” cannot form due to the rapid dissipation of detained solutes driven by stronger convection-dominated mass transfer around the RZs-MFZs interface, which in turn leads to the disappearance of multi-modal peaks and induces monotonically shortened tailings. This study fills the gap in the RZs evolution and the associated mass transfer process in the microscopic flow fields, which deepens our understanding of the anomalous transport mechanism.

Two-phase immiscible flow displacement in rock fractures is important for many subsurface engineering applications. When one fluid displaces another more viscous one, the displacement regime transitions from capillary fingering to viscous fingering with increasing flow rate. This displacement regime transition is widely investigated in rough fractures. However, how the variable aperture space of rough fractures influences the regime transition and associated mechanism remains unclear, especially for the anisotropic aperture space. In this work, we investigated the influence of the aperture field anisotropy on the regime transition via flow rate-controlled drainage experiments and simulations. By analyzing the pore-scale fluid–fluid interface advancements, we found that the frequency of the transverse pore-filling events (TPFEs) determined by the two-way coupling dynamics between flow rate and aperture anisotropy, and thus controlled the regime transition. For transverse-correlated aperture fields, the frequency of TPFEs was enhanced, which coupled with the non-monotonic effect of flow velocity on TPFEs, promoted the regime transition from capillary to viscous fingering with the increasing flow rate. Conversely, the longitudinal-correlated aperture fields reduced the frequency of TPFEs and suppressed the regime transition. Moreover, we obtained the theoretical models of the critical capillary number (Ca) corresponding to the regime transition as a function of aperture anisotropy. We also modified the phase diagram for fractured media by considering the aperture anisotropy. This work bridges the gap between fracture geometric attributes at the Darcy scale and pore-scale displacement processes, which provides the foundation for upscaling modeling of multiphase flow in the future.

Two-phase flow displacement in rock fractures is crucial for various subsurface mass transfer processes and engineering applications. In fractures, the displacement of a less viscous fluid by a more viscous one (i.e., viscosity ratio M > 1) involves viscous forces help stabilizing the displacement front in presence of capillary pressure fluctuations. Although previous studies have reported displacement patterns in isotropic fractures, the impact of anisotropic fractures on displacement patterns has not been systematically examined. In this study, we conducted flow-rate-controlled drainage experiments to examine how anisotropic aperture fields affect displacement patterns. We observed the transition of displacement patterns from capillary fingering (CF) to crossover zone (CZ) to compact displacement pattern (CD) based on variations in transverse pore-filling event (TPFE) frequency, which characterizes the competition between capillary and viscous forces. Increasing aperture correlation length in the transverse direction leads to increased TPFE frequency at a low flow rate, destabilizing displacement front. While the increasing aperture correlation length in longitudinal direction suppressed TPFE frequency, stabilizing displacement front. Therefore, the critical capillary number (CaCF-CZ), which indicates the onset of the CF-CZ transition, decreases as the aperture field varies from transversely to longitudinally correlated. At high flow rates, TPFEs almost disappeared, indicating that anisotropy did not affect CZ-CD transition (CaCZ-CD). Furthermore, we modified theoretical models of CaCF-CZ and CaCZ-CD by incorporating the aperture anisotropy factor, achieving a good fit with the experimental data. This study demonstrates the critical role of aperture field anisotropy in controlling two-phase displacement patterns and provides a theoretical framework for predicting multiphase flow behavior in natural fractures.

Spectral induced polarization (SIP) exhibits potential to be a nonintrusive approach to monitor bacterial activity in biological hotspots associated with the critical zone of the earth. The polarization of bacteria in a low-frequency electrical field is related to the polarization of their electrical double layer coating their surface. However, few studies have quantified the induced polarization responses on both gram-negative (GN) and gram-positive (GP) bacteria in soil column experiments. To address this gap, 17 experiments using two strains, Pseudomonas aeruginosa O1 (PAO1, GN) and Brevibacillus centrosporus (L3, GP) are conducted. Complex conductivity spectra are collected in the frequency range 10 mHz-10 kHz during bacterial growth and decay phases in soils. The complex conductivity spectra are fitted using a double Cole-Cole model to remove the effect of Maxwell-Wagner polarization. The change in the magnitude of the polarization (quadrature conductivity or normalized chargeability of the low-frequency contribution) is linearly related to the bacterial density, regardless of the type of bacteria. The changes in the normalized chargeability and Cole-Cole relaxation time are directly proportional to the density of bacteria. Furthermore, it is inferred that the thickness of microcolonies plays a critical role in the relaxation time rather than the diameter of individual bacteria. This study expands the potential of SIP for in situ monitoring of microbial activity in soils.

The induced polarization (IP) method holds a strong potential to better characterize the critical zone of our planet especially in areas characterized by multi-phase flow. Power-law relationships between the bulk, surface, and quadrature conductivities versus the pore water saturation are potentially useable to map the subsurface water content distribution. However, the saturation exponents n and p in these power-law relationships have been observed to vary with the texture of geomaterials and the wettabilities of pore fluids. Traditional experimental setups in the laboratory do not allow to independently visualize the pore fluid distribution. Therefore, the physical interpretations of the two saturation exponents have remained unclear. We developed a novel milli-fluidic micromodel using clay-coated glass beads that exhibit excellent visibility and high IP response. Through laboratory experiments, we simultaneously determined the micromodel complex conductivity and acquired the corresponding pore-scale fluid distributions generated by drainage and imbibition through such class of porous materials. Finite-element simulations of complex conductivity based on the upscaling of the complex surface conductance of grains were conducted to determine the saturation exponents under ideal pore fluid distributions. Results indicate that saturation exponents n and p vary depending on the ganglia size of the insulating fluids. The saturation exponents n and p exhibit power-law relationships with the change rate of pore water connectivity with saturation, which is calculated through the computation of the derivative of Euler characteristics. These findings provide a new physical explanation to the relationships between the saturation exponents and the microscopic fluid distributions within the geomaterials. Water saturation of porous bodies can be related to the complex conductivity through power-law relationships. The existence of these power-law relationships has been clearly documented in the literature. They are critical in the realm of hydrogeophysics to better characterize the critical zone of the solid Earth. However, the values of the saturation exponents in these power-laws have been observed to vary with different material textures and no underlying mechanisms have been able to explain these variations to date. This lack of physical understanding could limit the applicability of the induced polarization method to characterize the critical zone especially when immiscible fluid phases are present. To tackle this, we developed a milli-fluidic pore model that allows investigating the saturation exponents while monitoring the pore-scale fluid distributions. Together with numerical simulations for ideal fluid distribution cases, we found a relationship between the saturation exponents and a microscopic pore parameter called "change rate of the pore water connectivity with saturation." These findings suggest that when estimating subsurface water saturation from electrical parameters, taking the pore fluid distribution into account can significantly improve the estimation accuracy, therefore enhancing the efficiency of geo-electrical applications for a better characterization of hydrocarbon contaminated aquifers and oil reservoirs. A micromodel setup allows for simultaneous pore fluid visualization and complex electrical conductivity measurement The two saturation exponents for the in-phase and quadrature conductivities are correlated with the ganglia size of the insulating phase A power-law-type relationship between the saturation exponents and the change rate of water connectivity with saturation is demonstrated

The fracture-matrix system, in which water is stored and transported, has a significant impact on the hydraulic behaviors of fractured geologic media (FGM). However, it is very challenging to accurately simulate flow behavior in FGM due to the difficulty of characterizing dual-permeability media, including multiscale fractures with high permeability and porous matrix with low permeability. In this study, a multiscale fracture integrated equivalent porous medium (MFEPM) method is proposed for simulating fluid flow and solute transport in a fracture-matrix system. The synthetic enhanced matrix (SEM) is compounded by integrating small-scale fractures into the porous medium, and the medium-and large-scale fractures are mapped to refined grids of the MFEPM to increase the characterization of the fracture geometry and orientation. Then, the equivalent hydraulic properties are calculated according to the properties of the fractures. Finally, a finite difference method (FDM) based on the MFEPM is used to simulate flow and solute transport in coupled multiscale fractures and rock matrix. The case study illustrates that compared with the traditional equivalent porous medium (EPM) method, the MFEPM manifests an efficient effect in the characterization of preferential flow; compared with the discrete fracture network (DFN) method, the MFEPM can characterize flow and solute transport in the SEM and represent mass interactions between the fractures and matrix. The results of the numerical case studies of flow and solute transport also show that the MFEPM method has a high computational accuracy and good reliability, while maintaining an appropriate computational burden.

Fracture intersections are an essential element of fracture networks in energy extraction. The accurate evaluation of the flow behavior at intersections in the furcating fracture is the key to describing the flow behavior of fractured reservoirs, yet fracture intersections are often oversimplified or even ignored. To fill this knowledge gap, the effect of geometric characteristics internal to the intersection was investigated based on direct numerical simulation by solving Navier-Stokes equations. Two parameters, 491 and 492, which are the angles between the normal inlet branch direction and the two sides of a closed triangle, are introduced to characterize the geometry of the intersection in a furcating fracture, and the role of geometrical parameters 491 and 492 in the interference effect and flow rate distribution is investigated. The results demonstrate that intersection geometry interferes with the fluid flow passing through the intersection to a certain range. With an increasing hydraulic gradient J and 49, the flow interference range increases. Under the conditions of these simulations (J = 10-3 - 10-2, Re = 101-102), J is linearly correlated with flow rate q, but with an increasing J and 49, the area and percentage of the low-velocity zone caused by an intersection increases leading to the enhancement of nonlinear flow. Furthermore, based on the simulation results, an empirical formula for quantifying the flow redistribution at intersections is proposed, contributing to the effectiveness and accuracy of modeling seepage in the fracture networks. Finally, the physical meanings of geometric parameters 491 and 492 are determined. This study advances the understanding of flow behavior at intersections commonly observed in discrete fracture networks from a fresh perspective, despite the difficulty of monitoring.

The Izbash equation has been widely used in the subsurface applications. However, the Izbash equation is still empirical, and its coefficients (scaling factor λ and power exponent M) have not been systematically characterized and quantified. In this study, laboratory experiments and numerical simulations of fluid flow across a wide range of hydraulic gradients (J = 0–4) in horizontal rough fractures were conducted to comprehensively characterize and quantify the influence of fracture geometric attributes and fluid inertial effects on λ and M. The results showed that λ increased with fracture relative roughness (RSD). The fluid inertial effect (quantified by the non-Darcy effect factor E and Re) had a two-stage influence on λ. When the fluid flow was laminar, λ increased with E. However, when the fluid flow regime starts to transition from laminar flow to turbulent flow, λ decreased with increasing E. M is positively correlated with RSD and the fluid inertia effect E. We found that the transition of flow regime from laminar to turbulent flow depended on whether the recirculation zones are fully developed. The fully developed recirculation zones determine the distortions of the velocity field and flow field, which induced the turbulent flow. The quantitative models of λ and M were obtained based on numerical simulations, which quantified the coupling influence of the fracture geometric property and fluid inertial effect. The validity of quantitative models was verified by laboratory experiments. Our work provided a new understanding of the Izbash coefficients and laid a foundation for theoretical background exploration of the Izbash equation.

This study experimentally and numerically investigated the influence of the inertial effect of fluid flow on the non-Darcy coefficient (beta). The results showed that the non-Darcy coefficient (beta), apparent permeability (k(a)), and hydraulic aperture (e(h)) decreased with the increase of Reynolds number (Re), which indicated that the non-Darcy coefficient depended on the geometric properties of a single fracture and the fluid inertial effect. The quantification model of the non-Darcy coefficient was improved. Under the guidance of fluid inertial effect, the beta performed a positive power-law relationship with e(h); under the coupling effect of fluid inertial effect and fracture geometric properties, the beta behaved negative power-law relationship with e(h), which was similar to the non-Darcy coefficient quantitative model that ignored the fluid inertial effect in the previous study. Therefore, the geometric properties were the dominant ones for fracture flow when it was affected by the coupling effect of geometric properties and inertial effects. Then, the quantitative model of the non-Darcy coefficient was improved. Based on the improved non-Darcy coefficient quantitative model, a new critical Reynolds number prediction model (CRN model) was constructed. Compared with other CRN models, the new CRN model considering the coupling effect of fluid and media properties could more accurately predict the critical Reynolds number (Re-c) in rough single fractures, which further proved the existence and importance of the inertial dependence of the non-Darcy coefficient.

Hydraulic aperture ( e h ) which is independent of inertial effect has been widely used as an important parameter to study the geometric properties of single fractures in previous studies. However, the inertial effect dependence of hydraulic aperture was gradually revealed in recent studies. Therefore, hydraulic aperture could provide a new method to characterize the non‐Darcy effect. In this study, 2D numerical simulations based on Navier‐Stokes equation were carried out to explore the feasibility of characterizing non‐Darcy effect based on hydraulic aperture. The results showed that the onset of flow regimes transition from Darcy to non‐Darcy could be well characterized by the hydraulic aperture. Before the non‐Darcy flow appeared, the initial hydraulic aperture ( e h0 ) maintained constant regardless of the inertial effect (quantified by non‐Darcy effect factor E ), once the non‐Darcy effect emerged, the hydraulic aperture decreased with the increase of inertial effect. The increased volume of low‐velocity zone in the rough single fractures was considered as the mechanism of non‐Darcy effect. And the existence of critical relative roughness ( R c ) was found and it would affect the critical non‐Darcy effect factor ( E d ) or critical Reynolds number ( Re c ) by affecting the development of low velocity region and preferential flow paths. From the overall trend, the hydraulic aperture decreased with the increase of E , and it fluctuated due to the development of recirculation zones (RZs). Moreover, the predictive model of non‐Darcy effect factor ( E ) based on hydraulic aperture was obtained and the validity of the model was preliminarily verified.