Flow and Transport in Three-Dimensional Discrete Fracture Matrix Models using Mimetic Finite Difference on a Conforming Multi-Dimensional Mesh

by   Jeffrey D. Hyman, et al.
Los Alamos National Laboratory

We present a comprehensive workflow to simulate single-phase flow and transport in fractured porous media using the discrete fracture matrix approach. The workflow has three primary parts: (1) a method for conforming mesh generation of and around a three-dimensional fracture network, (2) the discretization of the governing equations using a second-order mimetic finite difference method, and (3) implementation of numerical methods for high-performance computing environments. A method to create a conforming Delaunay tetrahedralization of the volume surrounding the fracture network, where the triangular cells of the fracture mesh are faces in the volume mesh, that addresses pathological cases which commonly arise and degrade mesh quality is also provided. Our open-source subsurface simulator uses a hierarchy of process kernels (one kernel per physical process) that allows for both strong and weak coupling of the fracture and matrix domains. We provide verification tests based on analytic solutions for flow and transport, as well as numerical convergence. We also provide multiple expositions of the method in complex fracture networks. In the first example, we demonstrate that the method is robust by considering two scenarios where the fracture network acts as a barrier to flow, as the primary pathway, or offers the same resistance as the surrounding matrix. In the second test, flow and transport through a three-dimensional stochastically generated network containing 257 fractures is presented.



There are no comments yet.


page 1

page 2

page 5

page 10

page 12

page 13

page 15

page 23


An arbitrary order Mixed Virtual Element formulation for coupled multi-dimensional flow problems

Discrete Fracture and Matrix (DFM) models describe fractured porous medi...

A Finite-Volume Moving-Mesh Method for Two-phase Flow in Fracturing Porous Media

Flow in fractured porous media is modeled frequently by discrete fractur...

DuMu^x 3 – an open-source simulator for solving flow and transport problems in porous media with a focus on model coupling

We present version 3 of the open-source simulator for flow and transport...

An optimization approach for flow simulations in poro-fractured media with complex geometries

A new discretization approach is presented for the simulation of flow in...

Variable resolution Poisson-disk sampling for meshing discrete fracture networks

We present the near-Maximal Algorithm for Poisson-disk Sampling (nMAPS) ...

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

We present a topology-based method for mesh-partitioning in three-dimens...

Finite-Volume Simulation of Capillary-Dominated Flow in Matrix-Fracture Systems using Interface Conditions

In numerical simulations of multiphase flow and transport in fractured p...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

There are several computational models for simulating the flow of fluids and the associated transport of dissolved chemicals through subsurface fractured systems, which is of great importance for understanding the behavior of natural subsurface fracture systems and a variety of industrial and civil engineering endeavors, including the environmental restoration of contaminated fractured media national1996rock ; neuman2005trends ; vanderkwaak1996dissolution , aquifer storage and management kueper1991behavior , hydrocarbon extraction hyman2016understanding ; middleton2015shale , longterm storage of spent civilian nuclear fuel hadgu2017comparative ; joyce2014multiscale , and CO sequestration hyman2019characterizing ; jenkins2015state . The computational models currently used today can be loosely grouped into three categories berre2018flow : (1) continuum models - including stochastic continuum neuman2005trends ; neuman1988use ; tsang1996tracer and dual-porosity / dual-permeability models gerke1993dual ; lichtner2014modeling ; zimmerman1993numerical , (2) discrete fracture networks (DFN) models cacas1990modeling ; long1982porous ; maillot2016connectivity ; makedonska2016evaluating ; nordqvist1992variable , and (3) discrete fracture matrix (DFM) models ahmed2015control ; ahmed2015three ; antonietti2016mimetic ; arraras2019mixed ; boon2018robust ; berrone2018advanced ; Flemisch2016 ; fum ; ODS ; lagrange ; manzoor2018interior ; sandve2012efficient ; Schwenck2015 . The distinguishing feature between the categories is the relative fidelity with which the rock matrix and embedded fractures are represented. Fractures are mechanical breaks distinguished from the surrounding rock matrix by the following two fundamental properties: first, there is a stark contrast between their length and their width. The latter of these features, which we shall refer to as aperture, is typically much smaller than the former, i.e., a meter-long fracture can have an aperture of about m depending on the rock type SKB2010 . Second, the hydraulic resistance offered by a fracture can be quite different from the surrounding rock matrix. In turn, fractures can be hydraulic conductors, i.e., higher permeability than the rock matrix, or hydraulic barriers, i.e., lower permeability than the rock matrix.

On one end of the spectrum are continuum models where the fracture network properties are upscaled into an effective continuum and equivalent properties, e.g., bulk porosity and permeability, to represent the resistance offered by both the fracture network and matrix hadgu2017comparative ; hartley ; jackson ; neuman2005trends ; sweeney2019upscaled ; tsang1996tracer . In this case, all features are resolved in a single spatial dimension, typically a two or three dimensional spatial domain. This methodology is best suited for applications where the influence of the fracture network structure is minimal, and an effective permeability captures the bulk behavior of the system. Examples of this situation include a densely fractured medium that behaves similar to a porous media or when the rock matrix permeability offers roughly the same resistance to flow as the fractures. On the other end of the spectrum are DFN models where individual fractures are explicitly represented and the surrounding matrix is not considered berrone2013pde ; berrone2015parallel ; cacas1990modeling ; davy2013model ; davy2010likely ; de2004influence ; dershowitz1999derivation ; dreuzy2012influence ; erhel2009flow ; hyman2014conforming ; hyman2015dfnworks ; pichot2012generalized ; mustapha2007new . In the DFN methodology, the stark contrast in length scales of fractures leads to them being represented as dimensional objects embedded within an dimensional space, i.e., lines in two dimensions and planes in three dimensions. This methodology is best suited for applications where the influence of the fracture network structure is the primary control of flow and transport behavior, and short time-scale interactions with the surrounding matrix can be neglected. An example of this is a sparsely fractured medium with little to no flow in the surrounding matrix. In between the two methodologies is the discrete fracture matrix (DFM) model where both the fracture network and surrounding matrix are explicitly resolved ahmed2015control ; ahmed2015three ; antonietti2016mimetic ; berrone2018advanced ; Flemisch2016 ; fum ; ODS ; lagrange ; manzoor2018interior ; Schwenck2015 . Similar to DFN, the fractures are represented as dimensional objects, but the matrix is dimensional volume surrounding the fracture network that is also considered. DFM models are well suited for situations where the permeability of the matrix is sufficiently high that it cannot be ignored, and the structure of the fracture network plays a key role in determining the flow field, e.g., sandstones or sealed faults. For example, the case of geothermal systems, where flow is primarily through fractures and the rock matrix is commonly low permeability, is an application where DFMs are well suited because the quantity that flows is heat/energy between the two domains. This particular phenomenon of heat/energy transfer cannot be well represented by continuum or DFN models.

DFM models can be thought of as an extension of the DFN methodology where a DFN is first constructed, and then a mesh of the volume surrounding the network is generated. In the DFN methodology, each fracture within the network is assigned a shape, location, and orientation within the domain by sampling distributions whose parameters are determined by a site characterization. However, the size of the domains of interest, typically kilometers, and the cost of sufficiently sampling relevant quantities in the subsurface, both hydraulic and structural, results in a limited amount of data regarding the topological and geometry attributes of the fracture networks bonnet2001scaling ; national1996rock ; zimmerman1993numerical . In practice, these fracture networks are constructed stochastically to sample over the parameter space. This stochastic generation of a fracture network results in a wide range of feature length scales, which can be arbitrarily small if the generation is unconstrained. For these features to be resolved, the local mesh resolution must be smaller than the smallest length-scale, which is computationally impractical for arbitrarily small features. In turn, the automated mesh generation of an arbitrary stochastically generated DFM is relatively infeasible for an unconstrained three-dimensional network composed of many fractures.

Available DFM models can be loosely partitioned into two categories based on how they address the complications associated with mesh generation. DFM models either have a mesh representation of the matrix that aligns/coincides with the mesh of the fracture network, referred to as a conforming mesh method ahmed2017cvd ; bogdanov2003two ; helmig1997multiphase ; karimi2003efficient , or does not, referred to as a non-conforming mesh method berrone2018advanced ; berrone2017flow ; Flemisch2016 ; frih2012modeling ; fum ; Schwenck2015 . This computational geometry choice directs the selection or development of a numerical scheme to resolve flow in coupled fracture/matrix systems flemisch2018benchmarks

, which is complicated by the mixed-dimensional nature of the DFM, where fractures are co-dimension 1 of the matrix. Non-conforming methods bypass the issues associated with mesh generation. In these methods, neither the mesh of the fracture network needs to conform to intersections between fractures, nor the mesh of the matrix needs to conform to the mesh of the fracture network. Although mesh generation is simpler than when using a conforming method, the computational cost is pushed upward via additional degrees of freedom in the resulting discretized systems. To discretize the governing equations for flow and transport in the DFM, non-conforming methods have applied extended finite elements 

del2017well ; Flemisch2016 ; fum ; Schwenck2015 , embedded finite elements ODS , finite elements/boundary elements berrone2018advanced , mixed finite element boon2018robust , mortar-type methods frih2012modeling ; nordbotten2019unified , immersed boundaries berrone2017flow , and Lagrange multipliers lagrange . However, the inclusion of additional physical processes, such as multiphase flow or chemical transport and reactions, commonly requires additional novel mixed-dimensional discretizations for their respective governing equations rather than more straightforward extensions of their unidimensional counterparts. In contrast, the numerical methods for resolving flow and transport in the coupled system when using a conforming mesh are typically simpler. They can have fewer degrees of freedom compared to non-conforming mesh methods. The most commonly used conservative discretization methods for conforming unstructured meshes are control-volume distributed multi-point flux approximations ahmed2015control ; ahmed2015three ; ahmed2017cvd ; manzoor2018interior ; sandve2012efficient , mixed finite elements, and recently mimetic finite difference antonietti2016mimetic . These methods are widely used in non-DFM applications and have seen extensive development, including the construction of higher-order methods for flow and transport and multiphase flow, which is promising for their eventual application to DFM models abushaikha2020fully ; lipnikov2014mimetic . However, the trade-off is that the challenges associated with conforming mesh generation must be addressed. While it is not required that the mesh of the fracture network conforms to intersections between fractures, it considerably reduces the complexity of the discretized governing equations and volume mesh generation when compared to non-conforming methods. To date, most conforming DFM models in the literature are either two-dimensional ahmed2015control ; ahmed2017cvd , quasi two-dimensional ahmed2015three , or relatively simple three-dimensional geometries berre2018flow due to the complexity of mesh generation.

We present a conforming mesh methodology for resolving three-dimensional flow and transport that addresses the key challenges of the DFM approach in a unified manner via (1) a method for mesh generation of and around a complicated 3D fracture network, (2) the discretization of the governing equations for flow and transport using second-order accurate mimetic finite differencing, and (3) implementation of numerical methods for high-performance computing environments. In section 2, we present the governing equations along with our adopted discretization method and the methodology for the mixed dimensional coupling. We then discuss our hierarchical implementation of physical process kernels that allows for both strong and weak coupling of the fracture/matrix domains. In section 3, we address the issues of conforming mesh generation by extending the feature rejection algorithm for meshing (FRAM) presented in Hyman et al. hyman2014conforming to produce a conforming Delaunay tetrahedralization of the volume surrounding the fracture network where the triangular cells of the fracture mesh are faces in the volume mesh. In section 4, we provide two expositions of the method in complex fracture networks. In the first example, we demonstrate that the method is robust by considering three scenarios where the fracture is a barrier, a primary pathway, and offers the same resistance as the surrounding matrix. In the second example, flow and transport through a 3D stochastically generated network containing 257 fractures is presented. General remarks and discussion about the method are provided in the final section 5.

2 Numerical Methods: Governing Equations, Discretization Methods, and Implementation

In this section, we present the governing equations for flow and transport in two and three dimensions and their mixed-dimensional coupling. We then present the conservative discretization schemes and conclude with their computational implementation.

2.1 Governing Partial Differential Equations

The domain of interest is a fractured network surrounding by a porous media. We denote the domain as

. Our focus is on the laminar flow of an isothermal, incompressible, single-phase fluid passing through a porous media that can be described using Darcy’s Law. For simplicity, we present the equations for the case of steady flow. Under this model description, the partial differential equations (PDE) governing flow in



where is the Darcy flux [m/s], is the aqueous pressure [Pa],

is the absolute permeability tensor [m

], is the fluid viscosity [Pa s], is fluid density [kg/m],

is the gravity vector, and

denotes external sources. The conceptual PDE model of the advective transport of non-reactive solute in is


where is the porosity [-]. Equations (1) and (2) along with appropriate initial and boundary conditions describe flow and transport in .

We are interested in solving (1) and (2) when there are fractures present within . Due to the small size of fracture apertures relative to their length, we shall model fractures as lower dimensional objects. We will represent each fracture as a two-dimensional manifold in . We adopt a convention of subscripts and corresponding to the matrix and fracture networks, respectively. We denote the union of all fractures in the domain as and the surrounding rock matrix as . We assume that and are well-defined such that they make up the entire domain () and that they are disjoint (). We denote the surface of that is in contact with as . Note that both and can be disconnected sets.

Within , the mass balance equation and the constitutive law are similar to (1):


Within , the mass balance equation and the constitutive law are given by:


where is the unit normal vector pointing into the fracture domain, is the absolute permeability of the fracture [m], and is the flux between the matrix and fracture domains. The square brackets on the right-hand side denote the difference between values on two opposite sides of a fracture. Thus, the right-hand side contribution represents the mass of water going into/out of the fracture. We impose conventional boundary conditions for both matrix and fracture domains to close the problem. On the matrix-fracture interface the mass conservation law leads to the velocity continuity condition:


In the fracture, averaging the balance equation in normal direction leads to


where is the fracture aperture, and is the mean value of pressure in matrix on two sides of the fracture. A detailed derivation could be found in jaffre-coupled .

Similarly, the conceptual PDE model of the coupled transport of non-reactive solute is


subject to the Dirichlet boundary conditions on the inflow part of the computational domain. is concentration in fracture for outflow () and in matrix for inflow.

For the purpose of this work, we drop external sources and assume that the flow is driven by boundary conditions. Additional physical models supported by our open-source simulator include transport dispersion, energy transport, rock and fluid compressibility, and chemical reactions, but their implementation is not described here.

2.2 Discretization method

We discretize the system of equations (3)-(7) using the mimetic finite difference (MFD) method lipnikov2014mimetic . The MFD method combines flexibility of the finite volume and finite difference frameworks with the theoretical basis of the finite element framework. It allows us to use general unstructured polytopal, polygonal (2D) or polyhedral (3D), meshes to model complex geometric structures that appear in fractured porous media. We apply this method in both and domains.

To guarantee the local mass conservation property, the MFD method uses the mixed formulation and explicitly discretizes pressure and flux field unknowns. The MFD method operates with two types of pressure unknowns: (1) cell-based, denoted with a subscript in matrix and in fractures residing at cell centers (, ), and (2) face-based, denoted with a subscript residing at cell faces ( and ). In , there is one cell-based pressure unknown per volume cell and one face-based value per cell face. In , there is one cell-based pressure unknown per fracture cell, and one face-based value per fracture cell face. Cells of a fracture mesh in coincide with mesh faces of partition. To be able to approximate a pressure jump between two sides of a fracture surface MFD uses two face-based pressure unknowns, that approximates pressure from both sides of a fracture. The flux field is expressed in terms of normal components of the flux, and , which are collocated with the and unknowns, correspondingly. Similar to the two pressure degrees of freedom placed on both sides of fracture, we introduce two flux unknowns as well. Figure 1 provides a diagram with the positions of the unknowns on a matrix cell (tetrahedron right) and a fracture cell (triangle left).

Figure 1: Diagram with the positions of the unknowns on a matrix cell (tetrahedron right) and a fracture cell (triangle left). Note the the pressure unknown at the center of the fracture face is collocated with the unknown of pressure on one of the matrix cell faces. : Pressure at the center of matrix cell. : Pressure at the center of fracture cell. : Pressure on the boundary of the matrix cell. : Pressure on the boundary of the fracture cell. : Normal components of the flux on boundaries of the matrix cell. : Normal components of the flux on boundaries of the fracture cell.

The discrete equations for the flow equations (3)-(4) can be written locally for each matrix and fracture cell ( and , respectively) and then coupled and assembled into the global system. The local discrete equations have the following saddle-point forms:




where ,, , are vectors of degrees of freedom which correspond to a particular matrix cell, , or fracture cell, . The construction of the square mass matrices and is the key ingredient of the MFD construction, which is explained in detail in lipnikov2014mimetic . These matrices are constructed to guarantee the second order of accuracy of a numerical solution on unstructured polyhedral meshes with a full anisotropic permeability tensor. represents fluxes between and :


The local systems are coupled by the flux continuity conditions inside and and on :


From an implementation viewpoint, it is convenient to treat flux continuity conditions (11) as the flux boundary conditions for the matrix domain. This leads to uniform code implementation in interior and near-boundary mesh cells. Recall, that denote face-based degrees of freedom on the opposite sides of fracture cell . The last formula shows the convenience of using face-based pressure unknowns in the numerical scheme. The coupling with the fracture domain and the approximation of a pressure jump is included into the flux continuity conditions (11).

Using the static condensation algorithm MFD2006 the flux degrees of freedom can be excluded locally which converts local systems from the saddle-point form to symmetric positive definite (SPD) formulations:




where entries of are either zeros or one of the terms in .

Equations (10)-(11) can be assembled into the global coupled system with a symmetric positive definite matrix:


where the -vectors collect all respected local pressure unknowns. The final step is to apply the Dirichlet boundary conditions. This means to prescribe known values to and located on the domain boundary and to eliminate them from the global system.

To discretize the system of hyperbolic transport PDEs (7), we apply the finite volume method. The solute concentration is discretized by cell-centered degrees of freedom for volumetric matrix cells and surface fracture cells. The advective solute fluxes are approximated using the first-order upwind method (see, e.g. Barth:91 ). The solute flux is defined similarly in and . Consider and some mesh face with the corresponding Darcy flux which is one of the defined above. Then, the solute flux is approximated using the upwind scheme:




and are solute concentrations in the cells on both side of a face . Here, we assume that the normal points from cell to cell . We denote the right-hand side of equation (15) as by analogy with the flux definition.

Using the first-order time approximation, we replace equations (7) with discrete equations


Setting , we obtain the explicit scheme. The implicit scheme is obtained by setting . A second-order advection scheme is derived within the simulator using more accurate approximation of the solute flux based on a limited linear reconstruction of the concentration. A second-order accurate time integration is implemented via the Runge-Kutta scheme.

2.3 Implementation: Process Kernels

We implemented the coupled matrix-fracture flow system in the open-source code Amanzi ASCEM_Amanzi_design_1 which provides a suite of critical low-level components and is built on the flexible and extensible Arcos process management system ETCoon_JDMoulton_SLPainter_2016a . Amanzi provides a flexible parallel mesh infrastructure that supports the extraction of meshes on the matrix and fracture sub-domains, i.e., and . It provides an operators library that contains MFD and FV discretizations of the flow and transport equations. In addition, it includes a suite of nonlinear and linear solvers KLipnikov_DMoulton_DSvyatskiy_2016a , and provides interfaces to powerful third-party solvers such as Hypre’s algebraic multigrid.

Arcos, the underlying framework of Amanzi, provides services to manage process-rich systems with many coupled models, including both PDE and parametric models. At a high level, Arcos defines two key building blocks or objects, the process kernel (PK) and the multi-process coordinator (MPC). The process kernel provides an implementation of an API for a single process, e.g., a single PDE, on a specific mesh. The MPC provides an implementation of a particular coupling strategy using this same PK API to assemble the PKs and MPCs hierarchically in a PK-tree to represent the coupled system. A PK-tree for the coupled matrix-fracture flow and transport system is shown in Figure

2. The tight coupling of the matrix-fracture flow system shown there is realized, at each time step, by solving the global matrix formed by combining the associated matrix and fracture matrices. The vector and operator abstractions in Amanzi make this tight coupling relatively easy to implement. In contrast, the matrix-fracture transport can be weakly coupled, with implicit diffusion and explicit advective transport in the matrix and the fracture, or strongly coupled with the implicit solution of a single advection-diffusion system formed by combining the matrices from the matrix and fracture systems. The top-level or base MPC shows that these flow and transport MPCs are coupled sequentially. The flow time step is usually much larger than a stable and accurate transport time step, so the transport MPC is sub-cycled relative to the flow MPC.

At a lower level, Arcos dynamically manages a complete and self-consistent view of all variables in the system, dubbed the state. For example, the instantiation of PK or MPC registers its data requirements, i.e., pressure on cells of the matrix mesh, with the state data manager so that Arcos can construct a directed acyclic graph (DAG) of all data dependencies in the system. Acros uses this DAG to manage updates of all variables in the system, eliminating code from PKs and MPCs that would depend on their relative position in a hierarchy.

Figure 2: The Amanzi Process Kernel (PK) tree is shown for the matrix-fracture flow and transport system. Flow and Transport are coupled sequentially. Flow System: Flow Coupler - tight coupling forms and solves matrix-fracture flow as a single system. Flow PK on matrix mesh generates sub-system. Flow PK on fracture mesh generates sub-system. Transport System: Transport Coupler: tight coupling solves matrix-fracture transport implicitly. Can sub-cycle time-steps relative to flow. Transport PK on matrix mesh generates sub-system. Transport PK on fracture mesh generates sub-system.

Accuracy and efficiency of Amanzi framework was successfully tested in the benchmark studies berre2021verification . Amanzi code was one of a few participants that successfully simulated all proposed scenarios ranging from one fracture problem up to field scale scenarios based on DFM network with 52 stochastically generated fractures.

3 Conforming Mesh Generation

The discretization provided in Section 2 requires that the two-dimensional mesh of is the trace of the three-dimensional mesh of . In this section, we describe our methodology to accomplish this task.

3.1 DFN Meshing

We generate and mesh each fracture network using the Feature Rejection Algorithm for Meshing (FRAMhyman2014conforming implemented within the dfnWorks computational toolkit hyman2015dfnworks . FRAM produces a conforming Delauney triangulation of a DFN without user intervention or manual mesh manipulation by coupling the generation of the networks with meshing. The cornerstone of FRAM is a user-defined minimum length scale () that determines what geometric features will be represented in the network and is conceptually equivalent to the mesh resolution. FRAM constrains the network generation so that the local feature size on each fracture is greater than . This constraint ensures that pathological cases that degrade mesh quality do not exist and thus provides a firm lower bound on the required resolution of the mesh. In turn, the network can be meshed using cells with edges less than . The algorithm presented by Murphy et al., murphy2001point is used to produce a conforming Delaunay triangulation where contiguous edges in the mesh represent the lines of intersections between fractures. We refer the reader to Hyman et al. hyman2014conforming for details concerning fram. While fram and dfnWorks can generate meshes with either variable or uniform resolution, we only consider a uniform resolution mesh as it simplifies the extension to the creation of a conforming volume mesh.

3.2 DFM Meshing

The meshed DFN produced by fram is the input for our volume meshing strategy. The volume mesh generation workflow presented below is an automated approach to filling the volume around and between the triangulated planar polygons of the DFN with a tetrahedral mesh that conforms to the triangle facets of the DFN. While it is not a requirement for the discrete formulation presented in Sec. 2 the triangulations that tessellate each planar fracture polygon are Delaunay triangulations and the tetrahedral mesh that fills the volume is also a Delaunay mesh. The workflow described is a strategy to generate a distribution of vertices in 3D such that a Delaunay tetrahedralization of the vertices will result in a mesh with tetrahedron edges that conform to the lines of intersection between each planar fracture and a mesh with tetrahedron triangle facets that conform to the triangulation of the planar polygon fractures. The theoretical basis for the approach is proven in Murphy et al. murphy2001point . The method is presented for the case where the triangulated DFN has uniform resolution, the edge lengths of the triangles have similar length. There is no theoretical barrier for the method to be applied when the mesh resolution varies from coarse to fine.

The method is broken up into three primary steps outlined below. Figure 3 provides a visual depiction of the workflow and psuedo-code is provided in Algorithm 1 presented in A. We use a network composed of two rectangular fractures for visual exposition of the workflow, but the algorithm is general for an fracture DFN.

  1. Generate a uniform resolution conforming Delaunay triangulation of a DFN. Our example case of a two fracture DFN, one green and one blue, with the mesh overlaid in black is shown in Fig. 3 (a). The fractures are triangulated such that the triangulations are conforming, i.e., share edges, along the lines of intersection. This mesh is generated using the method provided in Sec. 3.1. Triangular edges in the mesh are of length .

  2. Create a background mesh, which is a directionally isotropic tetrahedrons mesh with cell edge length

    , which matches the size of the DFN triangles. We generate the background mesh by estimating the number of vertices necessary to produce a uniform sized mesh with the desired spacing and generate a distribution of vertices in a cubic pattern. This point distribution is connected to form Delaunay tetrahedrons, which is then smoothed using elliptic smoothing that guards against inverting cells 

    khamayseh1996anisotropic . Then, a final set of smoothing steps is performed using a three-dimensional implementation of the minimum error gradient adaption algorithm that creates a smoothed grid that is adapted to the standard function with constant Hessian bank1997mesh . The first smoothing step breaks the symmetry of the initial orthogonal mesh, which is a false local minimum of the smoothing algorithm, and the second smoothing step converges to a high quality, isotropic mesh. This methods also contains features that prevent inverting cells during vertex movement. However, these smoothing steps do not maintain a Delaunay mesh. Thus, the final step is to throw away the connectivity and then reconnect the vertices using a three-dimensional Delaunay tetrahedralization algorithm lagrit2011 . The resulting mesh is shown in Fig. 3-(b). This approach, while being effective at creating a high quality, isotropic, 3D Delaunay mesh, is inefficient with respect to use of computing resources kuprat1996adaptive . Other methods, such as Poisson disk sampling in 3D guo2016tetrahedral ; krotz2021maximal may prove to be a more efficient approach, but has not been implemented or utilized for the calculations presented here.

  3. Excavate the background mesh. We superimpose the triangle DFN and the tetrahedron mesh to compute the circumscribed sphere of each triangle in the DFN mesh, which are shown in Fig. 3 (c). If a vertex of the tetrahedron mesh falls within any circumscribed sphere it is removed along with any tetrahedron associated with the vertex. We call this the excavated mesh, which is shown in Fig. 3 (d). This criteria, that every edge or triangle facet of the input DFN that must have an empty circumsphere, is a sufficient condition so that the output DFM Delaunay mesh conforms to those edges and facets murphy2001point . In order to maximize mesh quality, the distance field between vertices of the DFN and the tetrahedron mesh is computed and any vertex within distance is also removed. This last step is not necessary for obtaining the desired conforming Delaunay mesh, but experience has shown that the mesh quality is improved in terms of tetrahedron aspect ratio.

  4. Create DFM mesh. We create a new point cloud that is the union of the DFN vertices and remaining vertices in matrix mesh. Vertices in this point cloud are then connected using a three-dimensional Delaunay tetrahedralization algorithm. The resulting mesh is shown in Figure 3 (e). These previous steps and the theory presented in Murphy et al. murphy2001point ensure that (i) all the lines of intersection between the fractures will exist as edges of the tetrahedron mesh and (ii) the triangles of the DFN mesh will exist as faces of the cells of the tetrahedron mesh. Figure 3 (f) is a close up of the mesh along the line of intersection between fractures to show that the tetrahedral mesh conforms to the DFN mesh as well as the line of intersection. The final step required for setup of the physics model is to determine this correspondence between the triangles in the DFN mesh and the tetrahedron faces in the DFM, which is a matter of book keeping.

Figure 3: Visual explanation of the DFM meshing workflow. (a) DFN composed of two fractures with the Delaunay Triangulation overlaid in black. (b) Isotropic tetrahedral mesh whose faces cells are the same resolution as the DFN mesh. (c) Red spheres with radii are placed onto every node in the mesh of DFN. (d) Tetrahedral mesh where vertices within the spheres have been removed. (e) Reconnected tetrahedral mesh where the DFN cells are now faces in the tetrahedral mesh. (f) Close up view of the hybrid DFM mesh which is a Delaunay tetrahedralization where the fracture triangulation is represented as faces of the tetrahedrons.

The workflow outlined above is robust by design, especially when the DFN is sparse and the angle between the fractures is sufficiently large. However, we acknowledge the method presented here is not perfect for any arbitrary networks. It is possible, albeit rare (typically ), for few facets of the original DFN to not be maintained along the fracture intersections. The reason for this is that while we ensure that no vertices of the DFM are within the circumsphere of the DFN triangles, we do not ensure that a vertex from a nearby intersecting fracture triangle is not within any DFN circumsphere. For most physical scenarios where a DFM model would be appropriate, this is a small part of the overall coupled fracture-matrix system, and the effect on flow, transport, and other physical quantities of interest is negligible. Nonetheless, we propose two optional extensions of methods to address these interference points in the meshing procedure. The first is a three-dimensional extension of FRAM and the second is a graph-based/upscaling method that reduces the network complexity while preserving key hydrological properties of the fracture and matrix system.

3.3 A Three-Dimensional Extension of the Feature Rejection Algorithm for Meshing

The version of FRAM presented in Hyman et al. hyman2014conforming used to create the DFN mesh does not consider features in three-dimensional space that affect the ability to automatically generate a conforming volume mesh of the surrounding matrix. Examples of three-dimensional features that influence volume meshing, but not that of the DFN, include the distance between two fractures that do not intersect and the angle of intersection between two fractures. Suppose two fractures are arbitrarily close but do not intersect. In that case, the volume cell between the two planes must have edge lengths smaller than that distance, which is impracticable if the minimum distance is unconstrained.

To address these issues, one can extend FRAM to constrain the generation of the network to adhere to a minimum feature size in three-dimensions and automatically produce a conforming volume mesh using the method presented above. While this methodology is straightforward to implement, it is computationally costly as it must be exhaustive. If the network is sufficiently sparse and contains a few primary features, e.g., large faults, then this method will produce a network that can be automatically meshed, i.e., without user manipulation, However, as networks become denser, it becomes impracticable to constrain the network so that all three-dimensional features are larger than while adhering to geological observations and keeping the number of mesh cells to a reasonable number.

3.4 Network Simplification Through Primary Network Identification

Flow channeling, isolated regions of high velocity within the flow field, is commonly observed in field observations, laboratory experiments, and numerical simulations in fracture media abelin1991large ; abelin1985final ; hyman2020flow ; dreuzy2012influence ; rasmuson1986radionuclide . These observations indicate the existence of primary sub-networks, also referred to as the network backbones, where the majority of flow and transport occurs. These backbones are commonly composed of a few large fractures connected by small fractures, and they tend to be structurally less complicated than the entire network due to their reduced size. The formation of flow channels suggests that in many applications, it may be unnecessary to explicitly resolve all fractures in the network to effectively capture the system’s bulk flow and transport behavior. The method proposed below exploits this feature by partitioning a network into primary and secondary subnetworks where the DFM volume mesh conforms to the mesh of the primary network and influence of the secondary network is accounted for using upscaled effective parameters.

The conceptual definition of a backbone is relatively simple but formulating an operational definition to partition a fracture network depends on various variables. Formally, we seek to partition in two disjoint subsets, the backbone and the secondary structure such that they make up the whole network and they are disjoint . Backbone membership has been shown to depend on the macro-scale structure of the network, e.g., orientations and density hyman2019linking ; sherman2020characterizing ; hyman2020flow , meso-scale hydrological attributes, e.g., fracture permeability de2002hydraulic ; hyman2016fracture as well as the flow domain boundary conditions grindrod1993channeling ; neuman2005trends . There are a number of methods to identify networks backbones aldrich2017analysis ; hyman2017predictions ; hyman2018identifying ; maillot2016connectivity ; wellman2009effects In principle, any of these methods can be applied to identify and its compliment. Once the backbone is identified, we shall generate a volume mesh that only conforms to the mesh of the backbone. This process reduces the complexity of generating a volume mesh around the entire network because fewer small local features that lead to degradation of mesh quality due to interference from nearby points exist. Once the backbone is identified and the matrix surrounding it meshed, we use the methods developed by Sweeney et al. sweeney2019upscaled to upscale the properties of the secondary network into equivalent properties of the matrix mesh. We refer the reader to the original paper for complete details, but we highlight the key aspects here. The basis for calculating equivalent permeabilities and porosities of the secondary fractures is to calculate their intersection areas with the matrix tetrahedron. Once the intersection areas are known, the porosity and permeability of each intersected cell can be calculated using the fracture apertures and orientations. This procedure is done for each tetrahedron cell in the matrix mesh and allows for properties of multiple fractures (intersecting fractures) to be upscaled to the same cell.

3.4.1 Network Simplification: Example

Figure 4: (Top Left) DFN composed of 257 Fractures. Mesh has 142841 cells and 74418 vertices. (Top Right) Single fracture from within network shown on left. 1790 cells and 949 vertices. (Bottom Left) Primary subnetwork/backbone. (Bottom Right) Secondary sub-network.

Figure 4 (top left) provides a example meshed DFN generated using fram, which we will use to demonstrate the network simplification method. The domain is a cube with sides of length 10 meters. The fractures are squares whose lengths are sampled from a truncated power-law distribution with an exponent of 2.4, lower cutoff of 1.5 meters, and upper cutoff of 4.5 meters. The fractures in the network are sampled from three different families whose mean orientations align with the principal Cartesian axis. There is slight variation around the mean direction in each family by sampling the normal vectors using a von-Misser Fisher distribution with dispersion parameter  wood1994simulation . The resulting network is made up of 257 fractures. In this network, is set to 0.25 m. Overlaid on the DFN is the conforming Delaunay triangulation, which contains 142841 triangular cells and 74418 vertices. Figure 4 (top right) shows the mesh of a single fracture extracted from the network. The mesh is isotropic and composed of uniformly size triangles. This particular fracture contains 1790 cells composed of 949 vertices. Vertices along the lines of intersection are blue spheres. Cells whose edges align with the line of intersections are yellow, and the remaining cells are colored orange. This particular fracture intersects eight other fractures and includes one point where there is an intersection of intersections on the fracture plane, i.e., a triple intersection. Note that the lines of intersections, including triple intersections, are represented in the triangulation.

Figure 4 (bottom left) shows network backbone which contains 135 fractures and the secondary structure is shown in Figure 4 (bottom right). This particular backbone was identified by recursively removing all dead-end fractures, i.e., isolating the 2-core of the graph based on the DFN hyman2017predictions . Notice that the backbone is much simpler than the entire network. The DFN backbone is meshed with 53103 vertices and 103068 cells, a reduction of about 30%. This method is one of the least aggressive backbone identification methods and retains a large portion of the network. Other methods have been used to reduce the network to less than 5% of the original number of fractures with good agreement in particular transport attributes hyman2017predictions ; hyman2018identifying . This network reduction/simplification facilitates the automated mesh generation of the surrounding matrix, which is shown in Fig. 5. The resulting volume mesh contains 2,902,906 tetrahedrons and 501,298 vertices. The blue cells in the mesh correspond to the variable permeability of the upscaled matrix cells, which effectively capture the secondary structure of the DFN. While we used the method presented by Sweeney et al. sweeney2019upscaled , any upscaling method that works with unstructured tetrahedrons can be employed. Note that the secondary structure does not need to be meshed, and there are no constraints on its size or complexity. This aspect thereby further reduces the computational burden of the overall method to create the hybrid mesh.

Figure 5: Resulting DFM mesh produced by extracting the primary subnetwork (red fractures) and upscaling the secondary structure whose affect on flow is captured via upscaling. The DFM volume mesh contains 2,902,906 tetrahedrons and 499,333 vertices. The primary subnetwork of the DFN is meshed with 53103 vertices and 103068 cells.

The principal parts of the methods are described in Algorithm 1.

  1. Generate a two-dimensional discrete fracture network

  2. is partitioned into primary and secondary subnetworks such that and they are disjoint .

  3. The primary fracture network is meshed with a uniform resolution Delauney triangulation that conforms to the fracture intersections denoted .

  4. The rock matrix is meshed using Delauney tetrahedrons, where the faces of the tetrahedrons conform to the mesh of the fractures denoted .

  5. The secondary network is upscaled into effective properties in the matrix volumes .

  6. - A multi-dimensional mesh representation of the domain is discretized using mimetic finite differences for the governing equations for flow and transport on the multidimensional system and .

Algorithm 1 DFM generation and meshing

4 Methodology Verification and Exposition

In this section, we provide three examples of the method. The numerical methods presented in this work have previously been verified for single-phase flow and transport problems in DFM systems in berre2021verification . In that study, four benchmark cases for single-phase flow and transport in three-dimensional fractured porous media were given to 11 participating groups, including the Amanzi team. In each of the problems, Amanzi successfully produced accurate results for the solutions. However, that analysis was completed using provided mesh files and not the meshing algorithm developed in this work. As a result, we first present verification for transport in a single fracture coupled to a matrix here at various resolutions for further completeness. Next, we consider a smaller network composed of five fractures. Thereon, we resolve flow and transport under a set of different fracture and matrix properties ratios to demonstrate the robust nature of the method. In the final example, we use the larger network provided in section 3.4.1 that originally contains 257 fractures before the backbone is calculated, which contains 135 fractures. Here, we demonstrate the scalability of the method to more complicated networks and the effectiveness of the backbone/upscaling method.

4.1 Example 1: Single Fracture

In this particular example, shown in Figure 6, the single fracture is parallel to the direction of flow. We prescribed Dirichlet pressure boundary conditions of Pa to the inlet face and Pa to the outlet face. In both cases, the fracture and matrix domain are given the same boundary conditions. The porosity of both the fracture and the matrix are arbitrarily chosen as 1.0 and 0.01, respectively. The permeability of the matrix is isotropic and given a value of m. The fracture aperture is m, hence the isotropic permeability is given by the cubic law, and is m. The domain is 1 m.

The transport of passive solute in the domain is modeled using the ADE equation (17) once a steady state pressure solution is achieved, which is shown in Figure 6a. We simulate a pulse of solute injected into both the fracture and matrix domains at s using a Heaviside function for s. The expected arrival time of the solute at the opposite boundary to the inlet face can be predicted using Darcy’s Law. In this case, we would expect two peaks of solute in the breakthrough PDFs, one at s due to the solute breakthrough from the fracture, and the other at s due to the solute breakthrough from the matrix.

Figure 6: Verification simulation of a single fracture coupled to matrix. (a) Steady state flow solution. (b) Early time transport of solute in fracture. (c) Late time transport of solute in matrix. (d) Breakthrough curves of solute in both fracture domain (red) and matrix domain (blue) . Analytical expected values calculated from Darcy’s Law are shown by vertical lines.

Figure 6 shows the flow solution (a), two snapshots of solute transport (b) and (c), and normalized breakthrough curves of the solute (d). In (d), the red line corresponds to transport in the fracture and the blue line corresponds to transport in the matrix. The curves are normalized to the peak breakthrough in the fracture. It is clear that the peaks agree with the analytical expected values (the vertical lines). In each case there is numerical diffusion due to the discretization and order of the transport scheme.

4.2 Example 2: Five fracture Network

Our next example network is made up of five square fractures within a 10-meter cube, shown in Figure 7. Each fracture has a side of length 4 meters shown in blue. The fracture orientations are not orthogonal to one another and their intersections form obtuse and acute angles. Therefore, they are not aligned with the background tetrahedral mesh. The mesh has a uniform resolution with triangle edge length of 0.5 meters, which corresponds to an value of 1 m. The mesh of the DFN has 1099 cells and 628 vertices and is colored black. The matrix mesh is made up of tetrahedrons with an edge length of 0.5 meters that is made up of 17861 cells and 3847 vertices. There are no regions of interference in this network, and the matrix mesh is perfectly conforming to the DFN mesh.

Figure 7: Example network composed of five fractures. There are no regions of interference in this network, and the matrix mesh is perfectly conforming to the DFN mesh even though the fractures are not aligned with the background hexahedron mesh.

4.2.1 Transport Verification

Note that the mesh here is not grid aligned as is the previous verification test. Specifics of the meshes are reported in table 1. To verify the second order-accuracy of the transport scheme, we compare the tracer breakthrough curves at the outlet of the domain. We set the fracture permeability to be greater than the matrix permeability: m, which corresponds to an aperture of m, and matrix permeability of m. We consider steady flow via a pressure difference across the x-direction of the domain (0.1 MPa) and a pulse injection into the inflow face. For each refinement level, we take the maximum difference between that level and the finest level of refinement (h = 0.25 m) over the entire transport breakthrough curve at the outlet face. The results are plotted in Fig. 8 along with lines showing first order and second order convergence. The data matches the second order convergence line thereby confirming the derivation and implementation of the numerical scheme for transport.

h # DFN Triangles # DFM Tetrahedrons Maximum Error from h = 0.25m
1 1088 44189 1.18e-6
0.9 1364 59111 1.04e-6
0.8 1695 88136 8.56e-7
0.7 2183 124345 6.80e-7
0.6 2970 202275 5.35e-7
0.5 4266 369361 4.12e-7
0.4 6674 716722 2.69e-7
0.3 11738 1686799 8.99e-8
0.25 16965 2991156 -
Table 1: Information about various mesh resolutions for the 5 fracture example shown in Fig.7. Convergence plots for transport though the DFM are shown in Fig.8

Figure 8: Transport convergence/second-order verification in the network shown in Fig.7. Additional mesh information is presented in Table 1.

4.2.2 Variable Hydraulic Properties

We perform three flow and transport simulations with various contrasts between fracture and aperture permeability. We perform these simulations on the h = 1 m case. In all cases, there is steady flow created by a pressure difference across the domain aligned with the x-axis. The fluid is set to be isothermal water. In all cases, we fix the fracture porosity at 1.0 and the matrix porosity at 0.1. In the first case, we set the fracture permeability to be equal to the matrix permeability: m, which corresponds to an aperture of m. In the second case, we set the fracture permeability to be greater than the matrix permeability: m, which corresponds to an aperture of m, and matrix permeability of m. In the third case, we set the fracture permeability to be less than the matrix permeability: m, which corresponds to an aperture of m, and matrix permeability of m. Pressure solutions in the three scenarios are shown in Fig. 9 (upper left) case 1: , (upper right) case 2: , and (lower left) case 3: . In case 1, there is a smooth pressure gradient through the entire domain, and no irregularities caused by the presence of the fractures are observed. In case 2, the higher permeability of the fractures results in a pressure conduit that forms through the network. The pressure in the matrix is altered accordingly to push flow into the fractures rather than continue through the lower permeable matrix. In case 3, the fractures act as a barrier. Specifically, observe the lower right corner of the domain where there is a pressure back up against the fractures and abrupt changes in the pressure through the fracture plane.

We also consider passive solute tracer transport of a plus injection along the inflow boundary. The normalized cumulative tracer breakthrough curves at the outlet plane for all three cases are shown in the lower right subfigure of Fig. 9. The line style indicates case number (case 1- solid, case 2-dashed, case 3-dotted), and the color corresponds to the region of the domain (whole face, matrix, or fractures) the tracer exits through. The curves are normalized so that the whole face values are equal to 100% of the tracer at the end of the simulation. For case 1, the tracer exits through both the fracture and matrix. The fracture tracer is slightly delayed because of the higher porosity, which results in lower transport velocity. For case 2, all of the tracer exits through the fracture, which is to be expected with the pressure solution and simulation design. For case 3, all of the tracer exits through the matrix, which is to be expected given the lower permeability of the fracture.

Figure 9: Network made up of five square fractures embedded within a permeable matrix. (Upper Left) Fractures and matrix has the same permeability. (upper right) Fractures are more permeable than the matrix. (lower left) Fractures are less permeable. (lower right) Solute Transport breakthrough curves in the three scenarios.

4.3 Example 3

Our final example is the case of the mesh previously shown in Section 3.4.1 that originally contains 257 fractures before the backbone is calculated, which contains 135 fractures. We consider two cases: one in which the 122 fractures not included in the backbone are upscaled to effective properties in the matrix mesh, and another in which they are not. In each case, the multi-dimensional DFM mesh is geometrically identical and we simulate the same single phase flow and transport simulation with the same boundary and initial conditions apart from the upscaled properties. Specifically, we consider isothermal flow and transport of water in the -direction driven by a pressure gradient where the left boundary is given a Dirichlet value of 1.1 MPa and the right boundary is given a Dirichlet value of 1.0 MPa. As in Section 4.2, transport is simulated by a pulse of non-reactive tracer.

Figure 10: (a) Steady state flow solution (b) transport of solute at s (c) breakthrough curve of solute exiting the fracture networks (d) breakthrough curve of solute exiting matrix (e) breakthrough curve of total solute exiting domain.

Fig. 10 a,b shows snapshots of the the steady-state flow solution and transport at s for the case where upscaled fractures are included. In Fig. 10 c,d,e the breakthrough curves of the solute through the right boundary face are plotted, comparing the two cases. It is clear that in each case, the breakthrough from the fracture network is equivalent, which is expected because the fracture networks are identical between the two simulations. However, there is a slight deviation in the matrix. When upscaled fractures are included in the matrix properties, there is a slightly earlier breakthrough of solute, as well as slightly more (10%) solute that passes through the outlet. This phenomenon can be attributed to disparate zones of higher permeability from the upscaled fractures. However, these differences in the matrix ultimately do not affect the total breakthrough since in this particular network the fracture network still dominates the flow and transport. How important it is to include the upscaled fractures will depend on the network and the hydraulic properties of the system and how aggressive of a backbone algorithm is chosen to include in the DFM.

5 Remarks

We have proposed a complete workflow to simulate single-phase flow and transport in fractured porous media using a discrete fracture matrix approach. The method is based on a method for mesh generation of and around a complicated 3D fracture network. We addressed the issues of conforming mesh generation to produce a conforming Delaunay tetrahedralization of the volume surrounding the fracture network where the triangular cells of the fracture mesh are faces in the volume mesh. We provide two extensions of the base method to address pathological cases that commonly arise and degrade mesh quality. The first is an extension of the feature rejection algorithm for meshing (FRAM) presented in Hyman et al. hyman2014conforming which ensures that all three-dimensional length scales are larger than a user-defined length-scale. The second method partitions the network into a primary structure, through which the majority of flow passes through and to which the mesh conforms, and the secondary structure, whose effects on flow and transport are upscaled into the spatially variable permeability field of the surrounding rock matrix.

This high quality conforming mesh facilitates the efficient discretization of the governing equations for flow and transport, which is carried out using second-order accurate mimetic finite differencing. The implementation of these numerical methods for high-performance computing environments is performed using a hierarchically designed multi-process coordinator and process kernel structure that allows for both strong and weak coupling of the fracture/matrix domains. Additional physical models beyond simulate single-phase flow and transport supported by the adopted open-source simulator include energy transport, rock and fluid compressibility, and chemical reactions.

We provided verification tests of the method compared with analytic solutions for flow and transport. We also provided two expositions of the method in complex fracture networks. In the first, we demonstrated that the method is robust by considering scenarios where a network of five fractures is a barrier to flow, a primary pathway, or offers the same resistance as the surrounding matrix. In the second example, flow and transport through a fully three-dimensional stochastically generated network containing 257 fractures are presented. The methods presented here will be made freely available in the open-source code dfnWorks in a future release.


J.D.H. and M.R.S. gratefully acknowledges support from the LANL LDRD program office Grant Number #20180621ECR. J.D.H, M.R.S., and J.D.M gratefully acknowledges support from the LANL LDRD program office Grant Number #20220019DR. M.R.S. would also like to thank support from the Center for Nonlinear Studies. J.D.H and M.R.S. thank the Department of Energy (DOE) Basic Energy Sciences program (LANLE3W1) for support. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001). LAUR # LA-UR-21-31458.

Appendix A Meshing Algorithm

DFN with intersection conforming Delaunay triangulation with cell edge length
DFM Delaunay tetrahedralization with cell edge length and triangle faces conforming to
Step 1: Create Background Mesh
Generate a directionally isotropic tetrahedrons mesh with cell edge length
Step 2: Excavate Background Mesh
for all  do
end for For each triangle in the DFN mesh, compute the circumscribed sphere  
if   then
end if If a vertex of the tetrahedron mesh falls on or within any circumscribed sphere, then it is removed along with any tetrahedron associated with the vertex to produce the excavated mesh
for all  do
end for Compute the distance field between vertices of and
Step 3: Create Conforming DFM Mesh
if   then
end if If any vertex in the excavated mesh is too close to the DFN, then it is removed.
Create point cloud that is the union of the DFN vertices and remaining points in matrix mesh
Reconnect the point cloud using a Delaunay algorithm to create the final mesh.
Algorithm 2 Algorithm to generate a conforming Delaunay tetrahedralization of a three-dimensional fracture network.


  • (1) National Research Council, Rock fractures and fluid flow: contemporary understanding and applications, National Academy Press, 1996.
  • (2) S. Neuman, Trends, prospects and challenges in quantifying flow and transport through fractured rocks, Hydrogeol. J. 13 (1) (2005) 124–147.
  • (3) J. VanderKwaak, E. Sudicky, Dissolution of non-aqueous-phase liquids and aqueous-phase contaminant transport in discretely-fractured porous media, J. Contam. Hydrol. 23 (1-2) (1996) 45–68.
  • (4) B. H. Kueper, D. B. McWhorter, The behavior of dense, nonaqueous phase liquids in fractured clay and rock, Ground Water 29 (5) (1991) 716–728.
  • (5) J. Hyman, J. Jiménez-Martínez, H. Viswanathan, J. Carey, M. Porter, E. Rougier, S. Karra, Q. Kang, L. Frash, L. Chen, et al., Understanding hydraulic fracturing: a multi-scale problem, Phil. Trans. R. Soc. A 374 (2078) (2016) 20150426.
  • (6) R. Middleton, J. Carey, R. Currier, J. Hyman, Q. Kang, S. Karra, J. Jiménez-Martínez, M. Porter, H. Viswanathan, Shale gas and non-aqueous fracturing fluids: Opportunities and challenges for supercritical CO, Appl. Energ. 147 (2015) 500–509.
  • (7) T. Hadgu, S. Karra, E. Kalinina, N. Makedonska, J. D. Hyman, K. Klise, H. S. Viswanathan, Y. Wang, A comparative study of discrete fracture network and equivalent continuum models for simulating flow and transport in the far field of a hypothetical nuclear waste repository in crystalline host rock, Journal of Hydrology 553 (2017) 59 – 70. doi:https://doi.org/10.1016/j.jhydrol.2017.07.046.
  • (8) S. Joyce, L. Hartley, D. Applegate, J. Hoek, P. Jackson, Multi-scale groundwater flow modeling during temperate climate conditions for the safety assessment of the proposed high-level nuclear waste repository site at Forsmark, Sweden, Hydrogeol. J. 22 (6) (2014) 1233–1249.
  • (9) J. D. Hyman, J. Jiménez-Martínez, C. W. Gable, P. H. Stauffer, R. J. Pawar, Characterizing the impact of fractured caprock heterogeneity on supercritical co injection, Transport in Porous Media (2019) 1–21.
  • (10) C. Jenkins, A. Chadwick, S. D. Hovorka, The state of the art in monitoring and verification—ten years on, Int. J. Greenh. Gas. Con. 40 (2015) 312–349.
  • (11) I. Berre, F. Doster, E. Keilegavlen, Flow in fractured porous media: A review of conceptual models and discretization approaches, Transport in Porous Media (Oct 2018). doi:10.1007/s11242-018-1171-6.
  • (12) S. P. Neuman, J. S. Depner, Use of variable-scale pressure test data to estimate the log hydraulic conductivity covariance and dispersivity of fractured granites near oracle, arizona, J. Hydrol. 102 (1-4) (1988) 475–501.
  • (13) Y. Tsang, C. Tsang, F. Hale, B. Dverstorp, Tracer transport in a stochastic continuum model of fractured media, Water Resour. Res 32 (10) (1996) 3077–3092.
  • (14) H. Gerke, M. T. Van Genuchten, A dual-porosity model for simulating the preferential movement of water and solutes in structured porous media, Water Rescour. Res. 29 (2) (1993) 305–319.
  • (15) P. Lichtner, S. Karra, Modeling multiscale-multiphase-multicomponent reactive flows in porous media: Application to co2 sequestration and enhanced geothermal energy using PFLOTRAN, in: Al-Khoury, R., Bundschuh, J. (eds.) Computational Models for CO2 Geo-sequestration & Compressed Air Energy Storage (http://www.crcnetbase.com/doi/pdfplus/10), CRC Press, 2014, pp. 81–136.
  • (16) R. W. Zimmerman, G. Chen, T. Hadgu, G. S. Bodvarsson, A numerical dual-porosity model with semianalytical treatment of fracture/matrix flow, Water Resour. Res 29 (7) (1993) 2127–2137.
  • (17) M. C. Cacas, E. Ledoux, G. De Marsily, A. Barbreau, P. Calmels, B. Gaillard, R. Margritta, Modeling fracture flow with a stochastic discrete fracture network: Calibration and validation: 2. The transport model, Water Resour. Res. 26 (3) (1990) 491–500.
  • (18) J. Long, J. Remer, C. Wilson, P. Witherspoon, Porous media equivalents for networks of discontinuous fractures, Water Resour. Res 18 (3) (1982) 645–658.
  • (19) J. Maillot, P. Davy, R. Le Goc, C. Darcel, J.-R. De Dreuzy, Connectivity, permeability, and channeling in randomly distributed and kinematically defined discrete fracture network models, Water Resour. Res. 52 (11) (2016) 8526–8545.
  • (20) N. Makedonska, J. D. D. Hyman, S. Karra, S. L. Painter, C. W. W. Gable, H. S. Viswanathan, Evaluating the effect of internal aperture variability on transport in kilometer scale discrete fracture networks, Adv. Water Resour. 94 (2016) 486–497.
  • (21) A. W. Nordqvist, Y. W. Tsang, C. F. Tsang, B. Dverstorp, J. Andersson, A variable aperture fracture network model for flow and transport in fractured rocks, Water Resources Research 28 (6) (1992) 1703–1713. doi:10.1029/92WR00216.
    URL http://dx.doi.org/10.1029/92WR00216
  • (22) R. Ahmed, M. G. Edwards, S. Lamine, B. A. Huisman, M. Pal, Control-volume distributed multi-point flux approximation coupled with a lower-dimensional fracture model, Journal of Computational Physics 284 (2015) 462–489.
  • (23) R. Ahmed, M. G. Edwards, S. Lamine, B. A. Huisman, M. Pal, Three-dimensional control-volume distributed multi-point flux approximation coupled with a lower-dimensional surface fracture model, Journal of Computational Physics 303 (2015) 470–497.
  • (24) P. F. Antonietti, L. Formaggia, A. Scotti, M. Verani, N. Verzott, Mimetic finite difference approximation of flows in fractured porous media, ESAIM: Mathematical Modelling and Numerical Analysis 50 (3) (2016) 809–832.
  • (25) A. Arraras, F. J. Gaspar, L. Portero, C. Rodrigo, Mixed-dimensional geometric multigrid methods for single-phase flow in fractured porous media, SIAM Journal on Scientific Computing 41 (5) (2019) B1082–B1114.
  • (26) W. M. Boon, J. M. Nordbotten, I. Yotov, Robust discretization of flow in fractured porous media, SIAM Journal on Numerical Analysis 56 (4) (2018) 2203–2233.
  • (27) S. Berrone, A. Borio, C. Fidelibus, S. Pieraccini, S. Scialo, F. Vicini, Advanced computation of steady-state fluid flow in discrete fracture-matrix models: Fem–bem and vem–vem fracture-block coupling, GEM-International Journal on Geomathematics 9 (2) (2018) 377–399.
  • (28) B. Flemisch, A. Fumagalli, A. Scotti, A Review of the XFEM-Based Approximation of Flow in Fractured Porous Media, Springer International Publishing, Cham, 2016, pp. 47–76. doi:10.1007/978-3-319-41246-7-3.
  • (29) A. Fumagalli, A. Scotti, A reduced model for flow and transport in fractured porous media with non-matching grids, in: A. Cangiani, R. L. Davidchack, E. Georgoulis, A. N. Gorban, J. Levesley, M. V. Tretyakov (Eds.), Numerical Mathematics and Advanced Applications 2011, Springer Berlin Heidelberg, Berlin, Heidelberg, 2013, pp. 499–507.
  • (30) L. H. Odsæter, T. Kvamsdal, M. G. Larson, A simple embedded discrete fracture–matrix model for a coupled flow and transport problem in porous media, Computer Methods in Applied Mechanics and Engineering 343 (2019) 572 – 601. doi:https://doi.org/10.1016/j.cma.2018.09.003.
  • (31) M. Köppel, V. Martin, J. Jaffré, J. E. Roberts, A lagrange multiplier method for a discrete fracture model for flow in porous media, Computational Geosciences 23 (2) (2019) 239–253. doi:10.1007/s10596-018-9779-8.
  • (32) S. Manzoor, M. G. Edwards, A. H. Dogru, T. M. Al-Shaalan, Interior boundary-aligned unstructured grid generation and cell-centered versus vertex-centered cvd-mpfa performance, Computational Geosciences 22 (1) (2018) 195–230.
  • (33) T. H. Sandve, I. Berre, J. M. Nordbotten, An efficient multi-point flux approximation method for discrete fracture–matrix simulations, Journal of Computational Physics 231 (9) (2012) 3784–3800.
  • (34) N. Schwenck, B. Flemisch, R. Helmig, B. I. Wohlmuth, Dimensionally reduced flow models in fractured porous media: crossings and boundaries, Computational Geosciences 19 (6) (2015) 1219–1230. doi:10.1007/s10596-015-9536-1.
  • (35) Svensk Kärnbränslehantering AB, Data report for the safety assessment SR-site (TR-10-52), Tech. rep., Svensk Kärnbränslehantering AB (2010).
  • (36) L. Hartley, S. Joyce, Approaches and algorithms for groundwater flow modeling in support of site investigations and safety assessment of the forsmark site, sweden, Journal of Hydrology 500 (2013) 200 – 216. doi:https://doi.org/10.1016/j.jhydrol.2013.07.031.
  • (37) C. P. Jackson, A. R. Hoch, S. Todman, Self-consistency of a heterogeneous continuum porous medium representation of a fractured medium, Water Resources Research 36 (1) (2000) 189–202. doi:10.1029/1999WR900249.
  • (38) M. R. Sweeney, C. W. Gable, S. Karra, P. H. Stauffer, R. J. Pawar, J. D. Hyman, Upscaled discrete fracture matrix model (udfm): an octree-refined continuum representation of fractured porous media, Computational Geosciences (2019) 1–18.
  • (39) S. Berrone, S. Pieraccini, S. Scialo, A pde-constrained optimization formulation for discrete fracture network flows, SIAM J. Sci. Comput. 35 (2) (2013) B487–B510.
  • (40) S. Berrone, S. Pieraccini, S. Scialò, F. Vicini, A parallel solver for large scale dfn flow simulations, SIAM J. Sci. Comput. 37 (3) (2015) C285–C306.
  • (41) P. Davy, R. Le Goc, C. Darcel, A model of fracture nucleation, growth and arrest, and consequences for fracture density and scaling, J. Geophys. Res.-Sol. Ea. 118 (4) (2013) 1393–1407.
  • (42) P. Davy, R. Le Goc, C. Darcel, O. Bour, J. R. de Dreuzy, R. Munier, A likely universal model of fracture scaling and its consequence for crustal hydromechanics, Journal of Geophysical Research: Solid Earth 115 (B10) (2010). doi:10.1029/2009JB007043.
  • (43) J.-R. de Dreuzy, C. Darcel, P. Davy, O. Bour, Influence of spatial correlation of fracture centers on the permeability of two-dimensional fracture networks following a power law length distribution, Water Resour. Res. 40 (1) (2004).
  • (44) W. Dershowitz, C. Fidelibus, Derivation of equivalent pipe network analogues for three-dimensional discrete fracture networks by the boundary element method, Water Rescour. Res. 35 (9) (1999) 2685–2691.
  • (45) J.-R. de Dreuzy, Y. Méheust, G. Pichot, Influence of fracture scale heterogeneity on the flow properties of three-dimensional discrete fracture networks, J. Geophys. Res.-Sol. Ea. 117 (B11) (2012).
  • (46) J. Erhel, J.-R. de Dreuzy, B. Poirriez, Flow simulation in three-dimensional discrete fracture networks, SIAM J. Sci. Comput. 31 (4) (2009) 2688–2705.
  • (47) J. D. Hyman, C. W. Gable, S. L. Painter, N. Makedonska, Conforming Delaunay triangulation of stochastically generated three dimensional discrete fracture networks: A feature rejection algorithm for meshing strategy, SIAM J. Sci. Comput. 36 (4) (2014) A1871–A1894.
  • (48) J. D. Hyman, S. Karra, N. Makedonska, C. W. Gable, S. L. Painter, H. S. Viswanathan, dfnWorks: A discrete fracture network framework for modeling subsurface flow and transport, Comput. Geosci. 84 (2015) 10–19.
  • (49) G. Pichot, J. Erhel, J.-R. de Dreuzy, A generalized mixed hybrid mortar method for solving flow in stochastic discrete fracture networks, SIAM J. Sci. Comput. 34 (1) (2012) B86–B105.
  • (50) H. Mustapha, K. Mustapha, A new approach to simulating flow in discrete fracture networks with an optimized mesh, SIAM J. Sci. Comput. 29 (2007) 1439.
  • (51) E. Bonnet, O. Bour, N. E. Odling, P. Davy, I. Main, P. Cowie, B. Berkowitz, Scaling of fracture systems in geological media, Reviews of Geophysics 39 (3) (2001) 347–383.
  • (52) R. Ahmed, M. G. Edwards, S. Lamine, B. A. Huisman, M. Pal, Cvd-mpfa full pressure support, coupled unstructured discrete fracture–matrix darcy-flux approximations, Journal of Computational Physics 349 (2017) 265–299.
  • (53) I. Bogdanov, V. Mourzenko, J.-F. Thovert, P. Adler, Two-phase flow through fractured porous media, Physical Review E 68 (2) (2003) 026703.
  • (54) R. Helmig, Multiphase flow and transport processes in the subsurface: a contribution to the modeling of hydrosystems., Springer-Verlag, 1997.
  • (55) M. Karimi-Fard, L. J. Durlofsky, K. Aziz, et al., An efficient discrete fracture model applicable for general purpose reservoir simulators, in: SPE Reservoir Simulation Symposium, Society of Petroleum Engineers, 2003.
  • (56) S. Berrone, S. Pieraccini, S. Scialò, Flow simulations in porous media with immersed intersecting fractures, Journal of Computational Physics 345 (2017) 768–791.
  • (57) N. Frih, V. Martin, J. E. Roberts, A. Saâda, Modeling fractures as interfaces with nonmatching grids, Computational Geosciences 16 (4) (2012) 1043–1060.
  • (58) B. Flemisch, I. Berre, W. Boon, A. Fumagalli, N. Schwenck, A. Scotti, I. Stefansson, A. Tatomir, Benchmarks for single-phase flow in fractured porous media, Advances in Water Resources 111 (2018) 239–258.
  • (59) M. Del Pra, A. Fumagalli, A. Scotti, Well posedness of fully coupled fracture/bulk darcy flow with xfem, SIAM Journal on Numerical Analysis 55 (2) (2017) 785–811.
  • (60) J. M. Nordbotten, W. M. Boon, A. Fumagalli, E. Keilegavlen, Unified approach to discretization of flow in fractured porous media, Computational Geosciences 23 (2) (2019) 225–237.
  • (61) A. S. Abushaikha, K. M. Terekhov, A fully implicit mimetic finite difference scheme for general purpose subsurface reservoir simulation with full tensor permeability, Journal of Computational Physics 406 (2020) 109194.
  • (62) K. Lipnikov, G. Manzini, M. Shashkov, Mimetic finite difference method, J. Comput. Phys. 257 (2014) 1163–1227.
  • (63) V. Martin, J. Jaffré, J. E. Roberts, Modeling fractures and barriers as interfaces for flow in porous media, SIAM J. Sci. Comput 26 (2005) 1667–1691.
  • (64) D. S. Konstantin Lipnikov, Mikhail Shashkov, The mimetic finite difference discretization of diffusion problem on unstructured polyhedral meshes, Journal of Computational Physics 2 (211) (2006) 473–491.
  • (65) T. Barth, A 3-D Upwind Euler Solver For Unstructured Meshes, in: AIAA 10th Computational Fluid Dynamics Conference, 1991, pp. 228–238, 10th Conf On Computational Fluid Dynamics, Honolulu, HI, Jun 24-27, 1991.
  • (66) J. D. Moulton, J. C. Meza, M. W. Buksas, M. Day, et al., High-level design of Amanzi, the multi-process high performance computing simualtor, Tech. Rep. ASCEM-HPC-2011-03-1 (LANL: LA-UR 12-22193), Office of Environmental Management, United States Department of Energy, Washington, DC (2012).
  • (67) E. T. Coon, J. D. Moulton, S. L. Painter, Managing complexity in simulations of land surface and near-surface processes, Environmental Modelling & Software 78 (2016) 134–149. doi:http://dx.doi.org/10.1016/j.envsoft.2015.12.017.
  • (68) K. Lipnikov, D. Moulton, D. Svyatskiy, New preconditioning strategy for jacobian-free solvers for variably saturated flows with richards’ equation, Adv. Water Resour. 94 (2016) 11–22. doi:http://dx.doi.org/10.1016/j.advwatres.2016.04.016.
  • (69) I. Berre, W. M. Boon, B. Flemisch, A. Fumagalli, D. Gläser, E. Keilegavlen, A. Scotti, I. Stefansson, A. Tatomir, K. Brenner, et al., Verification benchmarks for single-phase flow in three-dimensional fractured porous media, Advances in Water Resources 147 (2021) 103759.
  • (70) M. Murphy, D. M. Mount, C. W. Gable, A point-placement strategy for conforming Delaunay tetrahedralization, International Journal of Computational Geometry & Applications 11 (06) (2001) 669–682.
  • (71) A. Khamayseh, A. Kuprat, Anisotropic smoothing and solution adaption for unstructured grids, International Journal for Numerical Methods in Engineering 39 (18) (1996) 3163–3174.
  • (72) R. E. Bank, R. K. Smith, Mesh smoothing using a posteriori error estimates, SIAM Journal on Numerical Analysis 34 (3) (1997) 979–997.
  • (73) LaGriT, Los Alamos Grid Toolbox, (LaGriT) Los Alamos National Laboratory, http://lagrit.lanl.gov (2013).
  • (74) A. Kuprat, Adaptive smoothing techniques for 3-d unstructured meshes (4 1996). doi:10.2172/226042.
    URL https://www.osti.gov/biblio/226042
  • (75) J. Guo, D.-M. Yan, L. Chen, X. Zhang, O. Deussen, P. Wonka, Tetrahedral meshing via maximal poisson-disk sampling, Computer Aided Geometric Design 43 (2016) 186–199.
  • (76) J. Krotz, M. R. Sweeney, C. W. Gable, J. D. Hyman, J. M. Restrepo, Maximal poisson-disk sampling for variable resolution conforming delaunay mesh generation: Applications for three-dimensional discrete fracture networks and the surrounding volume, arXiv preprint arXiv:2105.10079 (2021).
  • (77) H. Abelin, L. Birgersson, L. Moreno, H. Widén, T. Ågren, I. Neretnieks, A large-scale flow and tracer experiment in granite: 2. results and interpretation, Water Resour. Res. 27 (12) (1991) 3119–3135.
  • (78) H. Abelin, I. Neretnieks, S. Tunbrant, L. Moreno, Final report of the migration in a single fracture: experimental results and evaluation, Nat. Genossenschaft fd Lagerung Radioaktiver Abfälle, 1985.
  • (79) J. D. Hyman, Flow channeling in fracture networks: Characterizing the effect of density on preferential flow path formation, Water Resources Research (2020).
  • (80) A. Rasmuson, I. Neretnieks, Radionuclide transport in fast channels in crystalline rock, Water Resources Research 22 (8) (1986) 1247–1256.
  • (81) J. D. Hyman, M. Dentz, A. Hagberg, P. Kang, Linking structural and transport properties in three-dimensional fracture networks, J. Geophys. Res. Sol. Ea. (2019).
  • (82) T. Sherman, J. D. Hyman, M. Dentz, D. Bolster, Characterizing the influence of fracture density on network scale transport, J. Geophys. Res. Sol. Ea. 125 (1) (2020) e2019JB018547, e2019JB018547 10.1029/2019JB018547. arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2019JB018547, doi:10.1029/2019JB018547.
    URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019JB018547
  • (83) J.-R. De Dreuzy, P. Davy, O. Bour, Hydraulic properties of two-dimensional random fracture networks following power law distributions of length and aperture, Water Rescour. Res. 38 (12) (2002) 12–1.
  • (84) J. D. Hyman, G. Aldrich, H. Viswanathan, N. Makedonska, S. Karra, Fracture size and transmissivity correlations: Implications for transport simulations in sparse three-dimensional discrete fracture networks following a truncated power law distribution of fracture size, Water Resour. Res. 52 (8) (2016) 6472–6489. doi:10.1002/2016WR018806.
  • (85) P. Grindrod, M. Impey, Channeling and fickian dispersion in fractal simulated porous media, Water Resour. Res. 29 (12) (1993) 4077–4089.
  • (86) G. Aldrich, J. D. Hyman, S. Karra, C. W. Gable, N. Makedonska, H. Viswanathan, J. Woodring, B. Hamann, Analysis and visualization of discrete fracture networks using a flow topology graph, IEEE transactions on visualization and computer graphics 23 (8) (2017) 1896–1909.
  • (87) J. D. Hyman, A. Hagberg, G. Srinivasan, J. Mohd-Yusof, H. Viswanathan, Predictions of first passage times in sparse discrete fracture networks using graph-based reductions, Phys. Rev. E 96 (2017) 013304. doi:10.1103/PhysRevE.96.013304.
  • (88) J. D. Hyman, A. Hagberg, D. Osthus, S. Srinivasan, H. Viswanathan, G. Srinivasan, Identifying backbones in three-dimensional discrete fracture networks: A bipartite graph-based approach, Multiscale Modeling & Simulation 16 (4) (2018) 1948–1968.
  • (89) T. P. Wellman, A. M. Shapiro, M. C. Hill, Effects of simplifying fracture network representation on inert chemical migration in fracture-controlled aquifers, Water Resour. Res. 45 (1) (2009).
  • (90) A. T. Wood, Simulation of the von Mises Fisher distribution, Commun. Stat. Simulat. 23 (1) (1994) 157–164.