Three-dimensional Discrete Fracture Network Research Articles

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.

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.

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.

Parameterization of a channel network model for groundwater flow in crystalline rock using geological and hydraulic test data

Groundwater flow in sparsely fractured crystalline rocks is highly channelized due to the existing complex hydrogeological heterogeneity in fracture networks. The impacts of multiscale hydraulic heterogeneity on channelized flow in a block of fractured rock under the unidirectional flow condition are investigated by three-dimensional (3D) discrete fracture network (DFN) modeling. A channel network (CN) model generated from DFN is used to model the channelized groundwater flow in fractured crystalline rock. An approach for parameterizing the channel conductance is proposed which uses information from the hydrogeological characterization data. The results show that the network scale hydrogeological heterogeneity dominates the distribution variability of flowrates. The proposed parameterization approach for the channel conductance is effective and robust. Based on the available hydrogeological characterization data, it is possible to compensate for the neglected heterogeneity in the CN model by enhancing the variability of assigned channel conductance. The findings from this work are useful for model simplification from 3D DFN to CN, and for overcoming the difficulty in the parameterization of CN models in applications.

Preferential flow in three-dimensional stochastic fracture networks: The effect of topological structure

This study explores how the topology-based measure of connectivity of a three-dimensional discrete fracture network (DFN) can be used to predict flow channelization, a commonly observed phenomenon from the individual fracture scale to the system scale. First, four-intersection types in fracture intersection classification: I, L, T, and X-intersections, in the DFN topological graph were examined to determine the ones most strongly correlated with preferential flow. T- and X-intersections boosted the formation of preferential flow, and thus, a method of combining their densities to obtain topological connectivity measures was designed. This parameter was tested on the flow behavior of numerous DFNs whose fracture lengths followed a power-law distribution. A method to measure flow channeling was adopted to identify the direct links between the flow channeling and connectivity parameters. The results showed that the network-scale degree of flow channeling decreased as network density increased. However, a strong preferential flow was found at a low DFN density, which was controlled by the network density. More homogenous flow fields were observed in high-density networks, where the main flow paths were determined by the T-intersection and X-intersection densities. Finally, an explicit formula (R-M model) for predicting the flow channeling related to the proposed connectivity parameter was developed. Comparing with published data, the proposed flow channeling estimation model was more reliable than the density-based prediction model.

Estimation of the representative elementary volume of three-dimensional fracture networks based on permeability and trace map analysis: A case study

A three-dimensional (3D) discrete fracture network (DFN) was generated to evaluate the representative elementary volume (REV) of granitic rock masses based on the statistical geometry parameters mapped from outcrop surveys at the Xinchang site, China. Then, extensive numerical simulations were performed to investigate the effect of fracture density on the side length of models of the REV size (LREV). The relationship between LREV of 3D models and the geometric parameters of two-dimensional (2D) cutting planes parallel and perpendicular to the direction of hydraulic gradient variation were analyzed. This relationship was determined by calculating fluid flow through 90 original 3D models, 810 3D submodels, and 270 2D cutting models. Finally, two multivariable regression equations were proposed for predicting the LREV of the 3D fracture network and equivalent permeability of the REV size. The results show that the clustering results obtained from the 3D DFN model are similar to those obtained from field measurements. The values of LREV of the 3D DFN model based on statistical geometry parameters mapped from outcrops are 18.89 m, 19.56 m, and 24.86 m in three flow directions, respectively, which clearly highlights the anisotropic characteristics of fractured rock masses. The LREV decreases with increasing fracture density following an exponential relationship. The multivariable regression function estimates the evolution of the LREV of 3D DFN models with a wide range of fracture densities from 0.53 to 1.33 m/m2. The predicted results are more accurate than those predicted by the fitting function between the fracture density of 2D cutting planes perpendicular to the z-axis and LREV. The proposed model provides a method to approximate the LREV of 3D DFN models and equivalent permeability of the REV size (KREV) using geometric parameters of the 2D cutting planes obtained from trace map analysis of fractured rock masses.

Laminar and turbulent groundwater flows in confined two- and three-dimensional discrete fracture networks

Using laminar or turbulent flow equations indiscriminately to describe groundwater flows in fracture networks may result in large errors. We propose a new method for simulating steady-state groundwater flows in two- and three-dimensional fracture networks by separating laminar flows and turbulent flows in individual fractures. Poiseuille's law is employed when the flow is laminar, while Swamee and Jain formula is utilized when the flow is turbulent. The model results are first compared to recently measured field groundwater flow data in a mining tunnel network. The effect of hydraulic gradient, aperture, fracture width in the third dimension, and fracture density is then examined. Using the laminar flow equation indiscriminately in all fractures overestimates groundwater flowrates, whereas using the turbulent flow equation underestimates them. The aperture contrast of large and small fractures has a significant impact on the potential errors of indiscriminately applying the laminar or turbulent flow equation and the flowrate difference between the three-dimensional and two-dimensional fracture networks. Total groundwater flowrates increase with fracture density in both three-dimensional and two-dimensional fracture networks, with the largest increase occurring when additional fractures are introduced parallel to the hydraulic gradient. The new approach may be utilized to estimate hydraulic head distribution and groundwater flowrate according to the crack patterns and geometric properties of the fractured rock.

Numerical Investigations on the Effect of Fracture Length Distribution on the Representative Elementary Volume of 3D Discrete Fracture Networks

Determination of the representative elementary volume (REV) of fractured rock masses based on equivalent permeability ( K ) is significantly dependent on the geometric characteristics of fractures. In this work, a series of numerical simulations were performed to analyze the relationship between geometric characteristics of fractures and the REV size, in which fracture length follows a power-law distribution. A method to evaluate the K of a three-dimensional (3D) discrete fracture network (DFN) by extracting the equivalent pipe network (EPN) model from the DFN model was utilized and verified. The results show that K of the 3D DFN model has an exponential relationship with the power exponent ( a ) of fracture length distribution and the evaluation of K agrees well with that reported in previous studies, confirming the reliability of the EPN model for calculating seepage properties of complex 3D DFN models. When the side length of submodels ( L n ) is small, the K varies significantly due to the influence of random number seeds used to generate fracture length, location, and orientation. The K holds a constant value after L n exceeds some specific value. The critical model scale is determined as the REV size, and the corresponding volume of the 3D DFN model is represented by V REV . The V REV varies within a narrow range when a ≤ 4.0 . When a = 4.5 , the V REV rapidly increases to more than 3.4 times than that when a = 4.0 . The fluid flow becomes more inhomogeneous due to the small nonpersistent fractures that dominate the preferential flow paths when a exceeds a certain value (i.e., 4.5). The K at the REV size decreases exponentially with the increment of a . This tendency can be explained by the decrease of the average intersection length ( L i ) with the increment of a , which is a geometric parameter for reflecting the connectivity of the fracture network.

Combined Effect of Contact Area, Aperture Variation, and Fracture Connectivity on Fluid Flow through Three-Dimensional Rock Fracture Networks

Abstract In order to investigate the combined effect of contact area, aperture variation, and fracture connectivity on the fluid flow through a fractured medium, a series of flow simulations were implemented on two types of three-dimensional discrete fracture network (3D DFN) models constituting fractures having spatially variable apertures and parallel plates, respectively. The flow tortuosity within the 3D DFN models was examined by changing the density, aperture distribution, and closure of fractures. The results show that compared with the 3D DFN models constituting parallel plates, the model with variable apertures provides more pronounced 3D preferential flow pathways. At the individual fracture scale, the preferential flow pathways mostly converge within the void spaces of large aperture, and at the network scale, they are located in the most transmissive fractures within the connected networks. The permeability of 3D DFNs depends not only on the contact area and aperture variation within individual fractures but also on the fracture connectivity and the contact at fracture intersections within the fracture network. Increasing the fracture connectivity tends to enhance the permeability, while increasing the contact at fracture intersections would significantly reduce the permeability. A correlation between the equivalent permeability of 3D DFNs constituting fractures with spatially variable apertures and parallel plates is proposed incorporating the effect of network-scale topology. A tortuosity factor for 3D DFNs is defined based on the proposed model, and it can account for two competing effects when the model is upscaled from individual fracture to fracture network: the permeability reduction induced by contact obstacles at fracture intersections and permeability enhancement induced by increasing the fracture connectivity.

An Equivalent Pipe Network Modeling Approach for Characterizing Fluid Flow through Three-Dimensional Fracture Networks: Verification and Applications

The equivalent pipe network (EPN) model is an effective way to model fluid flow in large-scale fractured rock masses with a complex fracture network due to its straightforwardness and computational efficiency. This study presents the EPN model for characterizing fluid flow through three-dimensional fracture networks using the Monte-Carlo method. The EPN model is extracted from an original three-dimensional discrete fracture network (DFN) model and is used to simulate the fluid flow processes. The validity of the proposed EPN modeling approach is verified via the comparisons of permeability (k) with analytical solutions and simulation results reported in the literature. The results show that the numerically calculated k using EPN models agrees well with the analytical values of simplified DFN models and the simulation results of complex DFN models. The k increases following an exponential function with the increment of mean length of exponentially distributed fractures (u), which is strongly correlated with fracture density (P32) and average intersection length (Li). The P32 increases in an exponential way with the increment of u. The Li increases as u increases, following a power-law function. The increment of u leads to the increment of a number of long fractures in three-dimensional DFN models. A larger u results in a denser fracture network and a stronger conductivity when the number and length distribution range of fractures remain the same. The representative elementary volumes (REVs) of three-dimensional DFN models with u = 9 m and P32 = 0.4 m2/m3 are determined as 2.36 × 104 m3, 9.16 × 103 m3, and 1.26 × 104 m3 in 3 flow directions, respectively.

A fast numerical method and optimization of 3D discrete fracture network considering fracture aperture heterogeneity

A novel numerical method was developed to evaluate the effect of the aperture heterogeneity of individual fractures on fluid flow through three-dimensional discrete fracture networks (DFNs). First, the rock fractures are modelled as circular discs with sizes, orientations, and positions following exponential, Fisher, and uniform distributions, respectively. Isolated fractures (clusters) that make no contribution to fluid flow through the DFN were then removed by optimization programming using breadth optimization and deep optimization algorithms. Subsequently, individual fractures with heterogeneous apertures were developed using successive random addition algorithms. After that, fracture networks were triangulated, and the Reynolds equations were solved for fluid dynamics computations. The verification showed that the developed method provides advantages in computational accuracy and efficiency. Finally, the influence of the aperture heterogeneity of individual fractures on fluid flow through a DFN was evaluated with the proposed method. For comparison, DFN models with identical mechanical apertures (DFN-I) and another model with heterogeneous apertures (DFNH) were developed. The results show that a preferential flow phenomenon appeared in both the DFN-I and DFNH models. However, owing to the aperture heterogeneity of the individual fractures, the results of the DFNH model were more substantial than those of the DFN-I model. This preferential flow phenomenon would be weakened by an increase in either the mechanical aperture or the fracture density. The results also showed that the permeability difference between the DFN-I and DFNH models was obvious for DFN models with small apertures and low fracture densities, but this permeability discrepancy decreased with increasing mechanical aperture or fracture density. Therefore, a critical mechanical aperture exists for a DFN model with a specified fracture density. Beyond this threshold, the influence of aperture heterogeneity within individual fractures on fluid flow can be neglected, and the DFN-I model can be selected to replace the DFNH model for efficient modelling and computation.

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.

Stability analysis of fractured rock mass around underground excavations based on a three-dimensional discrete fracture network

Stability analysis of fractured rock mass around underground excavations based on a three-dimensional discrete fracture network

Variable resolution Poisson-disk sampling for meshing discrete fracture networks

We present the near-Maximal Algorithm for Poisson-disk Sampling (nMAPS) to generate point distributions for variable resolution Delaunay triangular and tetrahedral meshes in two and three-dimensions, respectively. nMAPS consists of two principal stages. In the first stage, an initial point distribution is produced using a cell-based rejection algorithm. In the second stage, holes in the sample are detected using an efficient background grid and filled in to obtain a near-maximal covering. Extensive testing shows that nMAPS generates a variable resolution mesh in linear run time with the number of accepted points. We demonstrate nMAPS capabilities by meshing three-dimensional discrete fracture networks (DFN) and the surrounding volume. The discretized boundaries of the fractures, which are represented as planar polygons, are used as the seed of 2D-nMAPS to produce a conforming Delaunay triangulation. The combined mesh of the DFN is used as the seed for 3D-nMAPS, which produces conforming Delaunay tetrahedra surrounding the network. Under a set of conditions that naturally arise in maximal Poisson-disk samples and are satisfied by nMAPS, the two-dimensional Delaunay triangulations are guaranteed to only have well-behaved triangular faces. While nMAPS does not provide triangulation quality bounds in more than two dimensions, we found that low-quality tetrahedra in 3D are infrequent, can be readily detected and removed, and a high quality balanced mesh is produced.

Modeling of solute transport in a fracture-matrix system with a three-dimensional discrete fracture network

Understanding the fluid flow and solute transport mechanisms in fractured rocks is essential for many geoengineering applications. In this study, the fluid flow and solute transport in a fracture-matrix system with a three-dimensional (3-D) discrete fracture network (DFN) are modelled through an efficient numerical simulation workflow. The simulation approach is used to systematically investigate the effects of the rock matrix on the transport behaviors in a fracture-matrix system. The results show that the mass exchange between the DFN and the rock matrix can be accurately evaluated based on the conforming mesh at the interface between the fractures (using triangular elements) and the rock matrix (using tetrahedral elements). The complementary cumulative distribution function curves (CCDFs) for the physical processes that consider sorption and decay exhibit significant long tail characteristics, which suggests that the sorption and decay processes play an important role in retarding the migration of solutes in fractured rocks. It is also found that a larger matrix porosity enhances the mass exchange at the interface between the DFN and the rock matrix, which consequently promotes the matrix diffusion effects. The distribution of the concentration plumes in the matrix demonstrates in fracture-matrix systems with larger fracture densities could result in a better connection between the fracture networks and the larger interface (specific wetting) areas, which therefore, promotes the mass exchange. These findings are critical to understanding the migration behavior of radioactive nuclides in far field areas and for the deep geological disposal of nuclear waste.

Nonlinear fluid flow through three-dimensional rough fracture networks: Insights from 3D-printing, CT-scanning, and high-resolution numerical simulations

Nonlinear flow behavior of fluids through three-dimensional (3D) discrete fracture networks (DFNs) considering effects of fracture number, surface roughness and fracture aperture was experimentally and numerically investigated. Three physical models of DFNs were 3D-printed and then computed tomography (CT)-scanned to obtain the specific geometry of fractures. The validity of numerically simulating the fluid flow through DFNs was verified via comparison with flow tests on the 3D-printed models. A parametric study was then implemented to establish quantitative relations between the coefficients/parameters in Forchheimer's law and geometrical parameters. The results showed that the 3D-printing technique can well reproduce the geometry of single fractures with less precision when preparing complex fracture networks, numerical modeling precision of which can be improved via CT-scanning as evidenced by the well fitted results between fluid flow tests and numerical simulations using CT-scanned digital models. Streamlines in DFNs become increasingly tortuous as the fracture number and roughness increase, resulting in stronger inertial effects and greater curvatures of hydraulic pressure-low rate relations, which can be well characterized by the Forchheimer's law. The critical hydraulic gradient for the onset of nonlinear flow decreases with the increasing aperture, fracture number and roughness, following a power function. The increases in fracture aperture and number provide more paths for fluid flow, increasing both the viscous and inertial permeabilities. The value of the inertial permeability is approximately four orders of magnitude greater than the viscous permeability, following a power function with an exponent α of 3, and a proportional coefficient β mathematically correlated with the geometrical parameters.

Transport Upscaling under Flow Heterogeneity and Matrix-Diffusion in Three-Dimensional Discrete Fracture Networks

We investigate the combined effects of network scale flow variability and retention due to matrix-diffusion on the scaling behavior of transport through fractured media. Two of the principal mechanisms controlling the transport of solutes through fractured low-permeability media are broad distributions of flow velocities and retention times in the solid matrix. We study the relative impact of these two processes under different initial conditions using a set of three-dimensional discrete fracture network simulations. We use these simulations to develop and calibrate an upscaled continuous time random walk (CTRW) approach for advective transport based on an Ornstein-Uhlenbeck model for the particle velocities that accounts for the fracture-matrix coupling using a compound Poisson process. This CTRW model can be conditioned on the initial solute distribution and allows to observe late-time scaling behavior at distances beyond what is feasible to observe using high-fidelity direct numerical simulations. We determine that the initial distribution of particles leads to marked differences in the persistent long-term scale behavior in the solute travel time distributions, even those undergoing retention due to matrix diffusion through implementation and analysis of the model.