Flow In Discrete Fracture Networks Research Articles

Algorithmic approach to discrete fracture network flow modeling in consideration of realistic connections in large-scale fracture networks

Analyzing rock mass seepage using the discrete fracture network (DFN) flow model poses challenges when dealing with complex fracture networks. This paper presents a novel DFN flow model that incorporates the actual connections of large-scale fractures. Notably, this model efficiently manages over 20,000 fractures without necessitating adjustments to the DFN geometry. All geometric analyses, such as identifying connected fractures, dividing the two-dimensional domain into closed loops, triangulating arbitrary loops, and refining triangular elements, are fully automated. The analysis processes are comprehensively introduced, and core algorithms, along with their pseudo-codes, are outlined and explained to assist readers in their programming endeavors. The accuracy of geometric analyses is validated through topological graphs representing the connection relationships between fractures. In practical application, the proposed model is employed to assess the water-sealing effectiveness of an underground storage cavern project. The analysis results indicate that the existing design scheme can effectively prevent the stored oil from leaking in the presence of both dense and sparse fractures. Furthermore, following extensive modification and optimization, the scale and precision of model computation suggest that the proposed model and developed codes can meet the requirements of engineering applications.

Review of discrete fracture network characterization for geothermal energy extraction

Geothermal reservoirs are highly anisotropic and heterogeneous, and thus require a variety of structural geology, geomechanical, remote sensing, geophysical and hydraulic techniques to inform Discrete Fracture Network flow models. Following the Paris Agreement on reduction of carbon emissions, such reservoirs have received more attention and new techniques that support Discrete Fracture Network models were developed. A comprehensive review is therefore needed to merge innovative and traditional technical approaches into a coherent framework to enhance the extraction of geothermal energy from the deep subsurface. Traditionally, statistics extracted from structural scanlines and unmanned aerial vehicle surveys on analogues represent optimum ways to constrain the length of joints, bedding planes, and faults, thereby generating a model of the network of fractures. Combining borehole images with seismic attributes has also proven to be an excellent approach that supports the stochastic generation of Discrete Fracture Network models by detecting the orientation, density, and dominant trends of the fractures in the reservoirs. However, to move forward to flow modelling, computation of transmissivities from pumping tests, and the determination of hydraulically active fractures allow the computation of the hydraulic aperture in permeable sedimentary rocks. The latter parameter is fundamental to simulating flow in a network of discrete fractures. The mechanical aperture can also be estimated based on the characterization of geomechanical parameters (Poisson’s ratio, and Young’s modulus) in Hot Dry Rocks of igneous-metamorphic origin. Compared with previous review studies, this paper will be the first to describe all the geological and hydro-geophysical techniques that inform Discrete Fracture Network development in geothermal frameworks. We therefore envisage that this paper represents a useful and holistic guide for future projects on preparing DFN models.

The Criteria for Transition of Fluid to Nonlinear Flow for Fractured Rocks: The Role of Fracture Intersection and Aperture

Understanding the fluid pattern is of special significance for estimating the hydraulic conductivity of fractured rock masses. The nonlinearity of fluid flow in discrete fracture networks (DFNs) originates from inertial effects and is enhanced by complex geometric topologies, which produces additional viscous friction and is subject to inertia effects, consequently transitioning the fluid to the nonlinear flow regime. Therefore, it is important to obtain the critical conditions for the transition of a fluid from laminar to turbulent flow. To investigate the role of fracture aperture and fracture intersection on the onset of the transition of a fluid to nonlinear flow in fractured rocks, the fluid dynamic computation was performed by solving Navier–Stokes (N–S) equations in DFN models. The results show that the flow flux initially linearly correlates with the hydraulic gradient (J) and the permeability of DFNs initially remains constant. As the hydraulic gradient increases, the flow flux presents a strong nonlinear relationship with the hydraulic gradient, and the permeability decreases dramatically. In particular, significant inertial effects appear earlier with a large fracture aperture or a dense fracture intersection. A critical hydraulic gradient (Jc) is proposed to judge the onset of nonlinear flow. The mathematical expression of Jc and Forchheimer coefficients A and B involving the fracture aperture and fracture intersection density is established through a multiple regression algorithm. Finally, the reliability of the predictive model was verified by comparing the results of the prediction and fluid dynamic computation of a series of DFN models with well-known geometric distributions. The consistency of the fitted equations and a correlation coefficient greater than 0.9 between them indicate that the predictive model proposed in this study is reliable.

On the onset of nonlinear fluid flow transition in rock fracture network: Theoretical and computational fluid dynamic investigation

Fluid flow regimes affect the determination of hydraulic conductivity of fractured rocks, and the critical criteria for the onset of nonlinear fluid flow transitions in discrete fracture networks (DFNs) of rocks have yet to be established. First, the factors causing the fluid flow transition regime of fracture intersections and rough surface fractures are theoretically and numerically analyzed. This reveals that the fluid flow regime is governed by the fracture aperture, density of fracture intersections, surface roughness, and Reynolds number (Re). Then, these identified parameters are redefined in DFN models, and their influence on the onset of nonlinear fluid flow is further investigated by performing computational fluid dynamic analysis. The results show that the fracture intersection and aperture play a more significant role in the linear-to-nonlinear fluid flow transition than the fracture aperture heterogeneity. With the increase in the fracture aperture, unevenness of fracture surfaces, and connectivity of DFNs, the onset of the nonlinear fluid flow appeared at the lower flow velocity. With the Forchheimer equation, it is found that the critical hydraulic gradient Jc, defined as the hydraulic gradient at which inertial effects assume 10% of the total pressure loss, is highly correlated with the fracture aperture, fracture intersection, and roughness of the surface. Finally, the mathematical expression of Jc and the Forchheimer coefficients are formulated based on the regression analysis of fluid dynamic computation results, which provides an approach to determine whether the cubic law should be applied as governing equations for the computation of fluid flow in DFNs.

Efficient partitioning of conforming virtual element discretizations for large scale discrete fracture network flow parallel solvers

Discrete Fracture Network models are largely used for large scale geological flow simulations in fractured media. For these complex simulations, it is worth investigating suitable numerical methods and tools for efficient parallel solutions on High Performance Computing systems. In this paper we propose and compare different partitioning strategies, that result to be highly efficient and scalable, overperforming the classical mesh partitioning approach used to balance the workload of a conforming mesh among several processes. The proposed DFN-based partitioning strategies rely on the distribution of the fractures among parallel processes. The computational cost of the DFN-based partitionings is very small compared to the cost of classical mesh partitioning and the numerical results prove their effectiveness and good performances in solving linear systems for realistic DFN flow simulations.

Investigation of Nonlinear Flow in Discrete Fracture Networks Using an Improved Hydro-Mechanical Coupling Model

Fractures commonly exist in rock masses; the coalescence of fractures provides fluid flow pathways in a fractured rock mass and greatly increases the flow capacity of fractured rock. This work aims to study the characteristics of nonlinear flow in fractures. A series of tests were conducted and indicated that the Forchheimer law performed well when describing the nonlinear relationship between hydraulic gradient and flow. The test results also indicate that higher water pressure may induce stronger nonlinearity. Additionally, the linear and nonlinear coefficients of the Forchheimer law increase with a decrease in the particle size of the filling material in fractures. On the basis of the laboratory results, the classical Forchheimer law was modified by considering the influence of stress on the variation of fracture aperture. A hydro-mechanical coupling model for fractured rock masses was built and programmed with a subroutine through ABAQUS. Furthermore, a random discrete fracture network was generated and simulated to prove that a high flow velocity will result in a nonlinear flow, not only in a single fracture, but also in a fracture network. The numerical results from fractured rock masses show that a ratio of the flow to the hydraulc gradient will change the flow from linear to weak nonlinearity and, finally, to strong nonlinearity with an increase in the hydraulic gradient. It also shows that the linear and nonlinear coefficients increase with an increase in the confining pressure and that they decrease with an increase in the aperture. Due to the complexity of fracture channels, a nonlinear flow is likely to occur in a fractured rock mass. Finally, the developed model was applied to simulate the flow behavior of underground engineering; the results show that the smaller the hydraulic aperture is, the higher the water pressure is required to be in order to change the flow regime from linear to nonlinear.

Liquid Breakthrough Time in an Unsaturated Fracture Network

AbstractWe present an experimental and theoretical study of unsaturated flow in discrete fracture networks. The focus is on the breakthrough time of infiltrating liquid through the fracture networks as well as the spatial distribution of local flow status under a wide range of flow rates. Through visualized experiments, the fluid motion in the fracture network is recorded and analyzed, and gravity‐driven preferential pathways are observed. The location of the breakthrough path is found to be insensitive to the flow rate under the single‐inlet flow condition. Based on analogy with slug migration in a single fracture and considering liquid accumulation due to the capillary barrier effect at the fracture intersections, we propose a theoretical model of liquid breakthrough time through an initially dry fracture network. The theoretical predictions are found to be in excellent agreement with the experimental results. Using the model, we analyze the relative contribution of times needed for liquid migration and for local liquid accumulation to the breakthrough time. We also find that the spatial distribution of local flow status at the steady‐state varies significantly with the flow rate. The improved understanding and the proposed model are of significance for predicting the time required for a fluid/contaminant to reach a certain depth in unsaturated fractured media, and may provide new physical insights on the complex unsaturated flow dynamics.

Discontinuous boundary elements for steady-state fluid flow problems in discrete fracture networks

Modeling fluid flow in three-dimensional (3D) Discrete Fracture Networks (DFNs) is of relevance in many engineering applications, such as shale oil/gas production, geothermal energy extraction, nuclear waste disposal and CO2 sequestration. A new Boundary Element Method (BEM) technique with discontinuous quadratic elements, in conjunction with a parallel Domain Decomposition Method (DDM), is presented for the simulation of the steady-state fluid flow in DFNs, consisting of stochastically generated 3D planar fractures, arbitrarily oriented, and having differing hydraulic properties. Numerical examples characterized by DFNs of increasing complexity are proposed to show the accuracy and the efficiency of the presented technique, that provides good approximations of the fluid flow around domain interfaces, where the solution usually displays sharp gradients, like around intersections between traces (the segments originated by the intersection between two fractures), intersections between traces and fracture boundaries, or intersections between fractures and wellbores. The conjunction with a DDM approach is a promising strategy to speed up the computations, by also exploiting the advantages of parallel computing techniques. The technique is implemented in the code PyDFN3D, available at https://github.com/BinWang0213/PyDFN3D.

Modeling radial groundwater flow in fractured media using fracture continuum approach

Two modeling approaches are commonly utilized for simulating flow in fractured formations: the discrete fracture network (DFN) approach and the stochastic continuum (SC) approach. Although the DFN approach is the most accurate, it has computational and memory constraints. The SC approach ensures fast processing but results in system over-homogenization. The fracture continuum (FC) approach arises as an integrated technique that incorporates the merits of both approaches. The main objective of this research is to develop a computationally efficient technique based on the FC approach to simulate the radial groundwater flow towards wells through two-dimensional fractured media under both steady and transient conditions. A stochastic generation of the DFN is performed in a Monte Carlo framework taking into account wells positioning. The DFN flow system is solved by applying the mass balance equation at fracture intersections. Fracture segments are mapped onto grids of 1 × 1 m and 5 × 5 m resolution as conductivity and specific storage cells. The grid flow problem is solved via MODFLOW. Flow and head discrepancies between the proposed technique and the DFN approach (reference solution) are assessed in steady and transient conditions. A grid-conductivity correction is needed to preserve the DFN flow in the presence of wells. A porosity estimation is proposed to identify the grid-pressure transient response. Promising flow and head results are observed for fine and coarse grid models. Some of the studied cases show large discrepancies in the maximum drawdown obtained on the coarse grid model. Accordingly, a new technique is proposed to handle such discrepancies and is found efficient in transient simulations (e.g., 11% and 26.12% discrepancies are minimized to − 0.93% and − 1.03% for two studied cases). The adopted mapping technique is found efficient when the interest is to estimate the average drawdown over an aquifer as correlation coefficients of 0.89 and 0.97 are found for the coarse and fine grid models, respectively when compared to the DFN model. However, the technique has limitations in estimating the drawdown at locations of wells.

Investigating the effects of intersection flow localization in equivalent-continuum-based upscaling of flow in discrete fracture networks

Abstract. Predicting effective permeabilities of fractured rock masses is a crucial component of reservoir modeling. Its often realized with the discrete fracture network (DFN) method, whereby single-phase incompressible fluid flow is modeled in discrete representations of individual fractures in a network. Depending on the overall number of fractures, this can result in high computational costs. Equivalent continuum models (ECMs) provide an alternative approach by subdividing the fracture network into a grid of continuous medium cells, over which hydraulic properties are averaged for fluid flow simulations. While continuum methods have the advantage of lower computational costs and the possibility of including matrix properties, choosing the right cell size to discretize the fracture network into an ECM is crucial to provide accurate flow results and conserve anisotropic flow properties. Whereas several techniques exist to map a fracture network onto a grid of continuum cells, the complexity related to flow in fracture intersections is often ignored. Here, numerical simulations of Stokes flow in simple fracture intersections are utilized to analyze their effect on permeability. It is demonstrated that intersection lineaments oriented parallel to the principal direction of flow increase permeability in a process we term intersection flow localization (IFL). We propose a new method to generate ECMs that includes this effect with a directional pipe flow parameterization: the fracture-and-pipe model. Our approach is compared against an ECM method that does not take IFL into account by performing ECM-based upscaling with a massively parallelized Darcy flow solver capable of representing permeability anisotropy for individual grid cells. While IFL results in an increase in permeability at the local scale of the ECM cell (fracture scale), its effects on network-scale flow are minor. We investigated the effects of IFL for test cases with orthogonal fracture formations for various scales, fracture lengths, hydraulic apertures, and fracture densities. Only for global fracture porosities above 30 % does IFL start to increase the systems permeability. For lower fracture densities, the effects of IFL are smeared out in the upscaling process. However, we noticed a strong dependency of ECM-based upscaling on its grid resolution. Resolution tests suggests that, as long as the cell size is smaller than the minimal fracture length and larger than the maximal hydraulic aperture of the considered fracture network, the resulting effective permeabilities and anisotropies are resolution-independent. Within that range, ECMs are applicable to upscale flow in fracture networks.

A Reduced Basis Method for a PDE-constrained optimization formulation in Discrete Fracture Network flow simulations

In classic Reduced Basis (RB) framework, we propose a new technique for the offline greedy error analysis which relies on a residual-based a posteriori error estimator. This approach is as an alternative to classical a posteriori RB estimators, avoiding a discrete inf-sup lower bound estimate. We try to use less common ingredients of the RB framework to retrieve a better approximation of the RB error, such as the estimation of the distance between the continuous solution and the reduced one. In particular we focus on the application of the reduction model for the flow simulations in underground fractured media, in which high accurate simulations suffer for the complexity of the domain geometry. Finally, some numerical tests are assessed to confirm the viability and the efficacy of the technique proposed.

Correction to: Multilevel Graph Partitioning for Three-Dimensional Discrete Fracture Network Flow Simulations

The publication of this article unfortunately contained a mistake. The assignment to the affiliations of author Satish Karra was not correct

Multilevel Graph Partitioning for Three-Dimensional Discrete Fracture Network Flow Simulations

We present a topology-based method for mesh-partitioning in three-dimensional discrete fracture network (DFN) simulations that takes advantage of the intrinsic multi-level nature of a DFN. DFN models are used to simulate flow and transport through low-permeability fractured media in the subsurface by explicitly representing fractures as discrete entities. The governing equations for flow and transport are numerically integrated on computational meshes generated on the interconnected fracture networks. Modern high-fidelity DFN simulations require high-performance computing on multiple processors where performance and scalability depends partially on obtaining a high-quality partition of the mesh to balance work-loads and minimize communication across all processors. The discrete structure of a DFN naturally lends itself to various graph representations, which can be thought of as coarse-scale representations of the computational mesh. Using this concept, we develop two applications of the multilevel graph partitioning algorithm to partition the mesh of a DFN. In the first, we project a partition of the graph based on the DFN topology onto the mesh of the DFN and in the second, this DFN-based projection is used as the initial condition for further partitioning refinement of the mesh. We compare the performance of these methods with standard multi-level graph partitioning using graph-based metrics (cut, imbalance, partitioning time), computational-based metrics (FLOPS, iterations, solver time), and total run time. The DFN-based and the mesh-based partitioning methods are comparable in terms of the graph-based metrics, but the time required to obtain the partition is several orders of magnitude faster using the DFN-based partitions. The computation-based metrics show comparable performance between both methods so, in combination, the DFN-based partitions are several orders of magnitude faster than the mesh-based partition. Moreover, the method which uses the DFN-partition solution as the initial condition of the mesh partition provided cut and imbalance values that were close to the mesh-based partition but in a fraction of the time. In turn, this hybrid method outperformed both of the other methods in terms of the total run time.

Determination of fracture apertures via calibration of three-dimensional discrete-fracture-network models: application to Pahute Mesa, Nevada National Security Site, USA

The Pahute Mesa Corrective Action Unit on the Nevada National Security Site (USA) contains several fractured aquifers that can potentially provide high-permeability pathways for migration of radionuclides away from underground nuclear testing locations. Though geometric properties of fractures such as length and orientations, can generally be obtained by geophysical methods and borehole image analyses, their hydraulic properties (primarily influenced by their apertures) are often unknown or have a high degree of uncertainty. This study presents a novel approach to determine fracture apertures at a site by integrating numerical models of flow with data of certain geometrical fracture attributes and hydraulic response of the system because of long-term pumping. Discrete fracture network (DFN) models were used in this study to construct a three-dimensional (3-D) flow model for lava-flow aquifers of Western Pahute Mesa (WPM). Models built using a 3-D DFN code, dfnWorks, were calibrated to hydraulic drawdown observations recorded at the site of a forced-gradient experiment (BULLION FGE). Multiple conceptualizations of boundary conditions, fracture aperture distributions, and realizations of the fracture network were considered to simulate flow and migration of particles between an injection well and a pumping well. Calibrating the DFN flow models to hydraulic drawdown data constrained the ranges of aperture values and helped develop a realistic description of the properties of fractured rocks in WPM. The aperture values resulting from this study are expected to enhance understanding of radionuclide transport in WPM and support the development of large-scale flow and transport models.

Enhancing fracture-network characterization and discrete-fracture-network simulation with high-resolution surveys using unmanned aerial vehicles

A workflow is presented that integrates unmanned aerial vehicle (UAV) imagery with discrete fracture network (DFN) geometric characterization and quantification of fluid flow. The DFN analysis allows for reliable characterization and reproduction of the most relevant features of fracture networks, including: identification of orientation sets and their characteristics (mean orientation, dispersion, and prior probability); scale invariance in distributions of fracture length and spatial location/clustering; and the distribution of aperture values used to compute network-scale equivalent permeability. A two-dimensional DFN-generation approach honors field data by explicitly reproducing observed multi-scale fracture clustering using a multiplicative cascade process and power law distribution of fracture length. The influence of aperture on network-scale equivalent permeability is investigated using comparisons between a sublinear aperture-to-length relationship and constant aperture. To assess the applicability of the developed methodology, DFN flow simulations are calibrated to pumping test data. Results suggest that even at small scales, UAV surveys capture the essential geometrical properties required for fluid flow characterization. Both the constant and sublinear aperture scaling approaches provide good matches to the pumping test results with only minimal calibration, indicating that the reproduced networks sufficiently capture the geometric and connectivity properties characteristic of the granitic rocks at the study site. The sublinear aperture scaling case honors the directions of dominant fractures that play a critical role in connecting fracture clusters and provides a realistic representation of network permeability.

Uncertainty quantification analysis in discrete fracture network flow simulations

In the framework of underground flow simulations in fractured media, fractures may act as preferential paths, and may have a strong impact on the flow. The discrete fracture network (DFN) model allows for an explicit representation of the interconnected fractures. Flux on each fracture is assumed to be driven by Darcy Law, and suitable matching conditions are imposed along fracture intersections, ensuring flux balance and head continuity. The exact displacement of fractures in the network is usually not known, and networks are typically generated sampling geometrical and hydro-geological properties from probabilistic distributions; this stochastic generation is likely to generate geometrical configurations very challenging for the meshing process, which is a major issue in the framework of DFN simulations. Stochasticity in geometrical parameters may also result in a nonsmooth behavior of the quantity of interest with respect to the stochastic parameters. In order to quantify the influence of these stochastic parameters on the output of DFN models, we propose an approach based on the geometric Multi Level Monte Carlo (MLMC) method, applied in conjunction with a well assessed underlying solver for performing DFN flow simulations. Key properties of the solver are its capability of circumventing the need of conforming meshes, and its consequent extreme robustness with respect to geometrical complexities in the network. These features make the solver quite suitable to be used in conjunction with MLMC for the effective application of uncertainty quantification strategies, as they allow to tackle complex geometrical configurations, also with very coarse meshes.

A three-dimensional thermo-hydro-mechanical coupled model for enhanced geothermal systems (EGS) embedded with discrete fracture networks

Fractures play an important role in the geothermal energy recovery from enhanced geothermal systems (EGS) involving interacted multi-physical fields. A strongly coupled thermo-hydro-mechanical (THM) model is proposed to simulate the process of the long-term geothermal production in three-dimensional geothermal reservoirs containing arbitrary discrete fracture networks (DFNs). By introducing a strong discontinuity concept into the DFN model, the aperture variation of each fracture induced by the fluid pressure, external stresses and thermal expansion in the period of production can be captured. Non-isothermal fluid flow in DFN and the local thermal non-equilibrium (LTNE) between fluid and rock matrix are formulated and coupled with the fracture deformation model. Verification is carried out against an analytical solution, followed by a sensitivity and convergence analysis concerning time step and mesh size. This approach is then applied to Habanero EGS project in Australia to evaluate the geothermal productivity and efficiency for a period of 20 years with different injection and production pressures. The results demonstrate that the proposed model is robust and effective to simulate the THM coupled process in 3D fractured reservoirs.

Conforming, non-conforming and non-matching discretization couplings in discrete fracture network simulations

Simulations of fluid flow in naturally fractured rocks have implications for several subsurface applications, including energy storage and extraction, and waste disposal. We are interested in flow in discrete fracture networks, which explicitly represent flow in fracture surfaces, but ignore the impact of the surrounding host rock. Fracture networks, generated from observations or stochastic simulations, will contain intersections of arbitrary length, and intersection lines can further cross, forming a highly complex geometry. As the flow exchange between fractures, thus in the network, takes place in these intersections, an adequate representation of the geometry is critical for simulation accuracy. In practice, the intersection dynamics must be handled by a combination of the simulation grid, which may or may not resolve the intersection lines, and the numerical methods applied on the grid. In this work, we review different classes of numerical approaches proposed in recent years, covering both methods that conform to the grid, and non-matching cases. Specific methods considered herein include finite elements, mixed and virtual elements and control volume methods. We expose our methods to an extensive set of test cases, ranging from artificial geometries designed to test difficult configurations, to a network extruded from a real fracture outcrop. The main outcome is guidances for choice of simulation models and numerical discretization with a trade off on the computational cost and solution accuracy.

Segmented two-phase flow analysis in fractured geological medium based on the numerical manifold method

Abstract A segmented two-phase flow model is developed to describe immiscible two-phase flow in discrete fracture networks; i.e., there is an abrupt interface separating the two immiscible fluids. The capillary pressure on the fluid material interface and gravity are considered. For two-phase flow in fracture networks it is inevitable to form a mixed fluid in some fracture segments. The mixed fluid is described by a homogeneous model, and parameter γ is introduced as a critical value controlling new fluid formation in mixed fluid. The governing equations for two-phase flow are solved with a developed numerical manifold method (NMM) which is verified using the finite element and analytical methods with two examples. Two-phase flow through a single fracture intersection acting as one of the basic elements for complex fracture networks is then analysed. Different inlet and outlet configurations are considered including one-inlet/two-outlets, two-inlets/one-outlet and one-inlet/multi-outlets. For the second configuration, parametric studies of defined parameter γ and boundary conditions are carried out. In the third configuration, capillary pressure and gravity effects are discussed. The developed NMM is applied to the analysis of two-phase flow in a complicated fracture network showing that it is an efficient method for simulating segmented two-phase flow in fracture networks. Complex two-phase flow characteristics in the fracture networks that cannot be addressed in continuum models, e.g., the change in flow direction and transition from static to flow, are also revealed and analysed.

Identifying Non-Darcian Flow and Non-Fickian Pressure Propagation in Field-Scale Discrete Fracture Networks

Non-Darcian flow has been well documented for fractured media, while the potential non-Darcian flow and its driven factors in field-scale discrete fracture networks (DFNs) remain obscure. This study conducts Monte Carlo simulations of water flow through DFNs to identify non-Darcian flow and non-Fickian pressure propagation in field-scale DFNs, by adjusting fracture density, matrix hydraulic conductivity, and the general hydraulic gradient. Numerical simulations and analyses show that interactions of the fracture architecture with the hydraulic gradient affect non-Darcian flow in DFNs, by generating and adjusting complex pathways for water. The fracture density affects significantly the propagation of hydraulic head/pressure in the DFN, likely due to fracture connectivity and flow channeling. The non-Darcian flow pattern may not be directly correlated to the non-Fickian pressure propagation process in the regional-scale DFNs, because they refer to different states of water flow and their controlling factors may not be the same. Findings of this study improve our understanding of the nature of flow in DFNs.