Discrete Fracture Network Research Articles

Fracture modeling of carbonate rocks via radial basis interpolation and discrete fracture network

Fracture modeling of carbonate rocks via radial basis interpolation and discrete fracture network

Effects of in-situ stress on heat transfer in fracture networks

Stress-induced fracture deformation is the principal cause for permeability change in geothermal systems. This study focuses on the influence of the nonlinear deformation and dilation effect of fractures on the geothermal system under the action of in-situ stress. By adopting a nonlinear constitutive model of rock fractures and embedding discrete fracture networks, numerical studies are first conducted to investigate the effects of different in-situ stress schemes on fracture aperture evolution using a rigid-body spring method. Based on the anisotropic aperture field of the fracture network caused by the in-situ stress, a finite element method is then used to study the flow and heat transfer process. The effects of different stress schemes on the heat flow transfer process are analyzed. Numerical simulation results show that when the ratio of horizontal to vertical stresses is not sufficient to cause shear dilation effects, the nonlinear normal deformation is the main factor affecting flow and heat transfer. In this case, the heat extraction efficiency is reduced. As the stress ratio increases, the shear dilation gradually becomes the dominant mechanism, and the heat extraction performance is improved. The obtained results provide a practical guide for geothermal site siting and optimizing heat extraction efficiency in geothermal reservoirs.

Fracture Network Analysis of Karstified Subis Limestone Build-Up in Niah, Sarawak

Understanding complex carbonate fracture networks and karstification at various geological scales is challenging, especially with limited multi-scale datasets. This paper aims to reduce uncertainty in the fracture architecture of Central Luconia karstified reservoirs by narrowing observational gaps between seismic and well data by using the discrete fracture models of exposed limestone outcrops as analogues for the subsurface carbonate reservoir. An outcrop-based fracture network characterisation of a near-surface paleo-karst at Subis Limestone combined with lineament analysis was conducted to extract fracture parameters. The karst structure was first delineated using a digital elevation map and outcrop examination. Then, topology analysis was performed, following the creation of two-dimensional discrete fracture models. Two main fracture sets oriented northeast–southwest and northwest–southeast and 79 potential dolines were identified. Fracture intersections, northeast–southwest major orientations, and drainage systems highly influenced the karst features. The Subis Limestone fracture model revealed that the highest number of fractures and total length of fractures were concentrated in the northern part of the Subis Limestone build-up (X: 250–350, Y: 150–250) and became denser towards the northwest direction of the outcrop (X: 600–800). The fractures in the Subis paleo-karsts appear isolated, with I-nodes ranging from 0.74 to 0.94. Hence, it is crucial to incorporate matrix porosity into multiple scales of fracture network modelling to improve upscaling and the modelling of fracture–vug networks, as well as to minimise the underestimation of discrete fracture networks in fractured and karstified limestone.

A Numerical Method for Evaluating the Collapse of High-Steep Scarp Slopes Based on the Bonded Block Model–Discrete Fracture Network Model

Geotechnical engineering works in deep-incised valleys or open-pit mining areas often encounter high-steep scarp slopes with a slope angle greater than 75°. This type of slope directly threatens the safety of construction personnel, so assessing their stability is essential to ensure construction safety. The natural geometry of high-steep scarp slopes possesses complexity in terms of geometric morphology, structural features of rock mass, and occurrence mechanisms of collapse. There is little research and less emphasis on the evaluation of the collapse risk of high-steep scarp slopes. In particular, the fracture of intact rock or rock bridges is generally ignored in the analysis of collapse processes. A bonded block model (BBM)–discrete fracture network (DFN) coupling characterization model for the high-steep scarp slope is proposed based on a high-steep scarp slope containing dominant joint sets on the left bank of the dam site of the Huangzangsi Water Conservancy Project (Qinghai Province, China). By using the model, the complex geometric forms of the surface of the high-steep scarp slope are quantified, and the fracture process of falling rock masses as well as the controlling effect of dominant joints on the collapse of the scarp slope are revealed. A strength reduction method based on the BBM–DFN model is constructed, and the safety factor of the collapse-prone scarp slope is evaluated. The research results show that (1) the BBM–DFN model can be used to describe the local collapse process; (2) the occurrence of dominant joints plays an important part in controlling the collapse process; (3) there are differences in the safety factor of the scarp slope with different coupling methods; the collapse and failure modes also differ. For safety considerations, the safety factor of the scarp slope on the left bank of the dam site area is determined to be 1.85. The research findings can be used to guide the safety assessment of high-steep scarp slopes and the formulation of both collapse risk prevention and control measures to ensure construction safety in high-steep scarp slope areas.

Effects of Topological Properties with Local Variable Apertures on Solute Transport through Three-Dimensional Discrete Fracture Networks

In this study, the effects of topological properties with local variable apertures on fluid flow and solute transport through three-dimensional (3D) discrete fracture networks (DFNs) were investigated. A series of 3D DFNs with different fracture density, length, and aperture distribution were generated. The fluid flow and solute transport through the models were simulated by combining the MATLAB code and COMSOL Multiphysics. The effects of network topology and aperture heterogeneity on fluid flow and transport process were analyzed. The results show that the fluid flow and solute transport exhibit a strong channeling effect even in the DFNs with identical aperture, in which the areas of fast and slow migration fit well with the high- and low-flow regions, respectively. More obvious preferential paths of flow and migration are observed in individual fractures for the DFN with heterogeneous aperture than the model with identical aperture. Increasing the fracture length exponent reduces the available flow and transport paths for sparse fracture networks but does not significantly change the flow and transport channels for dense fracture networks. The breakthrough curves (BTCs) shift towards the right and slightly lag behind as the fracture density decreases or the aperture heterogeneity increases. The advection–diffusion equation (ADE) model cannot properly capture the evolution of BTCs for 3D DFNs, especially the long tails of BTCs. Compared to the ADE model, the mobile-immobile model (MIM) model separating the liquid phase into flowing and stagnate regions is proven to better fit the BTCs of 3D DFNs with heterogeneous aperture.

Study on the influence of statistical geometrical characteristics on the permeability tensor of fractured rock masses by the equivalent pipe network method

The permeability of the host rock of the nuclear waste underground repository is one of the most keys for estimating the potential risk for nuclide waste to escape. Therefore, a three-dimensional discrete fracture network (DFN) model using the Monte Carlo method is generated based on the statistical geometrical characteristics of fracture network in Beishan I area which is selected as the site for constructing the nuclear waste repository in China. And the permeability of Beishan I fractured rock masses is calculated through the equivalent pipe network (EPN) method and representative elementary volume (REV) of it is determined. A rotating algorithm is developed to calculate different direction permeability. Furthermore, the influence of fracture statistical geometrical characteristics on REV size is explored by creating a series of DFN models with varying log-normal distribution of side lengths. The results show that the REV size of the Beishan Ⅰ area is 1.66 × 105 m3 and the corresponding principal permeability direction is 111.29°/75.36°. The REV size and permeability both increase with the increment of the mean and variance of the log-normal distribution increases. And with the increase of mean and variance, the growth rate gradually accelerated.

Heat Production Capacity Simulation and Parameter Sensitivity Analysis in the Process of Thermal Reservoir Development

The development of a geothermal system involves changes in the temperature field (T), seepage field (H), stress field (M), and chemical field (C) and the influence among them and injecting the heat extraction working fluid into the injection well that flows (migrating) through natural fractures and exchanges heat with the geothermal high-temperature rock. At the same time, the injection of low-temperature working fluid will induce thermal stress, resulting in changes in the reservoir temperature field and stress field. To study the influence factors and influence degree of heat production performance and mining life under multi-field coupling in the process of thermal reservoir development, based on THMC multi-field coupling numerical simulation software, this paper deeply studies the control differential equations and boundary coupling conditions of rock mass (fracture) deformation, seepage, heat exchange, the chemical reaction, and other processes based on the numerical solution method of the discrete fracture network model, simulating heat production capacity during the deep geothermal resource extraction process. The reservoir geological model analysis and generalization, parameter setting, boundary conditions, initial condition settings, mesh generation, and other steps were carried out in turn. Two different heat extraction working fluids, water, and CO2 were selected for numerical simulation in the mining process. The changes in the thermal reservoir temperature, net heat extraction rate, and SiO2 concentration during the thirty years of systematic mining were compared. The results show that CO2 has a better heat extraction effect. Finally, the reservoir thermal conductivity, heat capacity, well spacing, injection temperature, fracture spacing, fracture permeability, fracture number, fracture length, and other parameters were set, respectively. The parameter variation range was set, and the parameter sensitivity analysis was carried out. The numerical simulation results show that the engineering production conditions (injection temperature, well spacing) have little effect on the thermal efficiency and mining life, and the properties of fractures (fracture permeability, fracture number, fracture length) have a great influence.

Extension of discontinuous deformation analysis method for coupled reactive transport problems

Fluid flow and the accompanying reactant transport in fractured rock masses are key concerns for understanding the performance and safety assessments of subsurface engineering problems, particularly in waste disposal and tunneling project. Nevertheless, predicting solute transport behavior remains a great challenge due to computational sophistication and numerical instability if numerous discontinuities and chemical reaction is taken into consideration. In this study, the partial differential equation describing the reactive transport of solute is solved and coupled with the Discontinuous Deformation Analysis (DDA) method. First, the extended DDA method is validated by simulating an indoor dissolution test on the single rock fracture. Then, it is applied to establish a discrete fracture network (DFN), and be further employed to investigate the reactive transport behavior and geochemical evolution in a case of acidic mine drainage in a fractured stratum. The study analyzes the effects of stress field, reaction speed, and diffusion coefficient on the reactive transport. Numerical results indicate that the conclusions reached are consistent with current findings, combined with its prominent features in simulating the mechanical behavior of discontinuities, demonstrating that the extended method is a promising tool for the coupled hydro-mechanical-chemo analysis.

An extended colloid filtration theory for modeling Escherichia coli transport in 3-D fracture networks

Microbial transport in fractured carbonate rock using enhanced solutions is a significant and neglected research topic in the literature. We propose an extended colloid filtration theory (CFT) combined with a particle-tracking following streamlines (PTFS) model for the rapid prediction of breakthrough curves (BTCs) and plumes of pathogens in three-dimensional (3-D) discrete fracture networks (DFNs). We adapted CFT in porous media to pathogen transport in fractures containing Terra Rossa (soil) deposits. As an example of the model capability, a simulation was used to predict the 3-D motion field and Escherichia coli count in groundwater originating from the Forcatella managed aquifer recharge (MAR) Facility (Brindisi, Italy) using a DFN composed of 3,900 fractures. In arid regions, MAR facilities are significant for sustaining basic human needs, such as freshwater supply for drinking and crop production. The Markov chain Monte Carlo (MCMC) technique was applied to E. coli counts in the collected water samples to increase data representativeness. The pathogen transport coefficients were further supported by batch filtration tests carried out in the CNR/IRSA Laboratory (Bari, Italy). The mean E. coli attachment rate coefficient of 0.15 × 10−8 m2 d–1 (sticking efficiency = 1.1 × 10−8 m) resulted in a 2.1 log10 removal in 600 m of reclaimed water filtration. The simulation output visualized the E. coli 3-D pathways in groundwater and the positions of contaminated groundwater spring outflows on Forcatella Beach. The simulation results agreed with the mean MCMC output of E. coli concentrations in bathing water under unperturbed geochemical and environmental flow and transport conditions. However, results indicate that concentrations of pathogenic strains, parasites, and enteric viruses may enter the marine environment of MAR sites during flood periods.

Impact of artificial topological changes on flow and transport through fractured media due to mesh resolution

We performed a set of numerical simulations to characterize the interplay of fracture network topology, upscaling, and mesh refinement on flow and transport properties in fractured porous media. We generated a set of generic three-dimensional discrete fracture networks at various densities, where the radii of the fractures were sampled from a truncated power-law distribution, and whose parameters were loosely based on field site characterizations. We also considered five network densities, which were defined using a dimensionless version of density based on percolation theory. Once the networks were generated, we upscaled them into a single continuum model using the upscaled discrete fracture matrix model presented by Sweeney et al. (2019). We considered steady, isothermal pressure-driven flow through each domain and then simulated conservative, decaying, and adsorbing tracers using a pulse injection into the domain. For each simulation, we calculated the effective permeability and solute breakthrough curves as quantities of interest to compare between network realizations. We found that selecting a mesh resolution such that the global topology of the upscaled mesh matches the fracture network is essential. If the upscaled mesh has a connected pathway of fracture (higher permeability) cells but the fracture network does not, then the estimates for effective permeability and solute breakthrough will be incorrect. False connections cannot be eliminated entirely, but they can be managed by choosing appropriate mesh resolution and refinement for a given network. Adopting octree meshing to obtain sufficient levels of refinement leads to fewer computational cells (up to a 90% reduction in overall cell count) when compared to using a uniform resolution grid and can result in a more accurate continuum representation of the true fracture network.

Discontinuity development patterns and the challenges for 3D discrete fracture network modeling on complicated exposed rock surfaces

Natural slopes usually display complicated exposed rock surfaces that are characterized by complex and substantial terrain undulation and ubiquitous undesirable phenomena such as vegetation cover and rockfalls. This study presents a systematic outcrop research of fracture pattern variations in a complicated rock slope, and the qualitative and quantitative study of the complex phenomena impact on three-dimensional (3D) discrete fracture network (DFN) modeling. As the studies of the outcrop fracture pattern have been so far focused on local variations, thus, we put forward a statistical analysis of global variations. The entire outcrop is partitioned into several subzones, and the subzone-scale variability of fracture geometric properties is analyzed (including the orientation, the density, and the trace length). The results reveal significant variations in fracture characteristics (such as the concentrative degree, the average orientation, the density, and the trace length) among different subzones. Moreover, the density of fracture sets, which is approximately parallel to the slope surface, exhibits a notably higher value compared to other fracture sets across all subzones. To improve the accuracy of the DFN modeling, the effects of three common phenomena resulting from vegetation and rockfalls are qualitatively analyzed and the corresponding quantitative data processing solutions are proposed. Subsequently, the 3D fracture geometric parameters are determined for different areas of the high-steep rock slope in terms of the subzone dimensions. The results show significant variations in the same set of 3D fracture parameters across different regions with density differing by up to tenfold and mean trace length exhibiting differences of 3–4 times. The study results present precise geological structural information, improve modeling accuracy, and provide practical solutions for addressing complex outcrop issues.

High‐Resolution Characterization of Excavation‐Induced Fracture Network Using Continuous and Discrete Inversion Schemes

AbstractThe Meuse/Haute‐Marne Underground Research Laboratory provides the location for the experiment designed to investigate the induced fracture network around open or sealed galleries and drifts. The objective of this study is to investigate and reconstruct the hydraulic properties and the geometry of the induced fracture network to improve the insights and validate the conceptual model of the induced fracture network due to the stress redistribution during tunnel excavations. Within the presented study, the cross‐hole responses of the pneumatic tests were analyzed in the first step with an equivalent porous media—3D travel time‐based tomographic approach. In the next step selected 2D profiles of the 3D model domain were inverted using a discrete fracture network inversion approach. The database of the tomographic analysis is based on 18 gas injection tests and 151 pressure interferences, which were recorded between nine closely spaced boreholes. The travel time‐based inversion approach allowed for the reconstruction of the 3D gas diffusivity distribution between nine boreholes with high‐resolution. The applied discrete fracture network inversion approach is based on a transdimensional Markov Chain Monte Carlo (MCMC) methodology and it operates with reversible model updates (jumps) that change the problem dimensions, that is, the number, length and position of fractures within the model domain after each iteration step. The synthesis of the results between the reconstructed 3D diffusivity tomogram and the 2D fracture tomograms improved the insights into the spatial geometry of the induced fracture network around galleries.

Existence analysis of hydraulic conductivity representative elementary volume in fractured rocks based on three-dimensional discrete fracture network method

The complexity of fractures in rocks results in heterogeneity and anisotropy of their hydraulic properties. The existence of hydraulic conductivity representative elementary volume (KREV) is a fundamental question in comprehending the hydraulic behavior of fractured rock masses. This article introduces an analytical procedure for determining the existence of KREV in fractured rock masses, which considers the distribution of underground apertures using a discrete fracture network (DFN) model. The procedure follows four basic steps: inverse modeling to generate a 3D DFN; calibration of the effective aperture by combining the results of packer tests and flow simulation using finite difference methods; calculation of the equivalent conductivities in 10 domains and 31 different flow directions; analysis of the fit error between the conductivity ellipsoid to determine the existence of KREV. The proposed approach is applied to the fractured rocks surrounding the 3013 exploration adit near the underground powerhouse of the Three Gorges Project as a case study. The results show that most fractures in the study area have apertures ranging from 0.06 to 0.25 mm and follow a lognormal distribution. The KREV size is 14–16 times the fracture spacing. Investigation of this model has significant implications for understanding flow through fracture networks in related fields.

Open fracture clustering: Integrating subsurface and outcrop analogues, Asmari Formation, SW Iran

Fracture networks play a significant role in controlling fluid flow in carbonate reservoirs. Here we have utilized a modern and rigorous technique known as normalized correlation count (NCC) to investigate the 1D spatial arrangement and size distribution of conductive opening-mode fractures in the Asmari Formation, a naturally fractured carbonate reservoir in the Zagros Fold and Thrust Belt. To overcome subsurface sampling limitations, here we integrate data from a subsurface fault-related fold called the Gachsaran anticline, and from adjacent outcrops. The Oligo-Miocene Asmari Formation is composed of three units with different lithology, thickness, and stratal geometry. Subsurface and outcrop data from the Asmari Formation have different patterns of spatial arrangements. While the fractures from outcrop anticlines show indistinguishable from random arrangements, the subsurface scanlines from the Gachsaran anticline have various clustering patterns, from self-organized regularly-spaced fracture clustering in the backlimb to externally-imposed clustering in the forelimb. The NW-SE striking fracture set in outcrops and subsurface show clustering. Aperture-size data for all fracture sets are best-fit by log-normal distributions, but sampling biases against the smallest apertures may obscure a power-law aperture-size distribution. Forcing the power-law fit to aperture-size distribution indicates steeper slopes in the backlimb and gentler slopes in the forelimb, which reflect lower and higher strain, respectively. This is consistent with higher tendency of clustering in the forelimb of the anticline which correlates with higher productivity index. Our results demonstrate the ability of NCC to identify hierarchical fracture clusters as well as to distinguish between self-organized and extrinsic organization. The results of quantitative analysis of fracture clustering and patterns in the Asmari Formation provides valuable insights that can aid in the refinement of Discrete Fracture Network (DFN) models and help make informed decisions regarding fractured reservoir development and production.

A Comprehensive Method for Similarity Evaluation in Discrete Fracture Network Modeling of Jointed Rock Masses

A Comprehensive Method for Similarity Evaluation in Discrete Fracture Network Modeling of Jointed Rock Masses

Fractured Zone Detection Using Petrophysical Logs by Rescaled Windowed R/S Analyses and Discrete Fracture Network (DFN) of Yanchang Formation in Mahuangshan-Hongliugou Area, Ordos Basin, China

Structural fractures generally develop in the upper crust strata and are usually distributed in a convergent pattern, forming structural fracture zones with a specific strike. Fracture zones control the reservoir seepage system and seriously affect the migration and accumulation of oil and gas in tight sandstone reservoirs. Therefore, characterizing the characteristics of the fracture zones for tight oil exploration and development effectively is essential. In this paper, the variable scale fractal method is introduced to calculate the petrophysical log, and a new curve H is built. An intensity log is to characterize the intensity of structural fracture development. The H curve is in a good linear relationship with the intensity curve after the comparison of the H curve and intensity curve in 32 wells. A quantitative relationship between H and the intensity curve is established. Based on the parameters obtained from the core and image logs, the discrete fracture network model was established using H curves from more than 300 wells, and the structural fracture zone was analyzed. The model shows that the fracture zones formed by structural fractures are in S-N and NW-SE directions in the study area. The orientation of the structural fracture zone is consistent with that of the fractured fault zone and fault, and the development of the fractured zone is consistent with the regional tectonic evolution characteristics. The characteristics of the fracture zone explain the distribution law of oil accumulation and groundwater salinity in the study area.

Coupled fluid flow, solute transport and dissolution processes in discrete fracture networks: An advanced Discontinuous Galerkin model

Modeling dissolution processes in discrete fracture networks (DFNs) is a challenging task. Challenges are related to the highly nonlinear coupling between flow, mass transport, and reactive processes associated with fracture aperture evolution by dissolution. Further, advection-dominated transport due to fast fluid flow in fractures renders the problem more complex from a computational point of view, as traditional numerical methods may introduce unphysical oscillations or excessive numerical diffusion. The Discontinuous Galerkin (DG) method is known to be suitable for the simulation of advection-dominated transport. In this work, an advanced DG model is developed to model transport with dissolution in DFNs. We propose an upwind formulation to deal with the upstream concentration at the intersection of several fractures. The upstream concentration at an intersection node is calculated based on the average nodal concentrations of all the fractures having an inflow at that node, weighted by the volumetric fluxes of these fractures. The dispersion term is discretized with the Mixed Finite Element (MFE) method, which ensures the continuity of the dispersive flux at the intersection of fractures with different apertures. The obtained nonlinear coupled flow-transport-dissolution equations are discretized in time with a high-order scheme via the method of lines (MOL). Numerical examples and comparisons with standard finite element (FE) and finite volume (FV) solutions are performed to investigate the correctness and efficiency of the developed model. Results show that the new DG-DFN model avoids unphysical oscillations encountered with the standard FE method and strongly reduces the numerical diffusion observed with the upwind FV scheme. The DG-DFN model is then used to investigate the effect of the dissolution rate on the flow, transport, and aperture evolution processes for a single fracture and for a DFN. A quasi-linear evolution of the fracture aperture is observed for low dissolution rates. For high dissolution rates, a funnel-shaped enlargement is observed with a significant widening for the fractures near the inlet and minor effects for those away from the injection location.

Borehole stability in naturally fractured rocks with drilling mud intrusion and associated fracture strength weakening: A coupled DFN-DEM approach

Borehole instability in naturally fractured rocks poses significant challenges to drilling. Drilling mud invades the surrounding formations through natural fractures under the difference between the wellbore pressure (Pw) and pore pressure (Pp) during drilling, which may cause wellbore instability. However, the weakening of fracture strength due to mud intrusion is not considered in most existing borehole stability analyses, which may yield significant errors and misleading predictions. In addition, only limited factors were analyzed, and the fracture distribution was oversimplified. In this paper, the impacts of mud intrusion and associated fracture strength weakening on borehole stability in fractured rocks under both isotropic and anisotropic stress states are investigated using a coupled DEM (distinct element method) and DFN (discrete fracture network) method. It provides estimates of the effect of fracture strength weakening, wellbore pressure, in situ stresses, and sealing efficiency on borehole stability. The results show that mud intrusion and weakening of fracture strength can damage the borehole. This is demonstrated by the large displacement around the borehole, shear displacement on natural fractures, and the generation of fracture at shear limit. Mud intrusion reduces the shear strength of the fracture surface and leads to shear failure, which explains that the increase in mud weight may worsen borehole stability during overbalanced drilling in fractured formations. A higher in situ stress anisotropy exerts a significant influence on the mechanism of shear failure distribution around the wellbore. Moreover, the effect of sealing natural fractures on maintaining borehole stability is verified in this study, and the increase in sealing efficiency reduces the radial invasion distance of drilling mud. This study provides a directly quantitative prediction method of borehole instability in naturally fractured formations, which can consider the discrete fracture network, mud intrusion, and associated weakening of fracture strength. The information provided by the numerical approach (e.g. displacement around the borehole, shear displacement on fracture, and fracture at shear limit) is helpful for managing wellbore stability and designing wellbore-strengthening operations.

Investigating the thermo-hydro-mechanical evolution of a UK geological disposal facility due to disposal of high-heat-generating wastes

The heat released from the disposal of High-Heat-Generating Waste (HHGW) in a Geological Disposal Facility (GDF) will result in an increase in temperature, and therefore thermal expansion of both the engineering materials and surrounding rocks. The changes to the stress state arising from the thermal expansion of these materials, restrained by the surrounding rock mass, could affect several of the Thermal, Hydraulic, and Mechanical (THM) processes that operate after closure of a GDF. In this paper a methodology is described which suggests how the influence of thermal stress on fracture aperture around a GDF can be evaluated. A set of calculations on an illustrative Discrete Fracture Network (DFN) are reported. These calculations estimate the responses of the fractures to changes in rock stress, considering changes in fracture aperture and the resultant changes to effective intrinsic permeability.

Seepage stability analysis of a deep-buried tunnel in fractured rocks based on a non-Darcy hydro-mechanical coupled method

Non-Darcy seepage commonly occurs in fractured rocks with an adequate groundwater system around deep-buried tunnels. This study develops a numerical method for capturing the non-linear hydro-mechanical coupled processes in a rock mass system comprised of a tunnel lining, grouting circle, and surrounding rocks with discrete fracture networks. The proposed model is validated by comparing it with analytical solutions and field data obtained from a diversion tunnel of the Jinping II Hydropower Station on the Yalong River. The investigation of the non-Darcy flow regions around the diversion tunnel indicates that fractures with relatively higher conductivities dominate the non-Darcy effect, and the highest Forchheimer number reaches 0.5 during the period of torrential rain. The seepage stability of the diversion tunnel lining under three different cases (failure of the seepage control capacity of the grouting circle, potential karst pipeline penetrating the diversion tunnel, and torrential rain attacking the area) is further studied. It is found that the pore pressure at the external edge of the lining increases with the enhancement of the grouting circle permeability, and reaches 2,472.4 kPa when the karst pipeline is connected. Additionally, it is observed that the pore pressures obtained based on the Darcian and non-Darcian assumptions differ by up to 468.8 kPa in Case I. However, the non-Darcy effect gradually decreases as the lining permeability decreases. The results of this study provide practically useful suggestions for the design and operation of diversion tunnels subjected to high water pressure.