Understanding and investigating the gas transport mechanisms in shale fractures is crucial for ensuring the safe operation of various underground engineering activities. In this study, experimental investigations were conducted on shale fractures with different matching degrees to examine their deformation and seepage characteristics under different confining pressures and temperatures. The results demonstrate that the permeability coefficient continuously decreases with increasing confining pressure and temperature, while the GBP shows a continuous increase. Notably, both parameters exhibit a higher rate of change during the initial stage of pressure application, which gradually slows down thereafter. Significant differences were observed in the seepage characteristics between fractures with different matching degrees. The flow rate in poorly matched fractures (S1) was approximately three times higher than that in well-matched fractures (S2), while the GBP of S1 was only about half of S2. This indicates that the matching degree not only affects the initial aperture of fractures but also influences their deformation behavior. Specifically, fractures with lower matching degrees (S1) exhibited larger initial apertures compared to well-matched fractures (S2). Under confining pressure, the primary deformation mechanism in poorly matched fractures was shear failure, whereas well-matched fractures predominantly underwent plastic deformation.
Characterizing fracture networks is critical for groundwater development, geothermal exploitation, hydrocarbon recovery, and geological CO2 sequestration, yet their complex and uncertain spatial distribution poses a persistent challenge. This study proposes an intelligent inversion framework that integrates a 3D-UNet surrogate model, reversible jump Markov Chain Monte Carlo (rjMCMC), and multi-source data fusion for threedimensional discrete fracture network (DFN) characterization at the field scale. Within this framework, a 3D-UNet model trained on large datasets of fracture configurations, hydraulic head, and electrical potential provides an efficient initial inversion of fracture parameters. Fracture geometries are then extracted with the RANSAC algorithm and iteratively refined via rjMCMC, where the surrogate 3D-UNet replaces conventional forward modeling. This innovation reduces computational costs by an order of magnitude, enabling efficient large-scale inversion. Furthermore, the fusion of electrical potential with hydraulic head data enhances inversion accuracy by about 10 %. Validation demonstrates that the framework reliably reconstructs the spatial distribution of fracture networks, capturing both low-density zones and dominant hydraulic pathways in highly heterogeneous domains. By combining computational efficiency with improved accuracy, this approach offers a practical and scalable solution for field-scale fracture network characterization in a wide range of hydro-geological and engineering applications.
The shear mechanical behavior of joint surfaces is crucial for the safety and stability of rock engineering. Under constant normal load (CNL) conditions, the boundary element method is employed to calculate the distribution of normal contact forces. By combining the “dip angle” with the Patton model, a shear force calculation program analyzes the shear results for each discrete element. This approach allows for the investigation of microscopic phenomena during the shear process and provides insights into the effects of variations in load and shear displacement. Changes in load and shear displacement influence both the shear forces and the shear zones. The shear zones primarily expand and contract around their original regions. At local points of fracture shear, the patterns of shear force exhibit a differential effect, indicating that regions with higher shear stress are more prone to shear failure. This study provides a reference for understanding the variation of shear forces and shear zones during the shear process and offers a new perspective for investigating the mechanisms and patterns of shear failure.
Extensive research has been conducted on chemical dissolution and nonlinear fluid flow in single fractures. This study investigates the dynamic coupling between the chemical dissolution and nonlinear fluid flow for a single fracture in fractured carbonate rock. This was achieved with a numerical model that incorporates a mesh variation technique to capture the geometric evolution of the fracture driven by surface dissolution, thereby simulating the bidirectional feedback between morphological changes and flow behavior. The results revealed that the initial aperture heterogeneity and curvature jointly lead to nucleation-like dissolution patterns. As the simulation time progresses, synchronous decreases in the dissolution rate and surface roughness, along with a reduction in the nonlinear coefficient (B), and an increase in the critical Reynolds number (Rec) collectively confirm fracture smoothing and transition toward a hydraulically efficient and stable flow conduit. This process is governed by a self-reinforcing feedback mechanism between the localized flow and chemical dissolution, which ultimately controls the long-term evolution of the fracture.
The acidic dissolution behavior of calcite is an important process associated with geological environments such as karst-related systems and carbonate weathering. The study investigates the fundamental coupling mechanisms among geometric characteristics, fluid transport, and dissolution kinetics through controlled numerical simulations. The core innovation lies in the systematic quantification of the independent impact of the macroscopic curvature (represented by the aspect ratio Rm) on dissolution kinetics by introducing the parameter "Curvature-Dissolution Rate Coupling Response Coefficient (RC)" for the first time, and in revealing the interplay between two-dimensional fracture structure and fluid dynamics. Two-dimensional geometric models of isolated calcite particles and fractured matrices were established coupled with a dynamic mesh approach. For elliptical particles, the dissolution rate initially decreases and then increases with increasing Rm, a trend consistent with the variation of the specific surface area analogy value (Bv). The variation of RC indicates that the reaction rate is highly sensitive to curvature changes when Rm < 1. Multiple nonlinear regression analysis (Sd = 0.037 Va0.763ac1.822Rm-0.283 (R2 = 0.972)) further reveals that among the acid injection rate (Va), acid concentration (ac), and aspect ratio (Rm), ac exerts the most dominant control on the shrinkage degree (Sd). At the fracture scale, acid etching drives the morphology toward channelization and significantly attenuates the nonlinear behavior of fluid flow, clarifying the dynamic feedback mechanism inherent in fluid-fracture interaction.
Non-Darcy flow in rough fractures is strongly influenced by average aperture and heterogeneity, yet its quantitative characterization and prediction remain insufficiently understood. This study employs the Navier-Stokes equations to simulate flow in 30 fractures parameterized by average aperture (b) and aperture standard deviation (6, representing the degree of heterogeneity). A systematic investigation was conducted on the identification of streamline transitions and the characterization of the non-Darcy coefficient (/1) within rough fractures. Results reveal that flow behavior is non-monotonically controlled by b and 6: smaller 6 produces a more uniform aperture field, concentrating flow into a limited number of high-velocity preferential pathways and thereby amplifying inertial effects and nonlinear behavior, while larger 6 disrupts small-channel connectivity and redistributes flow across low-velocity zones, reducing /1. The transition from linear to nonlinear flow is jointly governed by b and 6, with larger apertures exhibiting lower critical Reynolds numbers for the onset of nonlinearity. A semi-empirical model was developed to predict /1 based on b and 6, capturing both the initial increase at low 6 and subsequent decline at high 6. Model validation demonstrates its effectiveness in quantifying the coupled influence of aperture size and heterogeneity on nonlinear flow in rough fractures. Furthermore, comparative analysis reveals that hydraulic parameters are significantly dependent on the observation scale. Parameters including K, p, c, and r exhibit quantitative drift when comparing the 40 mm scale in this study with a 20 mm scale dataset. This study contributes to more accurate predictions of hydraulic properties in fractured rock masses.
The real area of contact between rock surfaces plays a key role in determining the deformation properties of rock joints. There has been a growing trend toward measuring the real area of contact of rock joints using pressure-sensitive films, which are commonly used in the mechanical engineering industry as pressure sensors. A rock joint closure test was conducted, in which the real area of contact was measured using the pressure-sensitive film. The boundary element method (BEM), considering plastic deformation, which takes into account the measured topographies of the rock surfaces, was used to predict the real area of contact. The current investigation provides valuable insights into the proper usage of pressure-sensitive films in rock joint tests. The results show that, without considering destructive deformation, the elastic–plastic model effectively simulates the contact distribution. Under the studied loading conditions, the proportion of plastic deformation decreases as the load increases for the contact area. The negative exponential relationship between the load and the average aperture based on the Hurst index (H) is proposed. The contact stress can be calculated by the BEM according to the contact area comparison. This provides a valuable reference for rock fracture contact deformation theory.
To reveal the size-dependent seepage characteristics and pressure response mechanisms of fractured media under high-pressure environments, marble and limestone were selected as research objects with three sample scales designed. Three sets of systematic experiments-permeability coefficient tests under different confining pressures, hydraulic fracturing experiments, and seepage characteristic tests under constant flow rates were conducted to investigate the coupling seepage effects of lithology, size, and high pressure. The permeability coefficients of both rocks decrease exponentially with increasing confining pressure, and marble is more sensitive to confining pressure. The fracture initiation pressure of both lithologies increases linearly with sample diameter. Larger-sized samples have higher fracture initiation pressure, with large-size limestone forming a more stable diversion network and marble featuring faster fracture propagation. Marble exhibits stable linear seepage across all sizes, while limestone presents significant nonlinear exponential seepage, with the nonlinearity alleviated by increased size. This study identifies the size effect on nonlinear seepage and validates linear and exponential models for marble and limestone. Compared with existing studies, this study further clarifies the incremental understanding of multi-scale fractured rock seepage characteristics and provides a scientific basis for safety assessment in engineering projects such as water conservancy and hydropower development in Western Sichuan.
Coalbed methane (CBM) development is strongly controlled by pore structure evolution in deformed coals and its influence on hydraulic fracturing behavior. To clarify the multifractal characteristics of cross-scale pores and their control on fracturing effectiveness, this study investigated eight different deformation coals from the Ordos Basin using low-temperature CO2/N-2 adsorption (LT-CO(2)A/LT-N(2)A) and high-pressure mercury intrusion porosimetry (HMIP). Micropores (<2 nm), mesopores (2-50 nm), and macropores (>50 nm) were systematically characterized, and their pore size distributions (PSDs) were quantitatively analyzed using the Coal Structure Index (CSI) and multifractal theory. The results indicate that the multifractal parameters of macropores are significantly distinct from those of mesopores and micropores, exhibiting lower H (0.824-0.893) and D-1 (0.766-0.853), and higher alpha(0) (1.422-1.541), Delta D (1.230-1.408), and Delta alpha (1.459-1.642). Macropores controlled by tectonic deformation exhibit stronger heterogeneity compared to mesopores and micropores in local parts of the coal mass; PSD varies significantly with deformation rising, derived from the differential pore structure evolution during brittle-ductile transition and the multi-scale synergistic effects including maturity and composition. Combined with field fracturing curves, the results further indicate that the alpha(0), Delta D, and Delta alpha of macropores are negatively correlated with breakdown pressure, with correlation coefficients of 0.51, 0.61, and 0.59, respectively, and that strong local heterogeneity of macropores favors fracture initiation and propagation and reduces breakdown pressure. Cataclastic coal is the most favorable for hydraulic fracturing, followed by undeformed coal, whereas granulated coal shows the poorest fracturing performance.
The friction factor (f) plays a vital role in pressure head loss from microscopic fractures to large-scale faults. However, its dependence on tortuosity and inertial effects remains insufficiently understood. This study proposes a model that separates pressure loss into internal frictional and tortuosity-induced components. Three tortuosity plate models, sin-type, conical-type, and oval-type, were constructed with varying magnitudes and sizes. In addition, 500 artificial fractures with self-affine fractal characteristics were generated, grouped into four sets with identical aperture distributions. Results show that non-Darcy flow is most significant in oval-type fractures, followed by sin-type and least in conical-type. Increasing tortuosity magnitude induces earlier non-Darcy onset, while tortuosity size generates oscillatory fluctuations in flow resistance. The f-Re relationships, expressed by parameters m and n, were applied to characterize frictional behavior. In Darcy flow regimes with negligible inertial effects, n approaches-1, with m values highest in oval-type, intermediate in sin-type, and lowest in conical-type fractures, and positively correlated with tortuosity magnitudes and sizes. In non-Darcy regimes dominated by inertia, n exceeds-1, while the m-tortuosity relationships remain consistent with Darcy conditions, though m further exhibits dual dependence on tortuosity size, combining oscillatory variation with positive correlation. A linear relationship exists between internal frictional and tortuosity resistance when Re < 20, but deviates significantly at higher Re. These findings provide new insights into the coupled influence of tortuosity and inertial effects on head loss in interconnected rock fractures, with implications for mass and energy transport in fractured aquifers.
The closure behavior of rock fractures is a critical factor influencing the stability and seepage properties of rock masses, governed by multiple parameters such as fracture surface morphology and dip angle. While the significant influence of dip angle on the macroscopic mechanical behavior of fractures is recognized, the mesoscale evolution of contact pressure and its theoretical description for low to moderate dip angles (0°–30°) remain poorly understood. This study integrates 3D laser scanning and pressure-sensitive film experiments with theoretical analysis to systematically investigate the closure mechanisms and contact characteristics of granite fractures with various dip angles under normal stresses of 1–4 MPa. Results show that as the dip angle increases, the contact pressure distribution evolves from a symmetrical circular pattern to a strip-like pattern. The average contact stress decreases significantly (from 7 to 8.5 MPa to 1–1.8 MPa), while the contact ratio increases exponentially (from 6
The theory of rock joint closure is a crucial research subject in rock mechanics. However, existing theoretical studies predominantly analyze contact processes between two rough surfaces, which has proven effective for studying joint stress and contact deformation. Yet this approach fails to adequately utilize relevant parameters to analyze stress-induced contact behavior, as both surfaces undergo geometric deformation. To address this, the equivalent theory is employed, transforming the interaction between two rough surfaces into that between a rough composite surface and a rigid plane. This method effectively describes the evolution of the composite surface during contact. Combining composite surface parameters with the boundary element method (BEM), the composite surface enables accurate elastoplastic simulation of joint stress variations by using conjugate gradient method. The roughness parameters of the composite surface can be more conveniently expressed relative to the two surfaces. The multi-peak structure of the composite surface reveals contact nuclei during joint closure, along with independent development and mutual fusion processes, while stress distribution from the nucleus center outward follows specific attenuation patterns. The composite surface model allows precise calculation of displacement at spatial points. The composite surface proposed in this study holds significant theoretical value for investigating joint deformation mechanisms.
Accurately identifying preferential flow paths in fractured geological media is critical for groundwater hazard prevention and pollution remediation. However, existing identification methods are hindered by the complexity of fracture geometry and network topology, leading to low efficiency and limited accuracy. A rapid identification method for three-dimensional preferential flow paths is proposed, based on a flow resistance model with topological search. A flow resistance coefficient is defined by integrating fracture length, aperture, and roughness (JRC) via a modified cubic law and Barton’s model. The fracture network is represented as a weighted graph, and the A-STAR algorithm with a K-shortest paths strategy efficiently identifies the top K paths with minimal flow resistance (i.e., preferential flow paths). The method is validated using three numerical examples of increasing complexity. The results demonstrate that the minimum-resistance flow paths identified using the graph-theory-based method are consistent with the preferential flow paths obtained from conventional finite-element-based numerical simulations; lower flow resistance of fractures corresponds to higher flow rates and more pronounced dominant flow; high fracture densities yield numerous flow pathways and dispersed flow, making it difficult to form distinct preferential flow channels; larger fractures enhance connectivity and facilitate dominant channels. This identification approach offers an effective tool for analyzing channelization patterns and supporting engineering applications.
Nonlinear flow in rock fractures is often described by the Forchheimer equation; however, there are still flaws and contradictions in the understanding of the Forchheimer equation and nonlinear evolution mechanisms, such as differing interpretations of various hydraulic apertures and the factors affecting the linear and nonlinear coefficients. This study utilizes wavelet analysis filtering and reconstruction techniques to investigate the impact of roughness at different scales on nonlinear seepage. It further provides new insights into the Forchheimer equation and hydraulic aperture by introducing the concepts of the ideal curve, the ideal Reynolds number, and process hydraulic aperture, along with the proposed criteria for determining them. A new mathematical model, incorporating geometric characteristics and Reynolds number effects, is developed and validated, effectively describing hydraulic aperture variations during nonlinear flow. The findings clarify the evolution of hydraulic aperture and flow coefficients, resolve previous contradictions, and advance understanding of nonlinear seepage mechanisms.
This study proposed a non-invasive electrical resistivity method to track the fracture remediation process driven by microbially induced calcium carbonate precipitation. A seven-cycle grouting experiment was conducted on a rough-walled single fracture with an average aperture of 1.02 mm. Resistivity was measured with a four-electrode method, and water pressure and precipitate thickness provided independent validation. The results indicated a non-uniform resistivity distribution across the fracture, with the resistivity increment decreasing from 26.6 Ω·m to 2.4 Ω·m along the flow path. The calcium carbonate precipitation and water pressure also exhibited an analogous spatial distribution pattern. The overall negative relationship between the electrical resistivity and the fracture aperture was well fitted by the power-law function with an R² value greater than 0.99. Incorporating pressure data into the established resistivity-aperture empirical relation yielded a reduction in absolute aperture estimation error, decreasing it from 0.103 mm to 0.075 mm. This study provided a promising non-invasive approach for real-time monitoring of the fracture remediation process, with potential applications in evaluating sealing efficiency in deep underground engineering.
Accurate characterization of fractured media is fundamental in the geological and geotechnical engineering applications such as coal mine production, deep geological disposal and enhanced geothermal systems (EGS). However, traditional inversion strategies are limited in their ability to characterize high-dimensional and nonGaussian fractured media. Furthermore, a significant amount of observation well was employed during the inversion process in the previous studies. In this work, we proposed a joint inversion framework based on deep learning technique to overcome the limitations of the traditional strategies and the challenge of excessive use of observation wells. The convolutional variational autoencoder (CVAE) network was trained to parameterize the fractured media. After that, the ensemble smoother with multiple data assimilation (ESMDA) combined with the CVAE to characterize fractured media assimilating the hydraulic tomography (HT) and thermal tracer tomography (TT) data. A numerical study using four observation points validates the framework's reliability. The characterization errors for single-data cases are 16.9 % (HT) and 18.1 % (TT), decreasing to 16.7 % when both types of data are incorporated, demonstrating the synergies of multisource data. Sequentially, the framework is extended to the real-world scenario. The results show that our framework can effectively characterize the fractured media, capturing more features while addressing the challenge posed by excessive use of observation wells through the integration of multisource data. Our framework provides valuable insights into the characterization of fractured media in the practical engineering applications and highlights the benefits of multisource data assimilation.
Fluid flow around irregular cylinders that are widespread in nature and engineering deserves attention, and its hydrodynamic coefficients can be predicted by regular cylinders. For this purpose, the flow fields around seven regular cylinders and five irregular cylinders were analyzed via computational fluid dynamics simulations at low Reynolds numbers (90 ≤ Re ≤ 150). When the surrounding fluid flow is unstable, vortex shedding will occur behind both irregular cylinder and regular cylinder. No vortex shedding behind the cylinder with small Ar at low Reynolds number: the elliptical cylinder with Ar = 2/5 at Re ≤ 150; the elliptical cylinders with Ar = 3/5 at Re ≤ 120; Irregular cylinder 4 with Ar = 2/5 at Re ≤ 110. The drag coefficient curve of the irregular cylinder has a large peak and a small peak, and the lift coefficient curve of the irregular cylinder deviates from CL = 0. Cylinders with the same geometric quantization parameters (Ar, area, and perimeter) have similar hydrodynamic characteristics.
Rarefied gas flow in rock fractures governed by Knudsen flow is employed in residual gas capture and storage. Previous studies mainly adopted the classical theory under the continuum hypothesis to describe the gas flow in fractures, and the research on the interaction between molecules and walls of rarefied gas flow in heterogeneous fracture channels was limited. The molecular flow module is established based on the equations of the kinetic theory of gases, and the mathematical particle tracking module is established based on Newton's laws of motion. 2D and 3D fracture models with different tortuosity and aperture were established to quantify the influence of surface heterogeneity on the transmission probability of rarefied gas flow. The research findings indicate that as the aperture standard deviation and tortuosity increases, the transmission probability correspondingly decreases. Furthermore, in the case of a 3D channel, the aperture standard deviation exerts a dominant influence on transmission probability, and the fitted relationship F = a - b & sdot;sigma(b)<^>c has been derived. This discovery underscores the limitations inherent in 2D models: these models exhibit anomalous molecular retention phenomena, which hinder their ability to accurately represent true three-dimensional transport mechanisms. In contrast, 3D models feature dominant channels that mitigate the impact of local surface roughness peaks. The findings offer universally applicable theoretical tools for regulating gas flow in a wide range of deep geological engineering applications, contributing to the resolution of the dual challenges of "efficient resource extraction" and "environmental risk prevention and control."
Fractures with different multiscale roughnesses are widely encountered in subsurface engineering, and their impact on non-Darcy flow is critical for predicting fluid transport. Although numerous studies have examined fluid through multiscale fractures, the specific effects of wavelet decomposed roughness wavelengths (L) and amplitudes (Delta) remains poorly understood. Here, a decomposed roughness model was introduced using multiscale wavelet decomposition to construct fracture geometries with controlled L and Delta. Relative velocity (V-r) field was employed to emphasize decomposed roughness disturbance and classify it into enhanced and reduced velocity region (EVE/RVR). The disturbance intensity (gamma) used to quantify disturbance. Sensitive analysis indicates that L exhibits a noteworthy role on the deviation of Darcy flow. Hydraulic parameters e(h), beta, and Re-c scale with L via power-law relationships. Specifically, beta exhibits periodic oscillations as L increases, which was modeled by adding a sinusoidal term. It observed that gamma increases with Re due to stronger viscous and inertial effects. However, at Re = 100-1000, gamma decreases as EVR expands and RVR contracts. The gamma variation with L using a logistic oscillation modulation function. Inertial effects and roughness also lead to recirculation zones (RZs). The recirculation zone area ratio (theta) is as follows: theta = A(1)Re/(A(2) + Re), indicating transitions from scale to inertial dominated stage. The relationship between gamma and Re-c transitions from a inertial dominated regime at small scales to a scale dominated regime at larger scales. This work offers new insight into multiscale roughness disturbance on the initial deviation of Darcy flow.
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.