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 dualporosity / dualpermeability 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 meterlong 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 timescale 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 lengthscale, 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 threedimensional 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 nonconforming 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 mixeddimensional nature of the DFM, where fractures are codimension 1 of the matrix. Nonconforming 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, nonconforming methods have applied extended finite elements
del2017well ; Flemisch2016 ; fum ; Schwenck2015 , embedded finite elements ODS , finite elements/boundary elements berrone2018advanced , mixed finite element boon2018robust , mortartype 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 mixeddimensional 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 nonconforming mesh methods. The most commonly used conservative discretization methods for conforming unstructured meshes are controlvolume distributed multipoint flux approximations ahmed2015control ; ahmed2015three ; ahmed2017cvd ; manzoor2018interior ; sandve2012efficient , mixed finite elements, and recently mimetic finite difference antonietti2016mimetic . These methods are widely used in nonDFM applications and have seen extensive development, including the construction of higherorder methods for flow and transport and multiphase flow, which is promising for their eventual application to DFM models abushaikha2020fully ; lipnikov2014mimetic . However, the tradeoff 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 nonconforming methods. To date, most conforming DFM models in the literature are either twodimensional ahmed2015control ; ahmed2017cvd , quasi twodimensional ahmed2015three , or relatively simple threedimensional geometries berre2018flow due to the complexity of mesh generation.We present a conforming mesh methodology for resolving threedimensional 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 secondorder accurate mimetic finite differencing, and (3) implementation of numerical methods for highperformance 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 mixeddimensional 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, singlephase 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
are(1) 
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 nonreactive solute in is(2) 
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 twodimensional 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 welldefined 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):
(3) 
Within , the mass balance equation and the constitutive law are given by:
(4) 
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 righthand side denote the difference between values on two opposite sides of a fracture. Thus, the righthand 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 matrixfracture interface the mass conservation law leads to the velocity continuity condition:
(5) 
In the fracture, averaging the balance equation in normal direction leads to
(6) 
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 jaffrecoupled .
Similarly, the conceptual PDE model of the coupled transport of nonreactive solute is
(7) 
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 opensource 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) cellbased, denoted with a subscript in matrix and in fractures residing at cell centers (, ), and (2) facebased, denoted with a subscript residing at cell faces ( and ). In , there is one cellbased pressure unknown per volume cell and one facebased value per cell face. In , there is one cellbased pressure unknown per fracture cell, and one facebased 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 facebased 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).
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 saddlepoint forms:
(8) 
and
(9) 
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 :
(10) 
The local systems are coupled by the flux continuity conditions inside and and on :
(11) 
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 nearboundary mesh cells. Recall, that denote facebased degrees of freedom on the opposite sides of fracture cell . The last formula shows the convenience of using facebased 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 saddlepoint form to symmetric positive definite (SPD) formulations:
(12) 
and
(13) 
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:
(14) 
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 cellcentered degrees of freedom for volumetric matrix cells and surface fracture cells. The advective solute fluxes are approximated using the firstorder 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:
(15) 
where
(16) 
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 righthand side of equation (15) as by analogy with the flux definition.
Using the firstorder time approximation, we replace equations (7) with discrete equations
(17) 
Setting , we obtain the explicit scheme. The implicit scheme is obtained by setting . A secondorder 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 secondorder accurate time integration is implemented via the RungeKutta scheme.
2.3 Implementation: Process Kernels
We implemented the coupled matrixfracture flow system in the opensource code Amanzi ASCEM_Amanzi_design_1 which provides a suite of critical lowlevel 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 subdomains, 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 thirdparty solvers such as Hypre’s algebraic multigrid.
Arcos, the underlying framework of Amanzi, provides services to manage processrich 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 multiprocess 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 PKtree to represent the coupled system. A PKtree for the coupled matrixfracture flow and transport system is shown in Figure
2. The tight coupling of the matrixfracture 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 matrixfracture 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 advectiondiffusion system formed by combining the matrices from the matrix and fracture systems. The toplevel 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 subcycled relative to the flow MPC.At a lower level, Arcos dynamically manages a complete and selfconsistent 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.
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 twodimensional mesh of is the trace of the threedimensional 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 (FRAM) hyman2014conforming 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 userdefined 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 psuedocode 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.

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 .

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 threedimensional 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 threedimensional 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. 
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.

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 threedimensional 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.
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 fracturematrix 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 threedimensional extension of FRAM and the second is a graphbased/upscaling method that reduces the network complexity while preserving key hydrological properties of the fracture and matrix system.
3.3 A ThreeDimensional 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 threedimensional space that affect the ability to automatically generate a conforming volume mesh of the surrounding matrix. Examples of threedimensional 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 threedimensions 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 threedimensional 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 subnetworks, 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 macroscale structure of the network, e.g., orientations and density hyman2019linking ; sherman2020characterizing ; hyman2020flow , mesoscale 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) 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 powerlaw 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 vonMisser 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 deadend fractures, i.e., isolating the 2core 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.
The principal parts of the methods are described in Algorithm 1.
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 singlephase flow and transport problems in DFM systems in berre2021verification . In that study, four benchmark cases for singlephase flow and transport in threedimensional 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 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 10meter 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.
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 orderaccuracy 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 xdirection 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.18e6 
0.9  1364  59111  1.04e6 
0.8  1695  88136  8.56e7 
0.7  2183  124345  6.80e7 
0.6  2970  202275  5.35e7 
0.5  4266  369361  4.12e7 
0.4  6674  716722  2.69e7 
0.3  11738  1686799  8.99e8 
0.25  16965  2991156   
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 xaxis. 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 2dashed, case 3dotted), 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.
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 multidimensional 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 nonreactive tracer.
Fig. 10 a,b shows snapshots of the the steadystate 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 singlephase 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 threedimensional length scales are larger than a userdefined lengthscale. 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 secondorder accurate mimetic finite differencing. The implementation of these numerical methods for highperformance computing environments is performed using a hierarchically designed multiprocess coordinator and process kernel structure that allows for both strong and weak coupling of the fracture/matrix domains. Additional physical models beyond simulate singlephase flow and transport supported by the adopted opensource 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 threedimensional stochastically generated network containing 257 fractures are presented. The methods presented here will be made freely available in the opensource code dfnWorks in a future release.
Acknowledgments
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 # LAUR2131458.
Appendix A Meshing Algorithm
References
 (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 nonaqueousphase liquids and aqueousphase contaminant transport in discretelyfractured porous media, J. Contam. Hydrol. 23 (12) (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énezMartínez, H. Viswanathan, J. Carey, M. Porter, E. Rougier, S. Karra, Q. Kang, L. Frash, L. Chen, et al., Understanding hydraulic fracturing: a multiscale 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énezMartínez, M. Porter, H. Viswanathan, Shale gas and nonaqueous 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, Multiscale groundwater flow modeling during temperate climate conditions for the safety assessment of the proposed highlevel nuclear waste repository site at Forsmark, Sweden, Hydrogeol. J. 22 (6) (2014) 1233–1249.
 (9) J. D. Hyman, J. JiménezMartí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/s1124201811716.
 (12) S. P. Neuman, J. S. Depner, Use of variablescale pressure test data to estimate the log hydraulic conductivity covariance and dispersivity of fractured granites near oracle, arizona, J. Hydrol. 102 (14) (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 dualporosity 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 multiscalemultiphasemulticomponent reactive flows in porous media: Application to co2 sequestration and enhanced geothermal energy using PFLOTRAN, in: AlKhoury, R., Bundschuh, J. (eds.) Computational Models for CO2 Geosequestration & 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 dualporosity 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, Controlvolume distributed multipoint flux approximation coupled with a lowerdimensional fracture model, Journal of Computational Physics 284 (2015) 462–489.
 (23) R. Ahmed, M. G. Edwards, S. Lamine, B. A. Huisman, M. Pal, Threedimensional controlvolume distributed multipoint flux approximation coupled with a lowerdimensional 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, Mixeddimensional geometric multigrid methods for singlephase 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 steadystate fluid flow in discrete fracturematrix models: Fem–bem and vem–vem fractureblock coupling, GEMInternational Journal on Geomathematics 9 (2) (2018) 377–399.
 (28) B. Flemisch, A. Fumagalli, A. Scotti, A Review of the XFEMBased Approximation of Flow in Fractured Porous Media, Springer International Publishing, Cham, 2016, pp. 47–76. doi:10.1007/97833194124673.
 (29) A. Fumagalli, A. Scotti, A reduced model for flow and transport in fractured porous media with nonmatching 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/s1059601897798.
 (32) S. Manzoor, M. G. Edwards, A. H. Dogru, T. M. AlShaalan, Interior boundaryaligned unstructured grid generation and cellcentered versus vertexcentered cvdmpfa performance, Computational Geosciences 22 (1) (2018) 195–230.
 (33) T. H. Sandve, I. Berre, J. M. Nordbotten, An efficient multipoint 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/s1059601595361.
 (35) Svensk Kärnbränslehantering AB, Data report for the safety assessment SRsite (TR1052), 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, Selfconsistency 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 octreerefined continuum representation of fractured porous media, Computational Geosciences (2019) 1–18.
 (39) S. Berrone, S. Pieraccini, S. Scialo, A pdeconstrained 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 twodimensional 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 threedimensional 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 threedimensional discrete fracture networks, J. Geophys. Res.Sol. Ea. 117 (B11) (2012).
 (46) J. Erhel, J.R. de Dreuzy, B. Poirriez, Flow simulation in threedimensional 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, Cvdmpfa full pressure support, coupled unstructured discrete fracture–matrix darcyflux approximations, Journal of Computational Physics 349 (2017) 265–299.
 (53) I. Bogdanov, V. Mourzenko, J.F. Thovert, P. Adler, Twophase 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., SpringerVerlag, 1997.
 (55) M. KarimiFard, 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 singlephase 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 3D 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 2427, 1991.
 (66) J. D. Moulton, J. C. Meza, M. W. Buksas, M. Day, et al., Highlevel design of Amanzi, the multiprocess high performance computing simualtor, Tech. Rep. ASCEMHPC2011031 (LANL: LAUR 1222193), 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 nearsurface 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 jacobianfree 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 singlephase flow in threedimensional fractured porous media, Advances in Water Resources 147 (2021) 103759.
 (70) M. Murphy, D. M. Mount, C. W. Gable, A pointplacement 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 3d 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 poissondisk 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 poissondisk sampling for variable resolution conforming delaunay mesh generation: Applications for threedimensional 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 largescale 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 threedimensional 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 twodimensional 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 threedimensional 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. MohdYusof, H. Viswanathan, Predictions of first passage times in sparse discrete fracture networks using graphbased 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 threedimensional discrete fracture networks: A bipartite graphbased 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 fracturecontrolled 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.
Comments
There are no comments yet.